<p><b>dwj07@fsu.edu</b> 2012-03-13 16:01:16 -0600 (Tue, 13 Mar 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding computation of deriv_two.<br>
Adding masks for advection order, which is required due to land cells under ocean cells.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-13 21:46:10 UTC (rev 1628)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-13 22:01:16 UTC (rev 1629)
@@ -158,6 +158,7 @@
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 - -
+var persistent integer advectionMasks ( TWO nVertLevels nEdges ) 0 - advectionMasks mesh - -
% !! NOTE: the following arrays are needed to allow the use
% !! of the module_advection.F w/o alteration
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F 2012-03-13 21:46:10 UTC (rev 1628)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F 2012-03-13 22:01:16 UTC (rev 1629)
@@ -9,7 +9,7 @@
contains
- subroutine ocn_initialize_advection_rk( grid )!{{{
+ subroutine ocn_initialize_advection_rk( grid, err )!{{{
!
! compute the cell coefficients for the polynomial fit.
@@ -19,6 +19,7 @@
implicit none
type (mesh_type), intent(in) :: grid
+ integer, intent(out) :: err
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
integer, dimension(:,:), pointer :: advCells
@@ -50,9 +51,7 @@
integer :: cell1, cell2
integer, parameter :: polynomial_order = 2
-! logical, parameter :: debug = .true.
logical, parameter :: debug = .false.
-! logical, parameter :: least_squares = .false.
logical, parameter :: least_squares = .true.
logical :: add_the_cell, do_the_cell
@@ -61,8 +60,15 @@
real (kind=RKIND) :: rcell, cos2t, costsint, sin2t
real (kind=RKIND), dimension(grid%maxEdges) :: angle_2d
-!---
+!---
+ err = 0
+ if(polynomial_order > 2) then
+ write (*,*) 'Polynomial for second derivitave can only be 2'
+ err = 1
+ return
+ end if
+
pii = 2.*asin(1.0)
advCells => grid % advCells % array
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-03-13 21:46:10 UTC (rev 1628)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-13 22:01:16 UTC (rev 1629)
@@ -138,7 +138,10 @@
block => domain % blocklist
do while (associated(block))
- call mpas_init_block(block, block % mesh, dt)
+ call mpas_init_block(block, block % mesh, dt, err)
+ if(err.eq.1) then
+ call mpas_dmpar_abort(dminfo)
+ endif
block % state % time_levs(1) % state % xtime % scalar = startTimeStamp
block => block % next
end do
@@ -217,20 +220,23 @@
end subroutine ocn_simulation_clock_init!}}}
- subroutine mpas_init_block(block, mesh, dt)!{{{
+ subroutine mpas_init_block(block, mesh, dt, err)!{{{
use mpas_grid_types
use mpas_rbf_interpolation
use mpas_vector_reconstruction
use mpas_ocn_tracer_advection
+ use ocn_advection
implicit none
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
+ integer, intent(out) :: err
integer :: i, iEdge, iCell, k
+ call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh)
call ocn_time_average_init(block % state % time_levs(1) % state)
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2012-03-13 21:46:10 UTC (rev 1628)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2012-03-13 22:01:16 UTC (rev 1629)
@@ -54,11 +54,12 @@
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 :: advectionMasks
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
integer, dimension(:), pointer :: cell_list, ordered_cell_list
- integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
+ integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell, k, nVertLevels
logical :: addcell
deriv_two => grid % deriv_two % array
@@ -67,16 +68,34 @@
cellsOnCell => grid % cellsOnCell % array
cellsOnEdge => grid % cellsOnEdge % array
advCellsForEdge => grid % advCellsForEdge % array
+ boundaryCell => grid % boundaryCell % array
+ avectionMasks => grid % advectionMasks % array
nEdgesOnCell => grid % nEdgesOnCell % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+ nVertLevels = grid % nVertLevels
+
allocate(cell_list(grid % maxEdges2 + 2))
allocate(ordered_cell_list(grid % maxEdges2 + 2))
+ avectionMasks = 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
+ advectionMasks(1, k, iEdge) = 1 ! Need to use low order advection for this edge
+ else
+ advectionMasks(2, k, iEdge) = 1 ! Can use high order advection for this edge
+ end if
+
+ write(*,*) iedge, k, advectionMasks(1,k,iEdge) + advectionMasks(2,k,iEdge)
+ end do
+
!
! do only if this edge flux is needed to update owned cells
!
@@ -202,6 +221,8 @@
deallocate(cell_list)
deallocate(ordered_cell_list)
+ write(*,*) 'Min max of sum, ', minval(advectionMasks(1,:,:)+advectionMasks(2,:,:)), maxval(advectionMasks(1,:,:) + advectionMasks(2,:,:))
+
end subroutine mpas_ocn_tracer_advection_coefficients!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-03-13 21:46:10 UTC (rev 1628)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-03-13 22:01:16 UTC (rev 1629)
@@ -62,6 +62,7 @@
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge
+ integer, dimension(:,:,:), pointer :: advectionMasks
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
@@ -89,6 +90,7 @@
adv_coefs_3rd => grid % adv_coefs_3rd % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ advectionMasks => grid % advectionMasks % array
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
@@ -179,7 +181,11 @@
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))
+ ! AdvectionMasks defines the mask for low or high order advection (1,:,:) is the low order mask, while (2,:,:) is the high order mask
+ tracer_weight = advectionMasks(1, k, iEdge) * 0.5 * dvEdge(iEdge) &
+ + advectionMasks(2, k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+
+ tracer_weight = uh(k,iEdge)*tracer_weight
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
end do ! k loop
end do ! i loop over nAdvCellsForEdge
</font>
</pre>