<p><b>mpetersen@lanl.gov</b> 2010-10-29 14:18:00 -0600 (Fri, 29 Oct 2010)</p><p>BRANCH COMMIT, topography branch<br>
<br>
Removing unnecessary if statements from time_integration module.<br>
See attached emails for description. <br>
<br>
Date: Thu, 28 Oct 2010 15:48:25 -0600 (MDT)<br>
From: Mark Petersen <mpetersen@lanl.gov><br>
To: MPAS developers list <mpas-developers@ucar.edu><br>
Subject: [mpas-developers] nVertices+1<br>
<br>
Michael and others,<br>
<br>
In the sw and ocean core, there is a section (pasted below) with the line<br>
if (verticesOnEdge(1,iEdge) <= nVertices) then<br>
<br>
I understand that a cell numbered nCells+1 is outside of the local domain.<br>
But all vertices are always within the local domain, right? That is, I<br>
think verticesOnEdge <= nVertices is always true.<br>
<br>
If you agree, I will remove the if lines. Actually, this was already done<br>
in core_hyd_atmos, so I can just copy it over to sw, and make similar<br>
changes to the ocean core.<br>
<br>
Date: Thu, 28 Oct 2010 16:48:37 -0600<br>
From: Michael Duda <duda@ucar.edu><br>
To: Mark Petersen <mpetersen@lanl.gov><br>
Cc: MPAS developers list <mpas-developers@ucar.edu><br>
Subject: Re: [mpas-developers] nVertices+1<br>
<br>
Hi, Mark.<br>
<br>
I agree -- it should never be the case that we have an edge in the local<br>
domain (owned + ghost edges) but not both of its endpoints; in other<br>
words, verticesOnEdge(:,1:nEdges) should never be greater than<br>
nVertices. Please feel free to remove these if-tests whenever it's<br>
convenient!<br>
<br>
Cheers,<br>
Michael<br>
<br>
<br>
Date: Thu, 28 Oct 2010 16:23:33 -0600 (MDT)<br>
From: Mark Petersen <mpetersen@lanl.gov><br>
To: MPAS developers list <mpas-developers@ucar.edu><br>
Subject: [mpas-developers] two more if statements<br>
<br>
<br>
Here are two more places to remove if statements, where they have remained<br>
in the sw and ocean core but have already been removed from the hyd_atm<br>
core. Again, I think they are unnecessary. If everyone agrees, I will<br>
update the sw and ocean cores<br>
<br>
Date: Thu, 28 Oct 2010 16:55:47 -0600<br>
From: Michael Duda <duda@ucar.edu><br>
To: Mark Petersen <mpetersen@lanl.gov><br>
Cc: MPAS developers list <mpas-developers@ucar.edu><br>
Subject: Re: [mpas-developers] two more if statements<br>
<br>
Hi, Mark.<br>
<br>
I also agree that the tests in both of the sections you've identified<br>
below should also be safe to remove. In both cases below, iEdge or eoe<br>
may indeed reference an edge outside the local block, but this reference<br>
would have been switched to the index of the "garbage edge" during<br>
domain setup, and we expect that the values from the garbage edge should<br>
never influence the values of owned cells/edges/vertices. Please feel free<br>
to eliminate these if-tests, too, and thanks for taking such a careful<br>
look through the code!<br>
<br>
Cheers,<br>
Michael<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-10-29 18:52:14 UTC (rev 587)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/Registry 2010-10-29 20:18:00 UTC (rev 588)
@@ -105,7 +105,7 @@
var persistent real h_s ( nCells ) 0 iro h_s mesh - -
# Space needed for advection
-var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
+var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 - deriv_two mesh - -
var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
# !! NOTE: the following arrays are needed to allow the use
@@ -152,31 +152,33 @@
var persistent real tend_tracer2 ( nVertLevels nCells Time ) 1 - tracer2 tend tracers testing
# Diagnostic fields: only written to output
-var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
+var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
-var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
-var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
-var persistent real h_vertex ( nVertLevels nVertices Time ) 2 o h_vertex state - -
+var persistent real vorticity ( nVertLevels nVertices Time ) 2 - vorticity state - -
+var persistent real pv_edge ( nVertLevels nEdges Time ) 2 - pv_edge state - -
+var persistent real h_edge ( nVertLevels nEdges Time ) 2 - h_edge state - -
+var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real ke_edge ( nVertLevels nEdges Time ) 2 o ke_edge state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
+var persistent real ke_edge ( nVertLevels nEdges Time ) 2 - ke_edge state - -
+var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 - pv_vertex state - -
var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
-var persistent real zMid ( nVertLevels nCells Time ) 2 o zMid state - -
-var persistent real zTop ( nVertLevelsP1 nCells Time ) 2 o zTop state - -
-var persistent real zMidEdge ( nVertLevels nEdges Time ) 2 o zMidEdge state - -
-var persistent real zTopEdge ( nVertLevelsP1 nEdges Time ) 2 o zTopEdge state - -
+var persistent real zMid ( nVertLevels nCells Time ) 2 - zMid state - -
+var persistent real zTop ( nVertLevelsP1 nCells Time ) 2 - zTop state - -
+var persistent real zMidEdge ( nVertLevels nEdges Time ) 2 - zMidEdge state - -
+var persistent real zTopEdge ( nVertLevelsP1 nEdges Time ) 2 - zTopEdge state - -
var persistent real p ( nVertLevels nCells Time ) 2 o p state - -
-var persistent real pTop ( nVertLevelsP1 nCells Time ) 2 o pTop state - -
-var persistent real pZLevel ( nVertLevels nCells Time ) 2 o pZLevel state - -
+var persistent real pTop ( nVertLevelsP1 nCells Time ) 2 - pTop state - -
+var persistent real pZLevel ( nVertLevels nCells Time ) 2 - pZLevel state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
var persistent real ssh ( nCells Time ) 2 o ssh state - -
+var persistent real pgrad ( nVertLevels nCells Time ) 2 o pgrad state - -
+var persistent real pgrad_edge ( nVertLevels nEdges Time ) 2 o pgrad_edge state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-10-29 18:52:14 UTC (rev 587)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_test_cases.F 2010-10-29 20:18:00 UTC (rev 588)
@@ -121,9 +121,9 @@
nVertLevels = block_ptr % mesh % nVertLevels
!mrp temp, to agree with previous test
-h(1,:) = 500.0
-h(2,:) = 1250.0
-h(3,:) = 3250.0
+!h(1,:) = 500.0
+!h(2,:) = 1250.0
+!h(3,:) = 3250.0
pi=3.1415
! Tracer blob in Central Pacific, away from boundaries:
@@ -162,12 +162,12 @@
! tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
! for x3, 25 layer test
- tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 10.0 ! temperature
- tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+ !tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 10.0 ! temperature
+ !tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
- (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+ !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
+ !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
+ ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
!tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = block_ptr % mesh % lonCell % array(iCell)
!tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = block_ptr % mesh % latCell % array(iCell)
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-10-29 18:52:14 UTC (rev 587)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/module_time_integration.F 2010-10-29 20:18:00 UTC (rev 588)
@@ -361,33 +361,62 @@
! for explanation of divergence operator.
!
! for z-level, only compute height tendency for top layer.
+
if (config_vert_grid_type.eq.'isopycnal') then
kmax = nVertLevels
elseif (config_vert_grid_type.eq.'zlevel') then
kmax = 1
endif
+
+ ! mrp efficiency note: I try to remove if statements from within
+ ! loops, but did not here. I could use a loop like
+ ! do k=1,min(1,maxLevelEdgeTop(iedge))
+ ! for z-level, but then cells in the outer halo ring do not get
+ ! the proper flux, as they do in the current code.
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,kmax
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,kmax
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell2) = tend_h(k,cell2) + flux
- end do
- end if
- end do
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,kmax
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
+ end do
+!mrp old begin
+if (1==2) then
+ ! mrp efficiency note: I try to remove if statements from within
+ ! loops, but did not here. I could use a loop like
+ ! do k=1,min(1,maxLevelEdgeTop(iedge))
+ ! for z-level, but then cells in the outer halo ring do not get
+ ! the proper flux, as they do in the current code.
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells) then
+ do k=1,kmax
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ end do
+ end if
+ if (cell2 <= nCells) then
+ do k=1,kmax
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
+ end if
+ end do
+endif
+!mrp old end
+
do iCell=1,nCells
do k=1,kmax
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
end do
+ !print *, 'nEdges,nEdgesSolve,nCells,size(cellsOnEdge),maxval(cellsOnEdge(1,1:nEdges)),maxval(cellsOnEdge(1,1:nEdgesSolve))', &
+ !nEdges,nEdgesSolve,nCells,size(cellsOnEdge),maxval(cellsOnEdge(1,1:nEdges)),maxval(cellsOnEdge(1,1:nEdgesSolve))
+
!
! height tendency: vertical advection term -d/dz(hw)
!
@@ -771,16 +800,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
- do iTracer=1,num_tracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
end do
- end if
+ end do
end do
else if (config_tracer_adv_order == 3) then
@@ -789,56 +816,52 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ !-- else u <= 0:
+ else
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
else if (config_tracer_adv_order == 4) then
@@ -847,46 +870,42 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
endif ! if (config_tracer_adv_order == 2 )
@@ -1126,7 +1145,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, rho0Inv
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
@@ -1136,7 +1155,7 @@
hZLevel, ssh, zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, ke_edge, &
+ circulation, vorticity, ke, ke_edge, pgrad, pgrad_edge, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -1180,6 +1199,8 @@
pTop => s % pTop % array
MontPot => s % MontPot % array
ssh => s % ssh % array
+ pgrad => s % pgrad % array
+ pgrad_edge => s % pgrad_edge % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1231,11 +1252,9 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
end do
else if (config_thickness_adv_order == 3) then
@@ -1244,50 +1263,46 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ 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, grid % nEdgesOnCell % array (cell1)
- 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, grid % nEdgesOnCell % array (cell2)
+ 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, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ !-- else u <= 0:
+ else
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else u <= 0:
- else
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
else if (config_thickness_adv_order == 4) then
@@ -1296,40 +1311,36 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ 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, grid % nEdgesOnCell % array (cell1)
- 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, grid % nEdgesOnCell % array (cell2)
+ 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, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ 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 if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
endif ! if(config_thickness_adv_order == 2)
@@ -1347,18 +1358,10 @@
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- if (vertex1 <= nVertices) then
- do k=1,maxLevelVertexBot(vertex1) !mrp remove if for efficiency
- circulation(k,vertex1) = circulation(k,vertex1) &
- - dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
- if (vertex2 <= nVertices) then
- do k=1,maxLevelVertexBot(vertex2) !mrp remove if for efficiency
- circulation(k,vertex2) = circulation(k,vertex2) &
- + dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
+ do k=1,maxLevelEdgetop(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)
@@ -1379,10 +1382,10 @@
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
+ r = 1.0 / areaCell(iCell)
+ do k = 1,maxLevelCell(iCell)
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
enddo
!
@@ -1408,11 +1411,9 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe <= nEdges) then
- do k = 1,maxLevelEdgeBot(iEdge)
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end if
+ do k = 1,maxLevelEdgeBot(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
end do
end do
@@ -1425,11 +1426,9 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
- ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
end do
!
@@ -1457,16 +1456,14 @@
!
pv_cell(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if (iCell <= nCells) then
- do k = 1,maxLevelCell(iCell) !mrp remove if for efficiency
- pv_cell(k,iCell) = pv_cell(k,iCell) &
- + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
- / areaCell(iCell)
- enddo
- endif
- enddo
+ do i=1,vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ do k = 1,maxLevelCell(iCell)
+ pv_cell(k,iCell) = pv_cell(k,iCell) &
+ + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
+ / areaCell(iCell)
+ enddo
+ enddo
enddo
!
@@ -1483,17 +1480,15 @@
!
! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
+ ! ( 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)
- if(iEdge <= nEdges) then
+ do i=1,vertexDegree
+ iEdge = edgesOnVertex(i,iVertex)
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
- endif
end do
end do
@@ -1512,13 +1507,11 @@
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
+ do k=1,maxLevelEdgeTop(iEdge) !mrp remove if for efficiency
gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
- pv_cell(k,cellsOnEdge(1,iEdge))) &
/ dcEdge(iEdge)
- enddo
- endif
+ enddo
enddo
!
@@ -1725,7 +1718,38 @@
endif
+ !
+ ! Compute pressure gradient in each cell
+ !
+ ! This is for output and testing only. We could eventually
+ ! remove the variables pgrad and pgrad_edge.
+ rho0Inv = 1.0/config_rho0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ pgrad_edge(k,iEdge) = - rho0Inv*( pZLevel(k,cell2) &
+ - pZLevel(k,cell1) )/dcEdge(iEdge)
+ end do
+ end do
+ pgrad(:,:) = 0.0
+ do iCell=1,nCells
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ do k=1,maxLevelEdgeTop(iEdge)
+ ! mrp first attempt, I think it is wrong:
+ !pgrad(k,iCell) = pgrad(k,iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * abs(pgrad_edge(k,iEdge))
+ pgrad(k,iCell) = pgrad(k,iCell) + abs(pgrad_edge(k,iEdge))
+ end do
+ end do
+ do k=1,maxLevelCell(iCell)
+ ! mrp first attempt, I think it is wrong:
+! pgrad(k,iCell) = pgrad(k,iCell) / areaCell(iCell)
+ pgrad(k,iCell) = pgrad(k,iCell) / nEdgesOnCell(iCell)
+ end do
+ end do
+
end subroutine compute_solve_diagnostics
Modified: branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-10-29 18:52:14 UTC (rev 587)
+++ branches/ocean_projects/topography2_mrp/src/core_ocean/mpas_interface.F 2010-10-29 20:18:00 UTC (rev 588)
@@ -62,7 +62,6 @@
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
maxLevelCell => block % mesh % maxLevelCell % array
- maxLevelCell => block % mesh % maxLevelCell % array
maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
@@ -80,20 +79,22 @@
if (config_vert_grid_type.eq.'zlevel') then
- ! hZLevel should be in the grid.nc and restart.nc file
+ ! hZLevel should be in the grid.nc and restart.nc file
+
+ zTopZLevel(1) = 0.0
+ do iLevel = 1,nVertLevels
+ zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
+ zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
+ end do
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
- enddo
-
endif
+ ! for z-grids, maxLevelCell should be in input state
! Isopycnal grid uses all vertical cells
if (config_vert_grid_type.eq.'isopycnal') then
- maxLevelCell = nVertLevels
+ maxLevelCell(1:nCells) = nVertLevels
endif
+ maxLevelCell(nCells+1) = 0
! mrp insert topography mesa for testing
! do iCell = 1,nCells
@@ -103,28 +104,28 @@
! lonCell(iCell)<210.0/180.*3.14) then
! maxLevelCell(iCell) = 10
! endif
- ! enddo
+ ! end do
! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- maxLevelEdgeTop(iEdge) = &
- min( maxLevelCell(cell1), &
- maxLevelCell(cell2) )
- else
- maxLevelEdgeTop(iEdge) = 0
- endif
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cell1), &
+ maxLevelCell(cell2) )
+ else
+ maxLevelEdgeTop(iEdge) = 0
+ endif
end do
+ maxLevelEdgeTop(nEdges+1) = 0
! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
-
if (cell1 <= nCells .and. cell2 <= nCells) then
maxLevelEdgeBot(iEdge) = &
max( maxLevelCell(cell1), &
@@ -132,8 +133,8 @@
else
maxLevelEdgeBot(iEdge) = 0
endif
-
end do
+ maxLevelEdgeBot(nEdges+1) = 0
! set boundary edge
boundaryEdge=1
@@ -153,34 +154,34 @@
boundaryCell(k,cell1) = 1
boundaryCell(k,cell2) = 1
endif
- enddo
- enddo
+ end do
+ end do
! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
do iVertex = 1,nVertices
- maxLevelVertexBot(iVertex) = 0
- do i=1,vertexDegree
- cell = cellsOnVertex(i,iVertex)
- if (cell <= nCells) then
- maxLevelVertexBot(iVertex) = &
- max( maxLevelVertexBot(iVertex), &
- maxLevelCell(cell) )
- endif
- end do
+ maxLevelVertexBot(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
end do
maxLevelVertexBot(nVertices+1) = 0
! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
do iVertex = 1,nVertices
- maxLevelVertexTop(iVertex) = 0
- do i=1,vertexDegree
- cell = cellsOnVertex(i,iVertex)
- if (cell <= nCells) then
- maxLevelVertexTop(iVertex) = &
- min( maxLevelVertexTop(iVertex), &
- maxLevelCell(cell) )
- endif
- end do
+ maxLevelVertexTop(iVertex) = 0
+ do i=1,vertexDegree
+ cell = cellsOnVertex(i,iVertex)
+ if (cell <= nCells) then
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cell) )
+ endif
+ end do
end do
maxLevelVertexTop(nVertices+1) = 0
@@ -196,27 +197,29 @@
end do
! update halos for maxLevel variables
- block => domain % blocklist
- do while (associated(block))
+! mrp 101028 actually, dont update halos. I want maxLevel* on the
+! outside edge of a halo to be zero, so we do not loop over those.
+! block => domain % blocklist
+! do while (associated(block))
- call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- block % mesh % maxLevelCell % array, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- block % mesh % maxLevelEdgeTop % array, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- block % mesh % maxLevelEdgeBot % array, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- block % mesh % maxLevelVertexTop % array, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- call dmpar_exch_halo_field1dInteger(domain % dminfo, &
- block % mesh % maxLevelVertexBot % array, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+! call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+! block % mesh % maxLevelCell % array, block % mesh % nCells, &
+! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+! block % mesh % maxLevelEdgeTop % array, block % mesh % nEdges, &
+! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+! block % mesh % maxLevelEdgeBot % array, block % mesh % nEdges, &
+! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+! call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+! block % mesh % maxLevelVertexTop % array, block % mesh % nVertices, &
+! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+! call dmpar_exch_halo_field1dInteger(domain % dminfo, &
+! block % mesh % maxLevelVertexBot % array, block % mesh % nVertices, &
+! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- block => block % next
- end do
+! block => block % next
+! end do
end subroutine mpas_init_domain
</font>
</pre>