<p><b>dwj07@fsu.edu</b> 2012-08-14 15:33:46 -0600 (Tue, 14 Aug 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding changes to ocean core that supposed bit reproducibility.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/option1_b4b_test/Makefile
===================================================================
--- branches/ocean_projects/option1_b4b_test/Makefile 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/Makefile 2012-08-14 21:33:46 UTC (rev 2101)
@@ -131,7 +131,7 @@
"DEBUG = $(DEBUG)" \
"SERIAL = $(SERIAL)" \
"USE_PAPI = $(USE_PAPI)" \
- "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -D_BITREPRODUCIBLE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
g95:
( $(MAKE) all \
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/Registry 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/Registry 2012-08-14 21:33:46 UTC (rev 2101)
@@ -209,12 +209,12 @@
var persistent real salinityRestore ( nCells ) 0 iro salinityRestore mesh - -
% mrp trying to figure out why these do not appear
-var persistent real windStressMonthly ( nMonths nEdges ) 0 iro windStressMonthly mesh - -
-var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 iro temperatureRestoreMonthly mesh - -
-var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 iro salinityRestoreMonthly mesh - -
+var persistent real windStressMonthly ( nMonths nEdges ) 0 - windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 - temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 - salinityRestoreMonthly mesh - -
% Prognostic variables: read from input, saved in restart, and written to output
-var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
+var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
var persistent real rho ( nVertLevels nCells Time ) 2 iro rho state - -
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
@@ -307,3 +307,7 @@
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 - -
+
+var persistent integer cellPermute ( nCells ) 0 - cellPermute mesh - -
+var persistent integer edgePermute ( nEdges ) 0 - edgePermute mesh - -
+var persistent integer vertexPermute ( nVertices ) 0 - vertexPermute mesh - -
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -396,7 +396,8 @@
type (state_type), intent(inout) :: s !< Input/Output: State information
type (mesh_type), intent(in) :: grid !< Input: Grid information
-
+ integer :: edgeIndex, vertexIndex
+ integer, dimension(:), pointer :: edgePermute, vertexPermute
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
integer :: boundaryMask, velMask, nCells, nEdges, nVertices, nVertLevels, vertexDegree, err
@@ -468,6 +469,8 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
+ edgePermute => grid % edgePermute % array
+ vertexPermute => grid % vertexPermute % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -489,7 +492,7 @@
h_edge = -1.0e34
coef_3rd_order = config_coef_3rd_order
- do iEdge=1,nEdges*hadv2nd
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
@@ -497,72 +500,72 @@
end do
end do
- do iEdge=1,nEdges*hadv3rd
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+! do iEdge=1,nEdges*hadv3rd
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+! d2fdx2_cell1 = 0.0
+! d2fdx2_cell2 = 0.0
- boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
+! boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
+! d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
+! d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
- !-- all edges of cell 1
- do i=1, nEdgesOnCell(cell1) * boundaryMask
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
+! !-- all edges of cell 1
+! do i=1, nEdgesOnCell(cell1) * boundaryMask
+! d2fdx2_cell1 = d2fdx2_cell1 + &
+! deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+! end do
- !-- all edges of cell 2
- do i=1, nEdgesOnCell(cell2) * boundaryMask
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+! !-- all edges of cell 2
+! do i=1, nEdgesOnCell(cell2) * boundaryMask
+! d2fdx2_cell2 = d2fdx2_cell2 + &
+! deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+! end do
- velMask = 2*(abs(transfer(u(k,iEdge) <= 0, velMask))) - 1
+! velMask = 2*(abs(transfer(u(k,iEdge) <= 0, velMask))) - 1
- h_edge(k,iEdge) = 0.5*(h(k,cell1) + h(k,cell2)) - (dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- + velMask * (dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+! h_edge(k,iEdge) = 0.5*(h(k,cell1) + h(k,cell2)) - (dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+! + velMask * (dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end do ! do k
- end do ! do iEdge
+! end do ! do k
+! end do ! do iEdge
- do iEdge=1,nEdges*hadv4th
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+! do iEdge=1,nEdges*hadv4th
+! cell1 = cellsOnEdge(1,iEdge)
+! cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+! do k=1,maxLevelEdgeTop(iEdge)
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+! d2fdx2_cell1 = 0.0
+! d2fdx2_cell2 = 0.0
- boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
+! boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
+! d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
+! d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
- !-- all edges of cell 1
- do i=1, nEdgesOnCell(cell1) * boundaryMask
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
+! !-- all edges of cell 1
+! do i=1, nEdgesOnCell(cell1) * boundaryMask
+! d2fdx2_cell1 = d2fdx2_cell1 + &
+! deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+! end do
- !-- all edges of cell 2
- do i=1, nEdgesOnCell(cell2) * boundaryMask
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+! !-- all edges of cell 2
+! do i=1, nEdgesOnCell(cell2) * boundaryMask
+! d2fdx2_cell2 = d2fdx2_cell2 + &
+! deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+! end do
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+! h_edge(k,iEdge) = &
+! 0.5*(h(k,cell1) + h(k,cell2)) &
+! -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
- end do ! do k
- end do ! do iEdge
+! end do ! do k
+! end do ! do iEdge
!
! set the velocity and height at dummy address
@@ -579,11 +582,16 @@
ke(:,:) = 0.0
v(:,:) = 0.0
do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ vertex1 = verticesOnEdge(1,edgeIndex)
+ vertex2 = verticesOnEdge(2,edgeIndex)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
invAreaTri1 = 1.0 / areaTriangle(vertex1)
invAreaTri2 = 1.0 / areaTriangle(vertex2)
@@ -592,9 +600,9 @@
invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0)
invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0)
- do k=1,maxLevelEdgeBot(iEdge)
+ do k=1,maxLevelEdgeBot(edgeIndex)
! Compute circulation and relative vorticity at each vertex
- r_tmp = dcEdge(iEdge) * u(k,iEdge)
+ r_tmp = dcEdge(edgeIndex) * u(k,edgeIndex)
circulation(k,vertex1) = circulation(k,vertex1) - r_tmp
circulation(k,vertex2) = circulation(k,vertex2) + r_tmp
@@ -602,23 +610,23 @@
vorticity(k, vertex2) = vorticity(k, vertex2) + r_tmp * invAreaTri2
! Compute the divergence at each cell center
- r_tmp = dvEdge(iEdge) * u(k, iEdge)
+ r_tmp = dvEdge(edgeIndex) * u(k, edgeIndex)
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)
+ r_tmp = r_tmp * dcEdge(edgeIndex) * u(k,edgeIndex)
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)
- eoe = edgesOnEdge(i,iEdge)
+ do i=1,nEdgesOnEdge(edgeIndex)
+ eoe = edgesOnEdge(i,edgeIndex)
! mrp 101115 note: in order to include flux boundary conditions,
! the following loop may need to change to maxLevelEdgeBot
- do k = 1,maxLevelEdgeTop(iEdge)
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ do k = 1,maxLevelEdgeTop(edgeIndex)
+ v(k,edgeIndex) = v(k,edgeIndex) + weightsOnEdge(i,edgeIndex) * u(k, eoe)
end do
end do
@@ -629,10 +637,15 @@
!
kev(:,:) = 0.0; kevc(:,:) = 0.0
do iEdge=1,nEdges*ke_vertex_flag
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
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
+ r_tmp = dcEdge(edgeIndex) * dvEdge(edgeIndex) * u(k, edgeIndex)**2
+ kev(k,verticesOnEdge(1,edgeIndex)) = kev(k,verticesOnEdge(1,edgeIndex)) + r_tmp
+ kev(k,verticesOnEdge(2,edgeIndex)) = kev(k,verticesOnEdge(2,edgeIndex)) + r_tmp
end do
end do
do iVertex = 1,nVertices*ke_vertex_flag
@@ -641,12 +654,17 @@
enddo
enddo
do iVertex = 1, nVertices*ke_vertex_flag
+#ifdef _BITREPRODUCIBLE
+ vertexIndex = vertexPermute(iVertex)
+#else
+ vertexIndex = iVertex
+#endif
do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
+ iCell = cellsOnVertex(i,vertexIndex)
!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
+ kevc(k,iCell) = kevc(k,iCell) + kiteAreasOnVertex(i, vertexIndex) * kev(k, vertexIndex) * invAreaCell1
enddo
enddo
enddo
@@ -693,9 +711,14 @@
Vor_cell(:,:) = 0.0
Vor_edge(:,:) = 0.0
do iVertex = 1,nVertices
+#ifdef _BITREPRODUCIBLE
+ vertexIndex = vertexPermute(iVertex)
+#else
+ vertexIndex = iVertex
+#endif
do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- iEdge = edgesOnVertex(i,iVertex)
+ iCell = cellsOnVertex(i,vertexIndex)
+ iEdge = edgesOnVertex(i,vertexIndex)
!dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
@@ -703,13 +726,13 @@
! 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
+ Vor_cell(k,iCell) = Vor_cell(k,iCell) + kiteAreasOnVertex(i, vertexIndex) * Vor_vertex(k, vertexIndex) * 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)
+ Vor_edge(k,iEdge) = Vor_edge(k,iEdge) + 0.5 * Vor_vertex(k,vertexIndex)
enddo
enddo
enddo
@@ -871,6 +894,8 @@
type (state_type), intent(inout) :: s2 !< Input/Output: State 2 information
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ integer :: edgeIndex
+ integer, dimension(:), pointer :: edgePermute
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
@@ -901,6 +926,7 @@
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
zstarWeight => grid % zstarWeight % array
+ edgePermute => grid % edgePermute % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -915,10 +941,15 @@
!
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)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
+ do k=1,maxLevelEdgeBot(edgeIndex)
+ flux = uTransport(k,edgeIndex) * dvEdge(edgeIndex) * h_edge(k,edgeIndex)
div_hu(k,cell1) = div_hu(k,cell1) + flux
div_hu(k,cell2) = div_hu(k,cell2) - flux
end do
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_thick_hadv.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -103,8 +103,10 @@
integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
integer :: iCell, nCells
+ integer :: edgeIndex
integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:,:), pointer :: cellsOnEdge
real (kind=RKIND) :: flux, invAreaCell1, invAreaCell2
@@ -129,12 +131,18 @@
cellsOnEdge => grid % cellsOnEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
+ edgePermute => grid % edgePermute % 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)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
+ do k=1,maxLevelEdgeBot(edgeIndex)
+ flux = u(k,edgeIndex) * dvEdge(edgeIndex) * h_edge(k,edgeIndex)
tend(k,cell1) = tend(k,cell1) - flux
tend(k,cell2) = tend(k,cell2) + flux
end do
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -99,6 +99,8 @@
real (kind=RKIND), dimension(:), allocatable:: uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+ integer :: edgeIndex
+
call mpas_timer_start("se timestep", .false., timer_main)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -409,21 +411,26 @@
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = block % mesh % edgePermute % array(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = block % mesh % cellsOnEdge % array(1,edgeIndex)
+ cell2 = block % mesh % cellsOnEdge % array(2,edgeIndex)
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)
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(edgeIndex)+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)) &
+ flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(edgeIndex) &
+ + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(edgeIndex)) &
* 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(edgeIndex)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(edgeIndex)
- block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ block % state % time_levs(1) % state % FBtr % array(edgeIndex) = block % state % time_levs(1) % state % FBtr % array(edgeIndex) &
+ FBtr_coeff*flux
end do
@@ -452,6 +459,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)
@@ -460,8 +469,9 @@
CoriolisTerm = 0.0
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) &
+ CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
+! * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ * uTemp(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -478,6 +488,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
@@ -503,8 +514,13 @@
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = block % mesh % edgePermute % array(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = block % mesh % cellsOnEdge % array(1,edgeIndex)
+ cell2 = block % mesh % cellsOnEdge % array(2,edgeIndex)
! SSH is a linear combination of SSHold and SSHnew.
sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
@@ -513,16 +529,16 @@
+ 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)
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(edgeIndex)+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)) &
+ flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(edgeIndex) &
+ + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(edgeIndex)) &
* 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(edgeIndex)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(edgeIndex)
- block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
+ block % state % time_levs(1) % state % FBtr % array(edgeIndex) = block % state % time_levs(1) % state % FBtr % array(edgeIndex) + flux
end do
! SSHnew = SSHold + dt/J*(-div(Flux))
@@ -684,6 +700,7 @@
! update halo for thickness and tracer tendencies
call mpas_timer_start("se halo h", .false., timer_halo_h)
call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
call mpas_timer_stop("se halo h", timer_halo_h)
block => domain % blocklist
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -61,8 +61,10 @@
type (mesh_type), intent(in) :: grid !< Input: Grid information
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
+ integer :: edgeIndex
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask
@@ -95,6 +97,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
highOrderAdvectionMask => grid % highOrderAdvectionMask % array
lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
+ edgePermute => grid % edgePermute % array
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
@@ -240,24 +243,29 @@
! 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)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
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
+ do k = 1, maxLevelEdgeTop(edgeIndex)
+ flux_upwind = dvEdge(edgeIndex) * (max(0.,uh(k,edgeIndex))*tracer_cur(k,cell1) + min(0.,uh(k,edgeIndex))*tracer_cur(k,cell2))
+ high_order_horiz_flux(k,edgeIndex) = high_order_horiz_flux(k,edgeIndex) - 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
+ flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,edgeIndex)) * invAreaCell1
+ flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,edgeIndex)) * invAreaCell1
+ flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,edgeIndex)) * invAreaCell2
+ flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,edgeIndex)) * invAreaCell2
end do ! k loop
end do ! iEdge loop
@@ -305,19 +313,24 @@
! Accumulate the scaled high order horizontal tendencies
do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
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
+ do k=1,maxLevelEdgeTop(edgeIndex)
+ tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, edgeIndex) * invAreaCell1
+ tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, edgeIndex) * 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
+ tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, edgeIndex) * invAreaCell1
+ tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, edgeIndex) * invAreaCell2
end if
end do ! k loop
end do ! iEdge loop
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -63,6 +63,9 @@
real (kind=RKIND), dimension(:,:), allocatable :: flux_arr
integer :: nVertLevels, num_tracers
+ integer :: edgeIndex
+ integer, dimension(:), pointer :: edgePermute
+
if (.not. hadvOn) return
cellsOnEdge => grid % cellsOnEdge % array
@@ -75,6 +78,7 @@
adv_coefs_3rd => grid % adv_coefs_3rd % array
highOrderAdvectionMask => grid % highOrderAdvectionMask % array
lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
+ edgePermute => grid % edgePermute % array
nVertLevels = grid % nVertLevels
num_tracers = size(tracers, dim=1)
@@ -83,15 +87,20 @@
! horizontal flux divergence, accumulate in tracer_tend
do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
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 i=1,nAdvCellsForEdge(edgeIndex)
+ iCell = advCellsForEdge(i,edgeIndex)
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.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ tracer_weight = lowOrderAdvectionMask(k, edgeIndex) * adv_coefs_2nd(i,edgeIndex) &
+ + highOrderAdvectionMask(k, edgeIndex) * (adv_coefs(i,edgeIndex) + coef_3rd_order*sign(1.,uh(k,edgeIndex))*adv_coefs_3rd(i,edgeIndex))
do iTracer=1,num_tracers
flux_arr(iTracer,k) = flux_arr(iTracer,k) + tracer_weight* tracers(iTracer,k,iCell)
end do
@@ -100,8 +109,8 @@
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)
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - uh(k,edgeIndex)*flux_arr(iTracer,k)/areaCell(cell1)
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + uh(k,edgeIndex)*flux_arr(iTracer,k)/areaCell(cell2)
end do
end do
end if
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -107,12 +107,14 @@
!
!-----------------------------------------------------------------
+ integer :: edgeIndex
integer :: iEdge, nEdges, nVertLevels, cell1, cell2
integer :: k, iTracer, num_tracers
integer, dimension(:,:), allocatable :: boundaryMask
integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask
real (kind=RKIND) :: invAreaCell1, invAreaCell2
@@ -144,25 +146,31 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
meshScalingDel2 => grid % meshScalingDel2 % array
+ edgePermute => grid % edgePermute % 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)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
invAreaCell1 = 1.0/areaCell(cell1)
invAreaCell2 = 1.0/areaCell(cell2)
- r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+ r_tmp = meshScalingDel2(edgeIndex) * eddyDiff2 * dvEdge(edgeIndex) / dcEdge(edgeIndex)
- do k=1,maxLevelEdgeTop(iEdge)
+ do k=1,maxLevelEdgeTop(edgeIndex)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="black">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
+ flux = h_edge(k,edgeIndex) * tracer_turb_flux * edgeMask(k, edgeIndex) * r_tmp
tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -107,10 +107,12 @@
!
!-----------------------------------------------------------------
+ integer :: edgeIndex
integer :: iEdge, nEdges, num_tracers, nVertLevels, nCells
integer :: iTracer, k, iCell, cell1, cell2
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2
@@ -146,6 +148,8 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
+ edgePermute => grid % edgePermute % array
+
edgeMask => grid % edgeMask % array
allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
@@ -154,18 +158,23 @@
! first del2: div(h </font>
<font color="gray">abla \phi) at cell center
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ invdcEdge = 1.0 / dcEdge(edgeIndex)
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k, iEdge)
+ do k=1,maxLevelEdgeTop(edgeIndex)
+ do iTracer=1,num_tracers * edgeMask(k, edgeIndex)
- r_tmp1 = dvEdge(iEdge) * h_edge(k,iEdge) * invdcEdge
+ r_tmp1 = dvEdge(edgeIndex) * h_edge(k,edgeIndex) * invdcEdge
r_tmp2 = r_tmp1 * tracers(iTracer,k,cell2)
r_tmp1 = r_tmp1 * tracers(iTracer,k,cell1)
@@ -178,21 +187,26 @@
! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = grid % cellsOnEdge % array(1,edgeIndex)
+ cell2 = grid % cellsOnEdge % array(2,edgeIndex)
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ invdcEdge = 1.0 / dcEdge(edgeIndex)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k,iEdge)
- tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
+ do k=1,maxLevelEdgeTop(edgeIndex)
+ do iTracer=1,num_tracers * edgeMask(k,edgeIndex)
+ tracer_turb_flux = meshScalingDel4(edgeIndex) * eddyDiff4 &
* (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) &
* invdcEdge
- flux = dvEdge (iEdge) * tracer_turb_flux
+ flux = dvEdge (edgeIndex) * tracer_turb_flux
tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -112,10 +112,12 @@
!
!-----------------------------------------------------------------
+ integer :: edgeIndex
integer :: iEdge, cell1, cell2, vertex1, vertex2, k
integer :: iCell, iVertex
integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexBot, &
maxLevelCell
integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
@@ -149,6 +151,7 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
edgeMask => grid % edgeMask % array
+ edgePermute => grid % edgePermute % array
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_vorticity(nVertLevels, nVertices+1))
@@ -157,11 +160,16 @@
delsq_divergence(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,edgeIndex)
+ vertex2 = verticesOnEdge(2,edgeIndex)
invAreaTri1 = 1.0 / areaTriangle(vertex1)
invAreaTri2 = 1.0 / areaTriangle(vertex2)
@@ -169,21 +177,21 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- invDcEdge = 1.0 / dcEdge(iEdge)
- invDvEdge = 1.0 / dvEdge(iEdge)
+ invDcEdge = 1.0 / dcEdge(edgeIndex)
+ invDvEdge = 1.0 / dvEdge(edgeIndex)
- do k=1,maxLevelEdgeTop(iEdge)
+ do k=1,maxLevelEdgeTop(edgeIndex)
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="black">abla vorticity
delsq_u = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
-viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
! vorticity using </font>
<font color="red">abla^2 u
- r_tmp = dcEdge(iEdge) * delsq_u
+ r_tmp = dcEdge(edgeIndex) * delsq_u
delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
! Divergence using </font>
<font color="gray">abla^2 u
- r_tmp = dvEdge(iEdge) * delsq_u
+ r_tmp = dvEdge(edgeIndex) * delsq_u
delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
end do
Modified: branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -427,9 +427,11 @@
!
!-----------------------------------------------------------------
+ integer :: edgeIndex
integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k
integer :: cell1, cell2
+ integer, dimension(:), pointer :: edgePermute
integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -453,6 +455,7 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
+ edgePermute => grid % edgePermute % array
allocate( &
drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
@@ -499,13 +502,18 @@
! interpolate du2TopOfEdge to du2TopOfCell
du2TopOfCell = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
+#ifdef _BITREPRODUCIBLE
+ edgeIndex = edgePermute(iEdge)
+#else
+ edgeIndex = iEdge
+#endif
+ cell1 = cellsOnEdge(1,edgeIndex)
+ cell2 = cellsOnEdge(2,edgeIndex)
+ do k=2,maxLevelEdgeBot(edgeIndex)
du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+ + 0.5 * dcEdge(edgeIndex) * dvEdge(edgeIndex) * du2TopOfEdge(k,edgeIndex)
du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
+ + 0.5 * dcEdge(edgeIndex) * dvEdge(edgeIndex) * du2TopOfEdge(k,edgeIndex)
end do
end do
do iCell = 1,nCells
Modified: branches/ocean_projects/option1_b4b_test/src/framework/mpas_block_creator.F
===================================================================
--- branches/ocean_projects/option1_b4b_test/src/framework/mpas_block_creator.F 2012-08-14 21:29:46 UTC (rev 2100)
+++ branches/ocean_projects/option1_b4b_test/src/framework/mpas_block_creator.F 2012-08-14 21:33:46 UTC (rev 2101)
@@ -1136,20 +1136,25 @@
cellIDSorted(2,i) = i
end do
call mpas_quicksort(block_ptr % mesh % nCells, cellIDSorted)
+ block_ptr % mesh % cellPermute % array(1:block_ptr % mesh % nCells) = cellIDSorted(2,:)
+ block_ptr % mesh % cellPermute % array(block_ptr % mesh % nCells + 1) = block_ptr % mesh % nCells + 1
do i=1,block_ptr % mesh % nEdges
edgeIDSorted(1,i) = block_ptr % mesh % indexToEdgeID % array(i)
edgeIDSorted(2,i) = i
end do
call mpas_quicksort(block_ptr % mesh % nEdges, edgeIDSorted)
+ block_ptr % mesh % edgePermute % array(1:block_ptr % mesh % nEdges) = edgeIDSorted(2,:)
+ block_ptr % mesh % edgePermute % array(block_ptr % mesh % nEdges + 1) = block_ptr % mesh % nEdges + 1
do i=1,block_ptr % mesh % nVertices
vertexIDSorted(1,i) = block_ptr % mesh % indexToVertexID % array(i)
vertexIDSorted(2,i) = i
end do
call mpas_quicksort(block_ptr % mesh % nVertices, vertexIDSorted)
+ block_ptr % mesh % vertexPermute % array(1:block_ptr % mesh % nVertices) = vertexIDSorted(2,:)
+ block_ptr % mesh % vertexPermute % array(block_ptr % mesh % nVertices + 1) = block_ptr % mesh % nVertices + 1
-
do i=1,block_ptr % mesh % nCells
do j=1,block_ptr % mesh % nEdgesOnCell % array(i)
k = mpas_binary_search(cellIDSorted, 2, 1, block_ptr % mesh % nCells, &
</font>
</pre>