<p><b>dwj07@fsu.edu</b> 2013-03-12 15:20:35 -0600 (Tue, 12 Mar 2013)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Updating shared advection routines.<br>
Creating a non-monotonic option (which might not be as efficient as it could be).<br>
<br>
Removing dead code from ocean core.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_operators/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/Makefile 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/Makefile 2013-03-12 21:20:35 UTC (rev 2593)
@@ -38,14 +38,6 @@
mpas_ocn_tendency.o \
mpas_ocn_diagnostics.o \
mpas_ocn_tracer_advection.o \
- mpas_ocn_tracer_advection_std.o \
- mpas_ocn_tracer_advection_std_hadv.o \
- mpas_ocn_tracer_advection_std_vadv.o \
- mpas_ocn_tracer_advection_std_vadv2.o \
- mpas_ocn_tracer_advection_std_vadv3.o \
- mpas_ocn_tracer_advection_std_vadv4.o \
- mpas_ocn_tracer_advection_mono.o \
- mpas_ocn_tracer_advection_helpers.o \
mpas_ocn_time_integration.o \
mpas_ocn_time_integration_rk4.o \
mpas_ocn_time_integration_split.o \
@@ -134,24 +126,8 @@
mpas_ocn_tracer_hmix_del4.o:
-mpas_ocn_tracer_advection.o: mpas_ocn_tracer_advection_std.o mpas_ocn_tracer_advection_mono.o
+mpas_ocn_tracer_advection.o:
-mpas_ocn_tracer_advection_std.o: mpas_ocn_tracer_advection_std_hadv.o mpas_ocn_tracer_advection_std_vadv.o
-
-mpas_ocn_tracer_advection_std_hadv.o: mpas_ocn_tracer_advection_helpers.o
-
-mpas_ocn_tracer_advection_std_vadv.o: mpas_ocn_tracer_advection_std_vadv2.o mpas_ocn_tracer_advection_std_vadv3.o mpas_ocn_tracer_advection_std_vadv4.o
-
-mpas_ocn_tracer_advection_std_vadv2.o: mpas_ocn_tracer_advection_helpers.o
-
-mpas_ocn_tracer_advection_std_vadv3.o: mpas_ocn_tracer_advection_helpers.o
-
-mpas_ocn_tracer_advection_std_vadv4.o: mpas_ocn_tracer_advection_helpers.o
-
-mpas_ocn_tracer_advection_mono.o: mpas_ocn_tracer_advection_helpers.o
-
-mpas_ocn_tracer_advection_helpers.o:
-
mpas_ocn_restoring.o:
mpas_ocn_vmix.o: mpas_ocn_vmix_coefs_const.o mpas_ocn_vmix_coefs_rich.o mpas_ocn_vmix_coefs_tanh.o
@@ -204,14 +180,6 @@
mpas_ocn_vmix_coefs_tanh.o \
mpas_ocn_restoring.o \
mpas_ocn_tracer_advection.o \
- mpas_ocn_tracer_advection_std.o \
- mpas_ocn_tracer_advection_std_hadv.o \
- mpas_ocn_tracer_advection_std_vadv.o \
- mpas_ocn_tracer_advection_std_vadv2.o \
- mpas_ocn_tracer_advection_std_vadv3.o \
- mpas_ocn_tracer_advection_std_vadv4.o \
- mpas_ocn_tracer_advection_mono.o \
- mpas_ocn_tracer_advection_helpers.o \
mpas_ocn_tendency.o \
mpas_ocn_diagnostics.o \
mpas_ocn_time_integration.o \
Modified: branches/ocean_projects/advection_operators/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/Registry 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/Registry 2013-03-12 21:20:35 UTC (rev 2593)
@@ -27,6 +27,7 @@
namelist character grid config_alter_ICs_for_pbcs 'zlevel_pbcs_off'
namelist real grid config_min_pbc_fraction 0.10
namelist logical grid config_check_ssh_consistency .true.
+namelist logical grid config_dzdk_positive .true.
namelist character decomposition config_block_decomp_file_prefix 'graph.info.part.'
namelist integer decomposition config_number_of_blocks 0
Modified: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -256,7 +256,7 @@
use mpas_grid_types
use mpas_rbf_interpolation
use mpas_vector_reconstruction
- use mpas_ocn_tracer_advection
+ use mpas_tracer_advection_helpers
use ocn_advection
implicit none
@@ -270,7 +270,7 @@
call ocn_setup_sign_and_index_fields(mesh)
call ocn_initialize_advection_rk(mesh, err)
- call mpas_ocn_tracer_advection_coefficients(mesh, err1)
+ call mpas_tracer_advection_coefficients(mesh, err1, mesh % maxLevelCell % array, mesh % highOrderAdvectionMask % array, mesh % boundaryCell % array)
err = ior(err, err1)
call ocn_time_average_init(block % state % time_levs(1) % state)
Modified: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -21,7 +21,7 @@
use mpas_sort
use mpas_hash
- use mpas_ocn_tracer_advection_std
+ use mpas_tracer_advection_std
use mpas_tracer_advection_mono
implicit none
@@ -29,7 +29,6 @@
save
public :: mpas_ocn_tracer_advection_init, &
- mpas_ocn_tracer_advection_coefficients, &
mpas_ocn_tracer_advection_tend
logical :: tracerAdvOn
@@ -39,196 +38,6 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
-! routine mpas_ocn_tracer_advection_coefficients
-!
-!> \brief MPAS ocean tracer advection coefficients
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine precomputes the advection coefficients for horizontal
-!> advection of tracers.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_coefficients( grid, err )!{{{
-
- implicit none
- type (mesh_type) :: grid !< Input: Grid information
- integer, intent(out) :: err
-
- 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, nCells
- logical :: addcell, highOrderAdvection
-
- type (hashtable) :: cell_hash
-
- deriv_two => grid % deriv_two % array
- adv_coefs => grid % adv_coefs % array
- adv_coefs_2nd => grid % adv_coefs_2nd % array
- adv_coefs_3rd => grid % adv_coefs_3rd % array
- cellsOnCell => grid % cellsOnCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- advCellsForEdge => grid % advCellsForEdge % array
- boundaryCell => grid % boundaryCell % array
- highOrderAdvectionMask => grid % highOrderAdvectionMask % array
- lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- maxLevelCell => grid % maxLevelCell % array
- nAdvCellsForEdge => grid % nAdvCellsForEdge % array
-
- nCells = grid % nCells
- nVertLevels = grid % nVertLevels
-
- allocate(cell_indices(grid % maxEdges2 + 2))
- allocate(sorted_cell_indices(2, grid % maxEdges2 + 2))
-
- err = 0
-
- highOrderAdvectionMask = 0
- lowOrderAdvectionMask = 0
-
- do iEdge = 1, grid % nEdges
- nAdvCellsForEdge(iEdge) = 0
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k = 1, nVertLevels
- if (boundaryCell(k, cell1) == 1 .or. boundaryCell(k, cell2) == 1) then
- highOrderAdvectionMask(k, iEdge) = 0
- lowOrderAdvectionMask(k, iEdge) = 1
- else
- highOrderAdvectionMask(k, iEdge) = 1
- lowOrderAdvectionMask(k, iEdge) = 0
- end if
- end do
-
- !
- ! do only if this edge flux is needed to update owned cells
- !
- if (cell1 <= grid % nCells .and. cell2 <= 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 <= 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 <= 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 <= 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 <= 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 <= 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 <= 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) then
- lowOrderAdvectionMask = 1
- highOrderAdvectionMask = 0
- end if
-
- if (maxval(highOrderAdvectionMask+lowOrderAdvectionMask) > 1) then
- write(*,*) "Masks don't sum to 1."
- err = 1
- endif
-
- end subroutine mpas_ocn_tracer_advection_coefficients!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
! routine mpas_ocn_tracer_advection_tend
!
!> \brief MPAS ocean tracer advection tendency
@@ -255,9 +64,9 @@
if(.not. tracerAdvOn) return
if(monotonicOn) then
- call mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)
+ call mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, grid, tend_h, tend, dt, grid % maxLevelCell % array, grid % maxLevelEdgeTop % array, grid % highOrderAdvectionMask % array)
else
- call mpas_ocn_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)
+ call mpas_tracer_advection_std_tend(tracers, uh, w, h, verticalCellSize, grid, tend_h, tend, dt, grid % maxLevelCell % array, grid % maxLevelEdgeTop % array, grid % highOrderAdvectionMask % array)
endif
end subroutine mpas_ocn_tracer_advection_tend!}}}
@@ -286,7 +95,7 @@
if(config_disable_tr_adv) tracerAdvOn = .false.
- call mpas_ocn_tracer_advection_std_init(err_tmp)
+ call mpas_tracer_advection_std_init(err_tmp)
call mpas_tracer_advection_mono_init(err_tmp)
err = ior(err, err_tmp)
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_helpers.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_helpers.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,68 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_helpers
-!
-!> \brief MPAS ocean tracer advection helper functions
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains helper functions tracer advection.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_helpers
-
- use mpas_kind_types
-
- implicit none
- save
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! function mpas_ocn_tracer_advection_vflux4
-!
-!> \brief MPAS ocean 4th order vertical tracer advection stencil
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This function provides the stencil for 4th order vertical advection of tracers.
-!
-!-----------------------------------------------------------------------
- real function mpas_ocn_tracer_advection_vflux4(q_im2, q_im1, q_i, q_ip1, w)!{{{
- real (kind=RKIND), intent(in) :: q_im2 !< Input: Tracer value at index i-2
- real (kind=RKIND), intent(in) :: q_im1 !< Input: Tracer value at index i-1
- real (kind=RKIND), intent(in) :: q_i !< Input: Tracer value at index i
- real (kind=RKIND), intent(in) :: q_ip1 !< Input: Tracer value at index i+1
- real (kind=RKIND), intent(in) :: w !< Input: vertical veloicity
- mpas_ocn_tracer_advection_vflux4 = w*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
- end function!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! function mpas_ocn_tracer_advection_vflux3
-!
-!> \brief MPAS ocean 3rd order vertical tracer advection stencil
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This function provides the stencil for 3rd order vertical advection of tracers.
-!
-!-----------------------------------------------------------------------
- real function mpas_ocn_tracer_advection_vflux3( q_im2, q_im1, q_i, q_ip1, w, coef)!{{{
- real (kind=RKIND), intent(in) :: q_im2 !< Input: Tracer value at index i-2
- real (kind=RKIND), intent(in) :: q_im1 !< Input: Tracer value at index i-1
- real (kind=RKIND), intent(in) :: q_i !< Input: Tracer value at index i
- real (kind=RKIND), intent(in) :: q_ip1 !< Input: Tracer value at index i+1
- real (kind=RKIND), intent(in) :: w !< Input: vertical veloicity
- real (kind=RKIND), intent(in) :: coef !< Input: Advection coefficient
-
- !dwj 02/21/12 flux3 is different in ocean and atmosphere
- !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
- mpas_ocn_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 function!}}}
-
-end module mpas_ocn_tracer_advection_helpers
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,401 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_mono
-!
-!> \brief MPAS ocean monotonic tracer advection with FCT
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains routines for monotonic advection of tracers using a FCT
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_mono
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_helpers
-
- implicit none
- private
- save
-
- real (kind=RKIND) :: coef_3rd_order
-
- public :: mpas_ocn_tracer_advection_mono_tend, &
- mpas_ocn_tracer_advection_mono_init
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_mono_tend
-!
-!> \brief MPAS ocean monotonic tracer advection tendency with FCT
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine computes the monotonic tracer advection tendencity using a FCT.
-!> Both horizontal and vertical.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, 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) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
- real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Tendency for thickness field
- real (kind=RKIND), intent(in) :: dt !< Input: Timestep
- type (mesh_type), intent(in) :: grid !< Input: Grid information
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
-
- 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, edgeSignOnCell, edgesOnCell
-
- 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
-
- ! 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_2nd => grid % adv_coefs_2nd % array
- adv_coefs_3rd => grid % adv_coefs_3rd % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- highOrderAdvectionMask => grid % highOrderAdvectionMask % array
- lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
- edgeSignOnCell => grid % edgeSignOnCell % array
- edgesOnCell => grid % edgesOnCell % array
-
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers,dim=1)
-
- allocate(high_order_horiz_flux_field)
- nullify(high_order_horiz_flux_field % next)
- high_order_horiz_flux_field % block => grid % block
- high_order_horiz_flux_field % sendList => grid % xEdge % sendList
- high_order_horiz_flux_field % recvList => grid % xEdge % recvList
- high_order_horiz_flux_field % copyList => 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 => 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 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_new is supposed to be the "new" tracer state. This allows bounds checks.
- if (config_check_tracer_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 = 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
- high_order_vert_flux(k,iCell) = mpas_ocn_tracer_advection_vflux3( 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)
- 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
- do i = 1, nAdvCellsForEdge(iEdge)
- iCell = advCellsForEdge(i,iEdge)
- do k = 1, maxLevelCell(iCell)
- tracer_weight = lowOrderAdvectionMask(k, iEdge) * adv_coefs_2nd(i,iEdge) &
- + highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,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
-
- ! 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 Ocean and Atmosphere are different in vertical
-! flux_upwind = max(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
- flux_upwind = min(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.0_RKIND,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.0_RKIND,high_order_vert_flux(k+1,iCell))-max(0.0_RKIND,high_order_vert_flux(k,iCell)))
-! flux_outgoing(k,iCell) = -(max(0.0_RKIND,high_order_vert_flux(k+1,iCell))-min(0.0_RKIND,high_order_vert_flux(k,iCell)))
-
- flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, 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.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell2))
- high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
- end do ! k loop
- end do ! iEdge loop
-
- do iCell = 1, nCells
- invAreaCell1 = 1.0 / areaCell(iCell)
- do i = 1, nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i, iCell)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1, maxLevelEdgeTop(iEdge)
- flux_upwind = dvEdge(iEdge) * (max(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell2))
-
- upwind_tendency(k,iCell) = upwind_tendency(k,iCell) + edgeSignOncell(i, iCell) * flux_upwind * invAreaCell1
-
- ! Accumulate remaining high order fluxes
- flux_outgoing(k,iCell) = flux_outgoing(k,iCell) + min(0.0_RKIND, edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge)) * invAreaCell1
- flux_incoming(k,iCell) = flux_incoming(k,iCell) + max(0.0_RKIND, edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge)) * invAreaCell1
- end do
- end do
- end do
-
- ! 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_RKIND, max( 0.0_RKIND, 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_RKIND, max( 0.0_RKIND, 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.0_RKIND,flux) * min(flux_outgoing(k,cell1), flux_incoming(k,cell2)) &
- + min(0.0_RKIND,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.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
-! + min(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
- flux = max(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
- + min(0.0_RKIND,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 iCell = 1, nCells
- invAreaCell1 = 1.0 / areaCell(iCell)
- do i = 1, nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i, iCell)
- do k = 1, maxLevelEdgeTop(iEdge)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
-
- if(config_check_tracer_monotonicity) then
- tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
- end if
- end do
- end do
- end do
-
- ! 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_tracer_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_tracer_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) < 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) > 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)
- end subroutine mpas_ocn_tracer_advection_mono_tend!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_mono_init
-!
-!> \brief MPAS ocean initialize monotonic tracer advection tendency with FCT
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine initializes the monotonic tracer advection tendencity using a FCT.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_mono_init(err)!{{{
- integer, intent(inout) :: err !< Input: Error Flags
-
- integer :: err_tmp
-
- err = 0
-
- 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
-
- end subroutine mpas_ocn_tracer_advection_mono_init!}}}
-
-end module mpas_ocn_tracer_advection_mono
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,100 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std
-!
-!> \brief MPAS ocean tracer advection driver (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains driver routine for tracer advection tendencies
-!> as well as the routines for setting up advection coefficients and
-!> initialization of the advection routines.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
- use mpas_timer
-
- use mpas_ocn_tracer_advection_std_hadv
- use mpas_ocn_tracer_advection_std_vadv
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_tend, &
- mpas_ocn_tracer_advection_std_init
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_tend
-!
-!> \brief MPAS ocean standard tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine is the driver routine for the standard computation of
-!> tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tracer tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer values
- real (kind=RKIND), dimension(:,:), intent(in) :: uh !< Input: Thickness weighted horizontal velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical Velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- call mpas_timer_start("tracer-hadv", .false.)
- call mpas_ocn_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
- call mpas_timer_stop("tracer-hadv")
- call mpas_timer_start("tracer-vadv", .false.)
- call mpas_ocn_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)
- call mpas_timer_stop("tracer-vadv")
-
- end subroutine mpas_ocn_tracer_advection_std_tend!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_init
-!
-!> \brief MPAS ocean standard tracer advection initialization
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine is the driver routine for the initializtion of the standard
-!> tracer advection routines.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_init(err)!{{{
- integer, intent(inout) :: err !< Input: Error Flags
-
- integer :: err_tmp
-
- err = 0
-
- call mpas_ocn_tracer_advection_std_hadv_init(err_tmp)
- err = ior(err, err_tmp)
- call mpas_ocn_tracer_advection_std_vadv_init(err_tmp)
- err = ior(err, err_tmp)
-
- end subroutine mpas_ocn_tracer_advection_std_init!}}}
-
-end module mpas_ocn_tracer_advection_std
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,140 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std_hadv
-!
-!> \brief MPAS ocean standard horizontal tracer advection (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains routines for horizontal tracer advection tendencies
-!> and initialization of the horizontal advection routines.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std_hadv
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_helpers
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_hadv_tend, &
- mpas_ocn_tracer_advection_std_hadv_init
-
- real (kind=RKIND) :: coef_3rd_order
-
- logical :: hadvOn
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_hadv_tend
-!
-!> \brief MPAS ocean standard horizontal tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine computes the tendency for 3rd order horizontal advection of tracers.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)!{{{
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/output: Tracer tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer values
- real (kind=RKIND), dimension(:,:), intent(in) :: uh !< Input: Thickness weighted horizontal velocity
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- 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, highOrderAdvectionMask, lowOrderAdvectionMask
- integer, dimension(:), pointer :: nAdvCellsForEdge
- real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
- 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_2nd => grid % adv_coefs_2nd % array
- adv_coefs_3rd => grid % adv_coefs_3rd % array
- highOrderAdvectionMask => grid % highOrderAdvectionMask % array
- lowOrderAdvectionMask => grid % lowOrderAdvectionMask % 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 = lowOrderAdvectionMask(k, iEdge) * adv_coefs_2nd(i,iEdge) &
- + highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,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_ocn_tracer_advection_std_hadv_tend!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_hadv_init
-!
-!> \brief MPAS ocean standard horizontal tracer advection initialization
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine initializes the 3rd order standard horizontal advection of tracers
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_hadv_init(err)!{{{
- integer, intent(inout) :: err !< Input/Output: Error flag
-
- err = 0
-
- hadvOn =.true.
-
- 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
- end subroutine mpas_ocn_tracer_advection_std_hadv_init!}}}
-
-end module mpas_ocn_tracer_advection_std_hadv
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,103 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std_vadv
-!
-!> \brief MPAS ocean vertical tracer advection driver (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains driver routines for vertical tracer advection tendencies
-!> and initialization of the advection routines.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std_vadv
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_std_vadv2
- use mpas_ocn_tracer_advection_std_vadv3
- use mpas_ocn_tracer_advection_std_vadv4
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_vadv_tend, &
- mpas_ocn_tracer_advection_std_vadv_init
-
- logical :: order2, order3, order4
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_vadv_tend
-!
-!> \brief MPAS ocean standard vertical tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine is the driver routine for the standard computation of
-!> vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)!{{{
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer Tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer Values
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical Velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of cell
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- if(order2) then
- call mpas_ocn_tracer_advection_std_vadv2_tend(tracers, w, verticalCellSize, grid, tend)
- else if(order3) then
- call mpas_ocn_tracer_advection_std_vadv3_tend(tracers, w, verticalCellSize, grid, tend)
- else if(order4) then
- call mpas_ocn_tracer_advection_std_vadv4_tend(tracers, w, verticalCellSize, grid, tend)
- endif
-
- end subroutine mpas_ocn_tracer_advection_std_vadv_tend!}}}
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_vadv_init
-!
-!> \brief MPAS ocean standard vertical tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine initializes the vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_vadv_init(err)!{{{
- integer, intent(inout) :: err !< Input/Output: Error flag
-
- err = 0
-
- order2 = .false.
- order3 = .false.
- order4 = .false.
-
- if (config_vert_tracer_adv_order == 2) then
- order2 = .true.
- else if (config_vert_tracer_adv_order == 3) then
- order3 = .true.
- else if (config_vert_tracer_adv_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_ocn_tracer_advection_std_vadv_init!}}}
-
-end module mpas_ocn_tracer_advection_std_vadv
-
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,96 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std_vadv2
-!
-!> \brief MPAS ocean 2nd order vertical tracer advection driver (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains routines for 2nd order vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std_vadv2
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_helpers
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_vadv2_tend
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_vadv2_tend
-!
-!> \brief MPAS ocean 2nd order standard vertical tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine is the driver routine for the 2nd order standard computation of
-!> vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_vadv2_tend(tracers, w, verticalCellSize, grid, tend)!{{{
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: tracer tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer values
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical Velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of cell
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- 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(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- 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
- do k = 2, maxLevelCell(iCell)
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*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))
- end do
- end do
- end do
-
- deallocate(vert_flux)
-
- end subroutine mpas_ocn_tracer_advection_std_vadv2_tend!}}}
-
-end module mpas_ocn_tracer_advection_std_vadv2
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,106 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std_vadv3
-!
-!> \brief MPAS ocean 3rd order vertical tracer advection driver (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains routines for 3rd order vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std_vadv3
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_helpers
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_vadv3_tend
-
- real (kind=RKIND) :: coef_3rd_order
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_vadv3_tend
-!
-!> \brief MPAS ocean 3rd order standard vertical tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine computes the 3rd order vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_vadv3_tend(tracers, w, verticalCellSize, grid, tend)!{{{
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer Tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer Values
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical Velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of cell
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-
- real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- integer :: nVertLevels, num_tracers
- integer, dimension(:), pointer :: maxLevelCell
-
- 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.
-
- do iCell=1,grid % nCellsSolve
-
- k = 2
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- enddo
-
- do k=3,maxLevelCell(iCell)-1
- do iTracer=1,num_tracers
- vert_flux(iTracer,k) = mpas_ocn_tracer_advection_vflux3( 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
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*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))
- end do
- end do
- end do
-
- deallocate(vert_flux)
-
- end subroutine mpas_ocn_tracer_advection_std_vadv3_tend!}}}
-
-end module mpas_ocn_tracer_advection_std_vadv3
Deleted: branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -1,105 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_tracer_advection_std_vadv4
-!
-!> \brief MPAS ocean 4th order vertical tracer advection driver (non-monotonic/fct)
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This module contains routines for 4th order vertical tracer advection.
-!
-!-----------------------------------------------------------------------
-module mpas_ocn_tracer_advection_std_vadv4
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
-
- use mpas_ocn_tracer_advection_helpers
-
- implicit none
- private
- save
-
- public :: mpas_ocn_tracer_advection_std_vadv4_tend
-
- contains
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! routine mpas_ocn_tracer_advection_std_vadv4_tend
-!
-!> \brief MPAS ocean 4th order standard vertical tracer advection tendency
-!> \author Doug Jacobsen
-!> \date 03/09/12
-!> \version SVN:$Id:$
-!> \details
-!> This routine computes the 4th order vertical tracer advection tendencies.
-!
-!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_std_vadv4_tend(tracers, w, verticalCellSize, grid, tend)!{{{
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: Tracer Values
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical Velocity
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of cell
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-
- real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- 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
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- enddo
- do k=3,nVertLevels-1
- do iTracer=1,num_tracers
- vert_flux(iTracer,k) = mpas_ocn_tracer_advection_vflux4( 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
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*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))
- end do
- end do
-
- end do
-
- deallocate(vert_flux)
-
- end subroutine mpas_ocn_tracer_advection_std_vadv4_tend!}}}
-
-end module mpas_ocn_tracer_advection_std_vadv4
Modified: branches/ocean_projects/advection_operators/src/operators/Makefile
===================================================================
--- branches/ocean_projects/advection_operators/src/operators/Makefile 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/operators/Makefile 2013-03-12 21:20:35 UTC (rev 2593)
@@ -4,6 +4,7 @@
mpas_vector_reconstruction.o \
mpas_spline_interpolation.o \
mpas_tracer_advection_mono.o \
+ mpas_tracer_advection_std.o \
mpas_tracer_advection_helpers.o
all: operators
@@ -12,6 +13,7 @@
ar -ru libops.a $(OBJS)
mpas_tracer_advection_mono.o: mpas_tracer_advection_helpers.o
+mpas_tracer_advection_std.o: mpas_tracer_advection_helpers.o
mpas_tracer_advection_helpers.o:
mpas_vector_reconstruction.o: mpas_rbf_interpolation.o
mpas_rbf_interpolation.o:
Modified: branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_helpers.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_helpers.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -12,7 +12,9 @@
!-----------------------------------------------------------------------
module mpas_tracer_advection_helpers
+ use mpas_grid_types
use mpas_kind_types
+ use mpas_configure
implicit none
save
@@ -21,6 +23,204 @@
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
+! routine mpas_tracer_advection_coefficients
+!
+!> \brief MPAS tracer advection coefficients
+!> \author Doug Jacobsen
+!> \date 03/09/12
+!> \version SVN:$Id:$
+!> \details
+!> This routine precomputes the advection coefficients for horizontal
+!> advection of tracers.
+!
+!-----------------------------------------------------------------------
+ subroutine mpas_tracer_advection_coefficients( grid, err, maxLevelCell_in, highOrderAdvectionMask_in, boundaryCell_in )!{{{
+ use mpas_hash
+ use mpas_sort
+
+ implicit none
+ type (mesh_type) :: grid !< Input: Grid information
+ integer, intent(out) :: err
+ integer, dimension(:), pointer, optional :: maxLevelCell_in !< Input - optional: Index to last active cell
+ integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !< Input - optional: Mask for high order advection
+ integer, dimension(:,:), pointer, optional :: boundaryCell_in
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge, highOrderAdvectionMask, 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, nCells
+ logical :: addcell, highOrderAdvection
+
+ type (hashtable) :: cell_hash
+
+ 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
+
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ if(present(maxLevelCell_in)) then
+ maxLevelCell => maxLevelCell_in
+ else
+ allocate(maxLevelCell(nCells))
+ maxLevelCell(:) = nVertLevels
+ end if
+
+ if(present(highOrderAdvectionMask_in)) then
+ highOrderAdvectionMask => highOrderAdvectionMask_in
+ else
+ allocate(highOrderAdvectionMask(nVertLevels, nCells))
+ end if
+
+ if(present(boundaryCell_in)) then
+ boundaryCell => boundaryCell_in
+ else
+ allocate(boundaryCell(nVertLevels, nCells))
+ boundaryCell = 0
+ end if
+
+ allocate(cell_indices(grid % maxEdges2 + 2))
+ allocate(sorted_cell_indices(2, grid % maxEdges2 + 2))
+
+ err = 0
+
+ highOrderAdvectionMask = 0
+
+ ! If 2nd order advection, set masks appropriately.
+ if(config_horiz_tracer_adv_order == 2) then
+ highOrderAdvectionMask = 0
+ return
+ end if
+
+ do iEdge = 1, grid % nEdges
+ nAdvCellsForEdge(iEdge) = 0
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ 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
+
+ !
+ ! do only if this edge flux is needed to update owned cells
+ !
+ if (cell1 <= grid % nCells .and. cell2 <= 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_3rd(:,iEdge) = 0.
+
+ k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell1))
+ if(k <= 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 <= 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 <= 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 <= 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 <= nAdvCellsForEdge(iEdge)) then
+ adv_coefs(k, iEdge) = adv_coefs(k, iEdge) + 0.5
+ end if
+
+ k = mpas_binary_search(sorted_cell_indices, 2, 1, nAdvCellsForEdge(iEdge), grid % indexToCellID % array(cell2))
+ if(k <= nAdvCellsForEdge(iEdge)) then
+ adv_coefs(k, iEdge) = adv_coefs(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_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)
+
+ end subroutine mpas_tracer_advection_coefficients!}}}
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
! function mpas_tracer_advection_vflux4
!
!> \brief MPAS 4th order vertical tracer advection stencil
Modified: branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_mono.F 2013-03-12 20:45:52 UTC (rev 2592)
+++ branches/ocean_projects/advection_operators/src/operators/mpas_tracer_advection_mono.F 2013-03-12 21:20:35 UTC (rev 2593)
@@ -23,6 +23,7 @@
private
save
+ logical :: vert2ndOrder, vert3rdOrder, vert4thOrder
real (kind=RKIND) :: coef_3rd_order
public :: mpas_tracer_advection_mono_tend, &
@@ -43,35 +44,33 @@
!> Both horizontal and vertical.
!
!-----------------------------------------------------------------------
- subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, grid, tend_h, tend, dt_in, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in, verticalDivergenceFactor_in, edgeSignOnCell_in)!{{{
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) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Tendency for thickness field
- real (kind=RKIND), intent(in) :: dt !< Input: Timestep
type (mesh_type), intent(in) :: grid !< Input: Grid information
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
+ real (kind=RKIND), intent(in), optional :: dt_in !< Input - optional: Timestep
+ integer, dimension(:), pointer, optional :: maxLevelCell_in !< Input - optional: Index to max level at cell center
+ integer, dimension(:), pointer, optional :: maxLevelEdgeTop_in !< Input - optional: Index to max level at edge with non-land cells on both sides
+ integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !< Input - optional: Mask for high order advection
+ real (kind=RKIND), dimension(:), pointer, optional :: verticalDivergenceFactor_in !< Input - optional: Denominator for vertical divergence. Given as an inverse which can be multiplied.
+ integer, dimension(:,:), pointer, optional :: edgeSignOnCell_in !< Input - optional: Sign of edges on cell, for bit-reproducibility.
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
- integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
+ integer :: i, j, iCell, iEdge, k, iTracer, cell1, cell2
+ integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve, maxEdges
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask, edgeSignOnCell, edgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, edgeSignOnCell, edgesOnCell
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) :: verticalWeightK, verticalWeightKm1, dt
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, verticalDivergenceFactor
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, 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
@@ -79,6 +78,14 @@
type (field2dReal), pointer :: high_order_horiz_flux_field
+ ! Get dimensions
+ nCells = grid % nCells
+ nCellsSolve = grid % nCellsSolve
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+ maxEdges = grid % maxEdges
+ num_tracers = size(tracers,dim=1)
+
! Initialize pointers
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -89,34 +96,65 @@
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
- adv_coefs_2nd => grid % adv_coefs_2nd % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- highOrderAdvectionMask => grid % highOrderAdvectionMask % array
- lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
- edgeSignOnCell => grid % edgeSignOnCell % array
edgesOnCell => grid % edgesOnCell % array
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers,dim=1)
+ if(present(edgeSignOnCell_in)) then
+ edgeSignOnCell => edgeSignOnCell_in
+ else
+ allocate(edgeSignOnCell(maxEdges, nCells+1))
+ do iCell = 1, nCells
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ if(cellsOnEdge(1, iEdge) == iCell) then
+ edgeSignOnCell(i, iCell) = -1
+ else
+ edgeSignOnCell(i, iCell) = 1
+ end if
+ end do
+ end do
+ endif
- allocate(high_order_horiz_flux_field)
- nullify(high_order_horiz_flux_field % next)
- high_order_horiz_flux_field % block => grid % block
- high_order_horiz_flux_field % sendList => grid % xEdge % sendList
- high_order_horiz_flux_field % recvList => grid % xEdge % recvList
- high_order_horiz_flux_field % copyList => 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 => high_order_horiz_flux_field % array
+ ! Setup arrays from optional arguments
+ if(present(maxLevelCell_in)) then
+ maxLevelCell => maxLevelCell_in
+ else
+ allocate(maxLevelCell(nCells+1))
+ maxLevelCell(:) = nVertLevels
+ end if
+ if(present(maxLevelEdgeTop_in)) then
+ maxLevelEdgeTop => maxLevelEdgeTop_in
+ else
+ allocate(maxLevelEdgeTop(nEdges+1))
+ maxLevelEdgeTop(:) = nVertLevels
+ end if
+
+ if(present(highOrderAdvectionMask_in)) then
+ highOrderAdvectionMask => 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
+
+ if(present(verticalDivergenceFactor_in)) then
+ verticalDivergenceFactor => verticalDivergenceFactor_in
+ else
+ allocate(verticalDivergenceFactor(nVertLevels))
+ verticalDivergenceFactor = 1.0_RKIND
+ end if
+
+ if(present(dt_in)) then
+ dt = dt_in
+ else
+ dt = 1.0_RKIND
+ end if
+
! allocate nCells arrays
-
allocate(tracer_new(nVertLevels, nCells+1))
allocate(tracer_cur(nVertLevels, nCells+1))
allocate(upwind_tendency(nVertLevels, nCells+1))
@@ -127,7 +165,7 @@
allocate(flux_outgoing(nVertLevels, nCells+1))
! allocate nEdges arrays
-! allocate(high_order_horiz_flux(nVertLevels, nEdges))
+ allocate(high_order_horiz_flux(nVertLevels, nEdges+1))
! allocate nVertLevels+1 and nCells arrays
allocate(high_order_vert_flux(nVertLevels+1, nCells))
@@ -162,7 +200,7 @@
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
+ 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))
@@ -170,14 +208,22 @@
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) = mpas_tracer_advection_vflux3( 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 )
+ if(vert4thOrder) then
+ high_order_vert_flux(k, iCell) = mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
+ 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), &
+ 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 = maxLevelCell(iCell)
+ 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))
@@ -195,11 +241,22 @@
! Compute the high order horizontal flux
do iEdge = 1, nEdges
+ ! Low order horizontal flux where required/requested
+ do i = 1, 2
+ iCell = cellsOnEdge(i, iEdge)
+ do k = 1, maxLevelCell(iCell)
+ tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1,1) * dvEdge(iEdge) * 0.5_RKIND
+
+ 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
+
+ ! High order horizontal flux where allowed/requested
do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k = 1, maxLevelCell(iCell)
- tracer_weight = lowOrderAdvectionMask(k, iEdge) * adv_coefs_2nd(i,iEdge) &
- + highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,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)
@@ -213,9 +270,11 @@
! Upwind fluxes are accumulated in upwind_tendency
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
- ! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! flux_upwind = max(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
- flux_upwind = min(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
+ if(config_dzdk_positive) then
+ flux_upwind = max(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
+ else
+ flux_upwind = min(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.0_RKIND,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
@@ -226,12 +285,13 @@
! 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.0_RKIND,high_order_vert_flux(k+1,iCell))-max(0.0_RKIND,high_order_vert_flux(k,iCell)))
-! flux_outgoing(k,iCell) = -(max(0.0_RKIND,high_order_vert_flux(k+1,iCell))-min(0.0_RKIND,high_order_vert_flux(k,iCell)))
-
- flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
+ if(config_dzdk_positive) then
+ flux_incoming (k,iCell) = -(min(0.0_RKIND,high_order_vert_flux(k+1,iCell))-max(0.0_RKIND,high_order_vert_flux(k,iCell)))
+ flux_outgoing(k,iCell) = -(max(0.0_RKIND,high_order_vert_flux(k+1,iCell))-min(0.0_RKIND,high_order_vert_flux(k,iCell)))
+ else
+ flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
+ flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
+ end if
end do ! k Loop
end do ! iCell Loop
@@ -303,11 +363,13 @@
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.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
-! + min(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
- flux = max(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
- + min(0.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell))
+ if(config_dzdk_positive) then
+ flux = max(0.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
+ + min(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
+ else
+ flux = max(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
+ + min(0.0_RKIND,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
@@ -318,10 +380,10 @@
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
do k = 1, maxLevelEdgeTop(iEdge)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
if(config_check_tracer_monotonicity) then
- tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ tracer_new(k, iCell) = tracer_new(k, iCell) + verticalDivergenceFactor(k) * edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
end if
end do
end do
@@ -368,7 +430,28 @@
deallocate(flux_outgoing)
deallocate(high_order_horiz_flux)
deallocate(high_order_vert_flux)
- deallocate(high_order_horiz_flux_field)
+
+ ! Deallocate arrays from optional arguments
+ 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
+
+ if(.not.present(verticalDivergenceFactor_in)) then
+ deallocate(verticalDivergenceFactor)
+ end if
+
+ if(.not.present(edgeSignOnCell_in)) then
+ deallocate(edgeSignOnCell)
+ end if
+
end subroutine mpas_tracer_advection_mono_tend!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
@@ -390,10 +473,19 @@
err = 0
- if ( config_horiz_tracer_adv_order == 3) then
+ vert2ndOrder = .false.
+ vert3rdOrder = .false.
+ vert4thOrder = .false.
+
+ if(config_horiz_tracer_adv_order == 2) then
+ coef_3rd_order = 0.0
+ vert2ndOrder = .true.
+ else 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
+ vert3rdOrder = .true.
+ else if(config_horiz_tracer_adv_order == 4) then
coef_3rd_order = 0.0
+ vert4thOrder = .true.
end if
end subroutine mpas_tracer_advection_mono_init!}}}
</font>
</pre>