<p><b>dwj07@fsu.edu</b> 2012-03-13 20:56:14 -0600 (Tue, 13 Mar 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Fixing coefficients for 3rd order advection.<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 22:56:30 UTC (rev 1634)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-14 02:56:14 UTC (rev 1635)
@@ -155,6 +155,7 @@
% Added for monotonic advection scheme
var persistent real adv_coefs ( FIFTEEN nEdges ) 0 - adv_coefs mesh - -
+var persistent real adv_coefs_2nd ( FIFTEEN nEdges ) 0 - adv_coefs_2nd mesh - -
var persistent real adv_coefs_3rd ( FIFTEEN nEdges ) 0 - adv_coefs_3rd mesh - -
var persistent integer advCellsForEdge ( FIFTEEN nEdges ) 0 - advCellsForEdge mesh - -
var persistent integer nAdvCellsForEdge ( nEdges ) 0 - nAdvCellsForEdge mesh - -
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 22:56:30 UTC (rev 1634)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-14 02:56:14 UTC (rev 1635)
@@ -235,9 +235,12 @@
real (kind=RKIND), intent(in) :: dt
integer, intent(out) :: err
integer :: i, iEdge, iCell, k
+ integer :: err1
call ocn_initialize_advection_rk(mesh, err)
- call mpas_ocn_tracer_advection_coefficients(mesh)
+ call mpas_ocn_tracer_advection_coefficients(mesh, err1)
+ err = ior(err, err1)
+
call ocn_time_average_init(block % state % time_levs(1) % state)
call mpas_timer_start("diagnostic solve", .false., initDiagSolveTimer)
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 22:56:30 UTC (rev 1634)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2012-03-14 02:56:14 UTC (rev 1635)
@@ -47,13 +47,14 @@
!> advection of tracers.
!
!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_coefficients( grid )!{{{
+ 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_3rd
+ 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
@@ -63,6 +64,7 @@
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
@@ -79,6 +81,8 @@
allocate(cell_list(grid % maxEdges2 + 2))
allocate(ordered_cell_list(grid % maxEdges2 + 2))
+ err = 0
+
highOrderAdvectionMask = 0
lowOrderAdvectionMask = 0
@@ -138,7 +142,7 @@
ordered_cell_list(i) = cell_list(j)
end if
end do
-! ordered_cell_list(i) = cell_list(j_in)
+! ordered_cell_list(i) = cell_list(j_in)
cell_list(j_in) = grid % nCells + 3
end do
@@ -150,6 +154,7 @@
! we have the ordered list, now construct coefficients
adv_coefs(:,iEdge) = 0.
+ adv_coefs_2nd(:,iEdge) = 0.
adv_coefs_3rd(:,iEdge) = 0.
! pull together third and fourth order contributions to the flux
@@ -202,17 +207,20 @@
if( ordered_cell_list(j) == cell1 ) j_in = j
enddo
adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+ adv_coefs_2nd(j_in,iEdge) = adv_coefs_2nd(j_in,iEdge) + 0.5
j_in = 0
do j=1, n
if( ordered_cell_list(j) == cell2 ) j_in = j
enddo
adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+ adv_coefs_2nd(j_in,iEdge) = adv_coefs_2nd(j_in,iEdge) + 0.5
! multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
do j=1,n
adv_coefs (j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs (j,iEdge)
+ adv_coefs_2nd(j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_2nd(j,iEdge)
adv_coefs_3rd(j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_3rd(j,iEdge)
end do
@@ -223,7 +231,10 @@
deallocate(cell_list)
deallocate(ordered_cell_list)
- write(*,*) 'Min max of sum, ', minval(highOrderAdvectionMask+lowOrderAdvectionMask), maxval(highOrderAdvectionMask+lowOrderAdvectionMask)
+ if (maxval(highOrderAdvectionMask+lowOrderAdvectionMask) > 1) then
+ write(*,*) "Masks don't sum to 1."
+ err = 1
+ endif
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 22:56:30 UTC (rev 1634)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-03-14 02:56:14 UTC (rev 1635)
@@ -68,7 +68,7 @@
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_3rd
+ 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
@@ -86,6 +86,7 @@
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
@@ -181,7 +182,7 @@
do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k = 1, maxLevelCell(iCell)
- tracer_weight = lowOrderAdvectionMask(k, iEdge) * 0.5 * dvEdge(iEdge) &
+ tracer_weight = lowOrderAdvectionMask(k, iEdge) * adv_coefs_2nd(i,iEdge) &
+ highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
tracer_weight = uh(k,iEdge)*tracer_weight
</font>
</pre>