<p><b>dwj07@fsu.edu</b> 2011-10-25 12:05:49 -0600 (Tue, 25 Oct 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Some simple optimizations are applied to the restoring tendency computation, as well as ocn_diagnostic_solve.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_restoring.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_restoring.F 2011-10-25 17:48:47 UTC (rev 1131)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_restoring.F 2011-10-25 18:05:49 UTC (rev 1132)
@@ -111,6 +111,7 @@
integer :: iCell, nCellsSolve, k
real (kind=RKIND), dimension(:), pointer :: temperatureRestore, salinityRestore
+ real (kind=RKIND) :: invTemp, invSalinity
err = 0
@@ -121,17 +122,14 @@
temperatureRestore => grid % temperatureRestore % array
salinityRestore => grid % salinityRestore % array
+ invTemp = 1.0 / (temperatureTimeScale * 86400.0)
+ invSalinity = 1.0 / (salinityTimeScale * 86400.0)
+
k = 1 ! restoring only in top layer
do iCell=1,nCellsSolve
+ tend(indexT, k, iCell) = tend(indexT, k, iCell) - h(k,iCell)*(tracers(indexT, k, iCell) - temperatureRestore(iCell)) * invTemp
+ tend(indexS, k, iCell) = tend(indexS, k, iCell) - h(k,iCell)*(tracers(indexS, k, iCell) - salinityRestore(iCell)) * invSalinity
- tend(indexT, k, iCell) = tend(indexT, k, iCell) &
- - h(k,iCell)*(tracers(indexT, k, iCell) - temperatureRestore(iCell)) &
- / (temperatureTimeScale * 86400.0)
-
- tend(indexS, k, iCell) = tend(indexS, k, iCell) &
- - h(k,iCell)*(tracers(indexS, k, iCell) - salinityRestore(iCell)) &
- / (salinityTimeScale * 86400.0)
-
! write(6,10) iCell, tracers(indexT, k, iCell), &
! temperatureRestore(iCell), tracers(indexT, k, iCell), &
! (tracers(indexT, k, iCell) - temperatureRestore(iCell)) &
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tendency.F 2011-10-25 17:48:47 UTC (rev 1131)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tendency.F 2011-10-25 18:05:49 UTC (rev 1132)
@@ -620,6 +620,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
+ real (kind=RKIND) :: r_tmp, invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, invLength
h => s % h % array
u => s % u % array
@@ -796,66 +797,45 @@
tracers(s % index_temperature,:,nCells+1) = -1e34
tracers(s % index_salinity,:,nCells+1) = -1e34
- !
- ! Compute circulation and relative vorticity at each vertex
- !
circulation(:,:) = 0.0
+ vorticity(:,:) = 0.0
+ divergence(:,:) = 0.0
+ ke(:,:) = 0.0
+ v(:,:) = 0.0
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
- circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- do k=1,maxLevelVertexBot(iVertex)
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
- !
- ! Compute the divergence at each cell center
- !
- divergence(:,:) = 0.0
- do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,maxLevelCell(iCell)
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
- enddo
- !
- ! Compute kinetic energy in each cell
- !
- ke(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ invAreaTri1 = 1.0 / areaTriangle(vertex1)
+ invAreaTri2 = 1.0 / areaTriangle(vertex2)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k=1,maxLevelEdgeBot(iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- enddo
- end do
- do iCell = 1,nCells
- do k = 1,maxLevelCell(iCell)
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- enddo
- enddo
+ ! 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
- !
- ! Compute v (tangential) velocities
- !
- v(:,:) = 0.0
- do iEdge = 1,nEdges
+ 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 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)
eoe = edgesOnEdge(i,iEdge)
! mrp 101115 note: in order to include flux boundary conditions,
@@ -864,6 +844,7 @@
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
+
end do
!
@@ -908,59 +889,51 @@
end do
end do
- !
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
pv_cell(:,:) = 0.0
+ pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
do i=1,vertexDegree
iCell = cellsOnVertex(i,iVertex)
+ iEdge = edgesOnVertex(i,iVertex)
+
+ invAreaCell1 = 1.0 / areaCell(iCell)
+
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
do k = 1,maxLevelCell(iCell)
- pv_cell(k,iCell) = pv_cell(k,iCell) &
- + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
- / areaCell(iCell)
+ pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) * invAreaCell1
enddo
- enddo
- enddo
- !
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
- pv_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
+ ! Compute pv at the edges
+ ! ( this computes pv_edge at all edges bounding real cells )
do k=1,maxLevelEdgeBot(iEdge)
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
- end do
- end do
+ enddo
+ enddo
- !
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
gradPVn(:,:) = 0.0
+ gradPVt(:,:) = 0.0
do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnEdge(2, iEdge)
+ vertex1 = verticesOnedge(1, iEdge)
+ vertex2 = verticesOnedge(2, iEdge)
+
+ invLength = 1.0 / dcEdge(iEdge)
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradPVn for all edges bounding real cells )
do k=1,maxLevelEdgeTop(iEdge)
- gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
- - pv_cell(k,cellsOnEdge(1,iEdge))) &
- / dcEdge(iEdge)
+ gradPVn(k,iEdge) = (pv_cell(k,cell2) - pv_cell(k,cell1)) * invLength
enddo
- enddo
- !
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
- do iEdge = 1,nEdges
+ invLength = 1.0 / dvEdge(iEdge)
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
do k = 1,maxLevelEdgeBot(iEdge)
- gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
- - pv_vertex(k,verticesOnEdge(1,iEdge))) &
- /dvEdge(iEdge)
+ gradPVt(k,iEdge) = (pv_vertex(k,vertex2) - pv_vertex(k,vertex1)) * invLength
enddo
+
enddo
!
@@ -991,6 +964,7 @@
! Pressure
! This section must be after computing rho
!
+ ! dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
if (config_vert_grid_type.eq.'isopycnal') then
! For Isopycnal model.
</font>
</pre>