<p><b>dwj07@fsu.edu</b> 2012-10-05 09:40:59 -0600 (Fri, 05 Oct 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Fixing bit reproducibility of element version of ocean core.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_mpas_core.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_mpas_core.F 2012-10-05 15:28:04 UTC (rev 2197)
+++ branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_mpas_core.F 2012-10-05 15:40:59 UTC (rev 2198)
@@ -1017,7 +1017,7 @@
iEdge = edgesOnVertex(i, iVertex)
! Vector points from vertex 1 to vertex 2
- if(iVertex == verticesOnEdge(1, iVertex)) then
+ if(iVertex == verticesOnEdge(1, iEdge)) then
edgeSignOnVertex(i, iVertex) = -1
else
edgeSignOnVertex(i, iVertex) = 1
Modified: branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_tendency.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_tendency.F 2012-10-05 15:28:04 UTC (rev 2197)
+++ branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_tendency.F 2012-10-05 15:40:59 UTC (rev 2198)
@@ -404,7 +404,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell, kiteIndexOnCell, verticesOnCell
+ 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
@@ -479,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
@@ -580,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)
Modified: branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_time_integration_split.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_time_integration_split.F 2012-10-05 15:28:04 UTC (rev 2197)
+++ branches/omp_blocks/openmp_test/src/core_ocean_elements/mpas_ocn_time_integration_split.F 2012-10-05 15:40:59 UTC (rev 2198)
@@ -477,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)
@@ -486,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
@@ -503,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
</font>
</pre>