<p><b>dwj07@fsu.edu</b> 2012-10-10 11:44:59 -0600 (Wed, 10 Oct 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Creating bit-reproducible ocean core.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/option3_b4b_test/src/core_ocean
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean 2012-10-10 17:44:59 UTC (rev 2202)
Property changes on: branches/ocean_projects/option3_b4b_test/src/core_ocean
___________________________________________________________________
Added: svn:mergeinfo
## -0,0 +1,24 ##
+/branches/atmos_physics/src/core_ocean:1672-1846
+/branches/cam_mpas_nh/src/core_ocean:1260-1270
+/branches/ocean_projects/ale_split_exp/src/core_ocean:1437-1483
+/branches/ocean_projects/ale_vert_coord/src/core_ocean:1225-1383
+/branches/ocean_projects/ale_vert_coord_new/src/core_ocean:1387-1428
+/branches/ocean_projects/gmvar/src/core_ocean:1214-1514,1517-1738
+/branches/ocean_projects/imp_vert_mix_error/src/core_ocean:1847-1887
+/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
+/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
+/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
+/branches/ocean_projects/split_explicit_mrp/src/core_ocean:1134-1138
+/branches/ocean_projects/split_explicit_timestepping/src/core_ocean:1044-1097
+/branches/ocean_projects/vert_adv_mrp/src/core_ocean:704-745
+/branches/ocean_projects/vol_cons_RK_imp_mix/src/core_ocean:1965-1992
+/branches/ocean_projects/zstar_restart_new/src/core_ocean:1762-1770
+/branches/omp_blocks/block_decomp/src/core_ocean:1374-1569
+/branches/omp_blocks/ddt_reorg/src/core_ocean:1301-1414
+/branches/omp_blocks/halo/src/core_ocean:1570-1638
+/branches/omp_blocks/io/src/core_ocean:1639-1787
+/branches/omp_blocks/multiple_blocks/src/core_ocean:1803-2084
+/branches/omp_blocks/openmp_test/src/core_ocean:2107-2144
+/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
+/branches/source_renaming/src/core_ocean:1082-1113
+/branches/time_manager/src/core_ocean:924-962
\ No newline at end of property
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry 2012-10-10 17:44:59 UTC (rev 2202)
@@ -307,3 +307,8 @@
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+
+% Sign fields, for openmp and bit reproducibility without branching statements.
+var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
+var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
+var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -249,6 +249,7 @@
integer :: i, iEdge, iCell, k
integer :: err1
+ call ocn_setup_sign_and_index_fields(mesh)
call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh, err1)
err = ior(err, err1)
@@ -960,6 +961,72 @@
end subroutine ocn_compute_mesh_scaling!}}}
+ subroutine ocn_setup_sign_and_index_fields(mesh)!{{{
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell, edgesOnVertex, cellsOnVertex, cellsOnEdge, verticesOnCell, verticesOnEdge
+ integer, dimension(:,:), pointer :: edgeSignOnCell, edgeSignOnVertex, kiteIndexOnCell
+
+ integer :: nCells, nEdges, nVertices, vertexDegree
+ integer :: iCell, iEdge, iVertex, i, j, k
+
+ nCells = mesh % nCells
+ nEdges = mesh % nEdges
+ nVertices = mesh % nVertices
+ vertexDegree = mesh % vertexDegree
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ edgesOnCell => mesh % edgeSOnCell % array
+ edgesOnVertex => mesh % edgesOnVertex % array
+ cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ verticesOnCell => mesh % verticesOnCell % array
+ verticesOnEdge => mesh % verticesOnEdge % array
+ edgeSignOnCell => mesh % edgeSignOnCell % array
+ edgeSignOnVertex => mesh % edgeSignOnVertex % array
+ kiteIndexOnCell => mesh % kiteIndexOnCell % array
+
+ edgeSignOnCell = 0.0_RKIND
+ edgeSignOnVertex = 0.0_RKIND
+ kiteIndexOnCell = 0.0_RKIND
+
+ do iCell = 1, nCells
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+
+ ! Vector points from cell 1 to cell 2
+ if(iCell == cellsOnEdge(1, iEdge)) then
+ edgeSignOnCell(i, iCell) = -1
+ else
+ edgeSignOnCell(i, iCell) = 1
+ end if
+
+ do j = 1, vertexDegree
+ if(cellsOnVertex(j, iVertex) == iCell) then
+ kiteIndexOnCell(i, iCell) = j
+ end if
+ end do
+ end do
+ end do
+
+ do iVertex = 1, nVertices
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+
+ ! Vector points from vertex 1 to vertex 2
+ if(iVertex == verticesOnEdge(1, iEdge)) then
+ edgeSignOnVertex(i, iVertex) = -1
+ else
+ edgeSignOnVertex(i, iVertex) = 1
+ end if
+ end do
+ end do
+
+ end subroutine ocn_setup_sign_and_index_fields!}}}
+
end module mpas_core
! vim: foldmethod=marker
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -404,7 +404,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell
+ verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell, kiteIndexOnCell, verticesOnCell, edgeSignOnVertex, edgeSignOnCell, edgesOnCell
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, coef_3rd_order, r_tmp, invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, invLength, h_vertex
@@ -454,6 +454,7 @@
verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
edgesOnVertex => grid % edgesOnVertex % array
dcEdge => grid % dcEdge % array
@@ -468,6 +469,8 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
+ kiteIndexOnCell => grid % kiteIndexOnCell % array
+ verticesOnCell => grid % verticesOnCell % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -477,7 +480,10 @@
boundaryCell => grid % boundaryCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
@@ -578,39 +584,67 @@
divergence(:,:) = 0.0
ke(:,:) = 0.0
v(:,:) = 0.0
+ ! DWJ 10/05/12 Bit Reproducible Versions
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexBot(iVertex)
+ r_tmp = dcEdge(iEdge) * u(k, iEdge)
+
+ circulation(k, iVertex) = circulation(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp
+ vorticity(k, iVertex) = vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp * invAreaTri1
+ end do
+ end do
+ end do
+
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ r_tmp = dvEdge(iEdge) * u(k, iEdge) * invAreaCell1
+
+ divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * r_tmp
+ ke(k, iCell) = ke(k, iCell) + 0.25 * r_tmp * dcEdge(iEdge) * u(k,iEdge)
+ end do
+ end do
+ end do
+
+ ! DWJ 10/05/12 Non Bit Reproducible Version
do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+! vertex1 = verticesOnEdge(1,iEdge)
+! vertex2 = verticesOnEdge(2,iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
+! invAreaTri1 = 1.0 / areaTriangle(vertex1)
+! invAreaTri2 = 1.0 / areaTriangle(vertex2)
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0)
- invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0)
+! !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
+! invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0)
+! invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0)
- do k=1,maxLevelEdgeBot(iEdge)
- ! Compute circulation and relative vorticity at each vertex
- r_tmp = dcEdge(iEdge) * u(k,iEdge)
- circulation(k,vertex1) = circulation(k,vertex1) - r_tmp
- circulation(k,vertex2) = circulation(k,vertex2) + r_tmp
+! do k=1,maxLevelEdgeBot(iEdge)
+! ! Compute circulation and relative vorticity at each vertex
+! r_tmp = dcEdge(iEdge) * u(k,iEdge)
+! circulation(k,vertex1) = circulation(k,vertex1) - r_tmp
+! circulation(k,vertex2) = circulation(k,vertex2) + r_tmp
- vorticity(k, vertex1) = vorticity(k, vertex1) - r_tmp * invAreaTri1
- vorticity(k, vertex2) = vorticity(k, vertex2) + r_tmp * invAreaTri2
+! vorticity(k, vertex1) = vorticity(k, vertex1) - r_tmp * invAreaTri1
+! vorticity(k, vertex2) = vorticity(k, vertex2) + r_tmp * invAreaTri2
- ! Compute the divergence at each cell center
- r_tmp = dvEdge(iEdge) * u(k, iEdge)
- divergence(k,cell1) = divergence(k,cell1) + r_tmp * invAreaCell1
- divergence(k,cell2) = divergence(k,cell2) - r_tmp * invAreaCell2
+! ! Compute the divergence at each cell center
+! r_tmp = dvEdge(iEdge) * u(k, iEdge)
+! divergence(k,cell1) = divergence(k,cell1) + r_tmp * invAreaCell1
+! divergence(k,cell2) = divergence(k,cell2) - r_tmp * invAreaCell2
- ! Compute kinetic energy in each cell
- r_tmp = r_tmp * dcEdge(iEdge) * u(k,iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * r_tmp * invAreaCell1
- ke(k,cell2) = ke(k,cell2) + 0.25 * r_tmp * invAreaCell2
- end do
+! ! Compute kinetic energy in each cell
+! r_tmp = r_tmp * dcEdge(iEdge) * u(k,iEdge)
+! ke(k,cell1) = ke(k,cell1) + 0.25 * r_tmp * invAreaCell1
+! ke(k,cell2) = ke(k,cell2) + 0.25 * r_tmp * invAreaCell2
+! end do
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
@@ -627,30 +661,53 @@
!
! Compute kinetic energy in each vertex
!
+ !DWJ 09/19/12 Bit Reproducible Version
kev(:,:) = 0.0; kevc(:,:) = 0.0
- do iEdge=1,nEdges*ke_vertex_flag
- do k=1,nVertLevels
- r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * u(k, iEdge)**2
- kev(k,verticesOnEdge(1,iEdge)) = kev(k,verticesOnEdge(1,iEdge)) + r_tmp
- kev(k,verticesOnEdge(2,iEdge)) = kev(k,verticesOnEdge(2,iEdge)) + r_tmp
- end do
+ do iVertex = 1, nVertices*ke_vertex_flag
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * 0.25 / areaTriangle(iVertex)
+ do k = 1, nVertLevels
+ kev(k, iVertex) = kev(k, iVertex) + r_tmp * u(k, iEdge)**2
+ end do
+ end do
end do
- do iVertex = 1,nVertices*ke_vertex_flag
- do k=1,nVertLevels
- kev(k,iVertex) = kev(k,iVertex) / areaTriangle(iVertex) * 0.25
- enddo
- enddo
- do iVertex = 1, nVertices*ke_vertex_flag
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
- do k=1,nVertLevels
- kevc(k,iCell) = kevc(k,iCell) + kiteAreasOnVertex(i, iVertex) * kev(k, iVertex) * invAreaCell1
- enddo
- enddo
- enddo
+ do iCell = 1, nCells*ke_vertex_flag
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, nVertLevels
+ kevc(k, iCell) = kevc(k, iCell) + kiteAreasOnVertex(j, iVertex) * kev(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges*ke_vertex_flag
+! do k=1,nVertLevels
+! r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * u(k, iEdge)**2
+! kev(k,verticesOnEdge(1,iEdge)) = kev(k,verticesOnEdge(1,iEdge)) + r_tmp
+! kev(k,verticesOnEdge(2,iEdge)) = kev(k,verticesOnEdge(2,iEdge)) + r_tmp
+! end do
+! end do
+! do iVertex = 1,nVertices*ke_vertex_flag
+! do k=1,nVertLevels
+! kev(k,iVertex) = kev(k,iVertex) / areaTriangle(iVertex) * 0.25
+! enddo
+! enddo
+! do iVertex = 1, nVertices*ke_vertex_flag
+! do i=1,grid % vertexDegree
+! iCell = cellsOnVertex(i,iVertex)
+! !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
+! invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
+! do k=1,nVertLevels
+! kevc(k,iCell) = kevc(k,iCell) + kiteAreasOnVertex(i, iVertex) * kev(k, iVertex) * invAreaCell1
+! enddo
+! enddo
+! enddo
+
!
! Compute kinetic energy in each cell by blending ke and kevc
!
@@ -692,28 +749,51 @@
Vor_cell(:,:) = 0.0
Vor_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- iEdge = edgesOnVertex(i,iVertex)
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iEdge = 1, nEdges
+ vertex1 = verticesOnEdge(1, iEdge)
+ vertex2 = verticesOnEdge(2, iEdge)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ Vor_edge(k, iEdge) = 0.5 * (Vor_vertex(k, vertex1) + Vor_vertex(k, vertex2))
+ end do
+ end do
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
- ! Compute pv at cell centers
- ! ( this computes Vor_cell for all real cells and distance-1 ghost cells )
- do k = 1,maxLevelCell(iCell)
- Vor_cell(k,iCell) = Vor_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
- enddo
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ Vor_cell(k, iCell) = Vor_cell(k, iCell) + kiteAreasOnVertex(j, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
- ! Compute pv at the edges
- ! ( this computes Vor_edge at all edges bounding real cells )
- do k=1,maxLevelEdgeBot(iEdge)
- Vor_edge(k,iEdge) = Vor_edge(k,iEdge) + 0.5 * Vor_vertex(k,iVertex)
- enddo
- enddo
- enddo
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iVertex = 1,nVertices
+! do i=1,vertexDegree
+! iCell = cellsOnVertex(i,iVertex)
+! iEdge = edgesOnVertex(i,iVertex)
+
+! !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
+! invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
+
+! ! Compute pv at cell centers
+! ! ( this computes Vor_cell for all real cells and distance-1 ghost cells )
+! do k = 1,maxLevelCell(iCell)
+! Vor_cell(k,iCell) = Vor_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
+! enddo
+
+! ! Compute pv at the edges
+! ! ( this computes Vor_edge at all edges bounding real cells )
+! do k=1,maxLevelEdgeBot(iEdge)
+! Vor_edge(k,iEdge) = Vor_edge(k,iEdge) + 0.5 * Vor_vertex(k,iVertex)
+! enddo
+! enddo
+! enddo
+
! gradVor_n(:,:) = 0.0
! gradVor_t(:,:) = 0.0
do iEdge = 1,nEdges
@@ -872,7 +952,7 @@
type (mesh_type), intent(in) :: grid !< Input: Grid information
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum, invAreaCell
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
@@ -885,7 +965,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
+ boundaryEdge, boundaryCell, edgeSignOnCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
@@ -895,8 +975,11 @@
uTransport => s2 % uTransport % array
wTop => s2 % wTop % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
@@ -913,21 +996,36 @@
! Compute div(h^{edge} u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- div_hu(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- div_hu(k,cell1) = div_hu(k,cell1) + flux
- div_hu(k,cell2) = div_hu(k,cell2) - flux
- end do
- end do
+ div_hu(:,:) = 0.0_RKIND
+ !DWJ 09/19/12 Bit Reproducible Version
do iCell=1,nCells
- div_hu_btr(iCell) = 0.0
+ invAreaCell = 1.0_RKIND / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = uTransport(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ div_hu(k, iCell) = div_hu(k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell
+ end do
+ end do
+ end do
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! do k=1,maxLevelEdgeBot(iEdge)
+! flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+! div_hu(k,cell1) = div_hu(k,cell1) + flux
+! div_hu(k,cell2) = div_hu(k,cell2) - flux
+! end do
+! end do
+
+ do iCell=1,nCells
+ div_hu_btr(iCell) = 0.0_RKIND
do k=1,maxLevelCell(iCell)
- div_hu(k,iCell) = div_hu(k,iCell) / areaCell(iCell)
+! div_hu(k,iCell) = div_hu(k,iCell) / areaCell(iCell)
div_hu_btr(iCell) = div_hu_btr(iCell) + div_hu(k,iCell)
end do
end do
@@ -938,13 +1036,13 @@
!dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
- wTop=0.0
+ wTop=0.0_RKIND
else ! zlevel or zstar type vertical grid
do iCell=1,nCells
- hSum = 0.0
+ hSum = 0.0_RKIND
do k=1,maxLevelCell(iCell)
h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
hSum = hSum + zstarWeight(k)*h(k,iCell)
@@ -956,8 +1054,8 @@
! Vertical velocity through layer interface at top and
! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ wTop(1,iCell) = 0.0_RKIND
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND
do k=maxLevelCell(iCell),2,-1
wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
end do
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -101,13 +101,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, i
integer :: iCell, nCells
- integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell
- real (kind=RKIND) :: flux, invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: flux, invAreaCell, invAreaCell1, invAreaCell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
!-----------------------------------------------------------------
@@ -130,21 +130,38 @@
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux
- tend(k,cell2) = tend(k,cell2) + flux
- end do
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ tend(k, iCell) = tend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell
+ end do
+ end do
end do
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
- end do
- end do
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! do k=1,maxLevelEdgeBot(iEdge)
+! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+! tend(k,cell1) = tend(k,cell1) - flux
+! tend(k,cell2) = tend(k,cell2) + flux
+! end do
+! end do
+! do iCell=1,nCells
+! do k=1,maxLevelCell(iCell)
+! tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
+! end do
+! end do
+
!--------------------------------------------------------------------
end subroutine ocn_thick_hadv_tend!}}}
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -408,6 +408,31 @@
! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1, iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2, iEdge)
+
+ sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+
+ flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOncell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -420,8 +445,8 @@
+ config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
+! block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
+! block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ FBtr_coeff*flux
@@ -452,6 +477,8 @@
block => domain % blocklist
do while (associated(block))
+ allocate(utemp(block % mesh % nEdges+1))
+ uTemp(:) = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -461,7 +488,8 @@
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
- * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ !* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ * uTemp(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -478,6 +506,7 @@
+ dt/config_n_btr_subcycles *(CoriolisTerm - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
end do
+ deallocate(uTemp)
block => block % next
end do ! block
@@ -502,6 +531,35 @@
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ ! SSH is a linear combination of SSHold and SSHnew.
+ sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+ sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOnCell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -519,8 +577,8 @@
+ config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
+! block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
+! block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
end do
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -64,7 +64,7 @@
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
+ 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
@@ -95,6 +95,8 @@
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
@@ -230,8 +232,8 @@
! flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
! flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
- flux_incoming (k, iCell) = max(0.0, high_order_vert_flux(k+1, iCell)) - min(0.0, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0, high_order_vert_flux(k+1, iCell)) - max(0.0, high_order_vert_flux(k, iCell))
+ 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
@@ -239,6 +241,7 @@
! 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
+ !DWJ 09/19/12 Bit Reproducible Version
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -249,18 +252,52 @@
do k = 1, maxLevelEdgeTop(iEdge)
flux_upwind = dvEdge(iEdge) * (max(0.,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.,uh(k,iEdge))*tracer_cur(k,cell2))
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
-
- upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
- upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
-
- ! Accumulate remaining high order fluxes
- flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
- flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
end do ! k loop
end do ! iEdge loop
+ !DWJ 09/19/12 Bit Reproducible Version
+ 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) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ flux_incoming(k,iCell) = flux_incoming(k,iCell) - edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ end do
+ end do
+ end do
+
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge = 1, nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+
+! invAreaCell1 = 1.0 / areaCell(cell1)
+! invAreaCell2 = 1.0 / areaCell(cell2)
+
+! do k = 1, maxLevelEdgeTop(iEdge)
+! flux_upwind = dvEdge(iEdge) * (max(0.,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.,uh(k,iEdge))*tracer_cur(k,cell2))
+! high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
+
+! upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
+! upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
+
+! ! Accumulate remaining high order fluxes
+! flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
+! flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
+! flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
+! flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
+! end do ! k loop
+! end do ! iEdge loop
+
! Build the factors for the FCT
! Computed using the bounds that were computed previously, and the bounds on the newly updated value
! Factors are placed in the flux_incoming and flux_outgoing arrays
@@ -304,24 +341,40 @@
end do ! iCell loop
! Accumulate the scaled high order horizontal tendencies
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ !DWJ 09/19/12 Bit Reproducible Version
+ 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
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+ if(config_check_monotonicity) then
+ tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ end if
+ end do
+ end do
+ end do
- if (config_check_monotonicity) then
- !tracer_new holds a tendency for now.
- tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
- end if
- end do ! k loop
- end do ! iEdge loop
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge = 1, nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! invAreaCell1 = 1.0 / areaCell(cell1)
+! invAreaCell2 = 1.0 / areaCell(cell2)
+! do k=1,maxLevelEdgeTop(iEdge)
+! tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+! tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+
+! if (config_check_monotonicity) then
+! !tracer_new holds a tendency for now.
+! tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+! tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+! end if
+! end do ! k loop
+! end do ! iEdge loop
+
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -107,13 +107,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, nVertLevels, cell1, cell2
- integer :: k, iTracer, num_tracers
+ integer :: iCell, iEdge, nCells, nEdges, nVertLevels, cell1, cell2
+ integer :: i, k, iTracer, num_tracers
integer, dimension(:,:), allocatable :: boundaryMask
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2
real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp
@@ -134,6 +134,7 @@
if (.not.del2On) return
nEdges = grid % nEdges
+ nCells = grid % nCells
nVertLevels = grid % nVertLevels
num_tracers = size(tracers, dim=1)
@@ -145,31 +146,61 @@
dcEdge => grid % dcEdge % array
meshScalingDel2 => grid % meshScalingDel2 % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ !DWJ 09/19/12 Bit Reproducible Version
+ 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)
- r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers
+ r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers
! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)
+ tracer_turb_flux = tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1)
! div(h \kappa_2 </font>
<font color="red">abla \phi) at cell center
- flux = h_edge(k,iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
+ flux = h_edge(k, iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
- end do
- end do
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell1
+ end do
+ end do
+ end do
end do
+
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! invAreaCell1 = 1.0/areaCell(cell1)
+! invAreaCell2 = 1.0/areaCell(cell2)
+
+! r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+
+! do k=1,maxLevelEdgeTop(iEdge)
+! do iTracer=1,num_tracers
+! ! \kappa_2 </font>
<font color="blue">abla \phi on edge
+! tracer_turb_flux = tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)
+
+! ! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
+! flux = h_edge(k,iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
+
+! tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
+! tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
+! end do
+! end do
+
+! end do
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -108,10 +108,10 @@
!-----------------------------------------------------------------
integer :: iEdge, nEdges, num_tracers, nVertLevels, nCells
- integer :: iTracer, k, iCell, cell1, cell2
+ integer :: iTracer, k, iCell, cell1, cell2, i
- integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
- integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2
@@ -148,58 +148,107 @@
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
delsq_tracer(:,:,:) = 0.0
! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ invdcEdge = dvEdge(iEdge) / dcEdge(iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ r_tmp1 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell1)
+ r_tmp2 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k, iEdge)
+ delsq_tracer(iTracer, k, iCell) = delsq_tracer(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * (r_tmp2 - r_tmp1) * invAreaCell1
+ end do
+ end do
+ end do
+ end do
- r_tmp1 = dvEdge(iEdge) * h_edge(k,iEdge) * invdcEdge
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
- r_tmp2 = r_tmp1 * tracers(iTracer,k,cell2)
- r_tmp1 = r_tmp1 * tracers(iTracer,k,cell1)
+! invdcEdge = 1.0 / dcEdge(iEdge)
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) + (r_tmp2 - r_tmp1) * invAreaCell1
- delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) - (r_tmp2 - r_tmp1) * invAreaCell2
- end do
- end do
- end do
+! invAreaCell1 = 1.0 / areaCell(cell1)
+! invAreaCell2 = 1.0 / areaCell(cell2)
- ! second del2: div(h </font>
<font color="red">abla [delsq_tracer]) at cell center
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
+! do iTracer=1,num_tracers * edgeMask(k, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+! r_tmp1 = dvEdge(iEdge) * h_edge(k,iEdge) * invdcEdge
- invdcEdge = 1.0 / dcEdge(iEdge)
+! r_tmp2 = r_tmp1 * tracers(iTracer,k,cell2)
+! r_tmp1 = r_tmp1 * tracers(iTracer,k,cell1)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k,iEdge)
- tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
- * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) &
- * invdcEdge
+! delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) + (r_tmp2 - r_tmp1) * invAreaCell1
+! delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) - (r_tmp2 - r_tmp1) * invAreaCell2
+! end do
+! end do
+! end do
- flux = dvEdge (iEdge) * tracer_turb_flux
+ ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
+ !DWJ 09/19/12 Bit Reproducible Version
+ 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)
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
- enddo
- enddo
+ invdcEdge = meshScalingDel4(iEdge) * dvEdge(iEdge) * eddyDiff4 / dcEdge(iEdge)
+
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
+ tracer_turb_flux = (delsq_tracer(iTracer, k, cell2) - delsq_tracer(iTracer, k, cell1))
+
+ flux = tracer_turb_flux * invdcEdge
+
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * invAreaCell1
+ end do
+ end do
+ end do
end do
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,grid % nEdges
+! cell1 = grid % cellsOnEdge % array(1,iEdge)
+! cell2 = grid % cellsOnEdge % array(2,iEdge)
+
+! invAreaCell1 = 1.0 / areaCell(cell1)
+! invAreaCell2 = 1.0 / areaCell(cell2)
+
+! invdcEdge = 1.0 / dcEdge(iEdge)
+
+! do k=1,maxLevelEdgeTop(iEdge)
+! do iTracer=1,num_tracers * edgeMask(k,iEdge)
+! tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
+! * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) &
+! * invdcEdge
+
+! flux = dvEdge (iEdge) * tracer_turb_flux
+
+! tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
+! tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
+! enddo
+! enddo
+! end do
+
deallocate(delsq_tracer)
!--------------------------------------------------------------------
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-10-10 16:49:23 UTC (rev 2201)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-10-10 17:44:59 UTC (rev 2202)
@@ -427,13 +427,13 @@
!
!-----------------------------------------------------------------
- integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k
+ integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k, i
integer :: cell1, cell2
- integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOncell, edgeSignOnCell
- real (kind=RKIND) :: coef
+ real (kind=RKIND) :: coef, invAreaCell
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: drhoTopOfCell, du2TopOfCell, &
drhoTopOfEdge, du2TopOfEdge
@@ -453,6 +453,9 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
allocate( &
drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
@@ -498,22 +501,35 @@
! interpolate du2TopOfEdge to du2TopOfCell
du2TopOfCell = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- end do
+ !DWJ 09/19/12 Bit Reproducible Version
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+
+ do k = 2, maxLevelEdgeBot(iEdge)
+ du2TopOfCell(k, iCell) = du2TopOfCell(k, iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k, iEdge) * invAreaCell
+ end do
+ end do
end do
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
- end do
- end do
+ !DWJ 09/19/12 Non Bit Reproducible Version
+! do iEdge=1,nEdges
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
+! do k=2,maxLevelEdgeBot(iEdge)
+! du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
+! + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+! du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
+! + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+! end do
+! end do
+! do iCell = 1,nCells
+! do k = 2,maxLevelCell(iCell)
+! du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
+! end do
+! end do
+
! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
! coef = -g/rho_0/2
RiTopOfEdge = 0.0
</font>
</pre>