<p><b>dwj07@fsu.edu</b> 2012-05-31 13:28:35 -0600 (Thu, 31 May 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
General code clean up and commenting in block_creator. Added doxygen comments to each routine, and a header comment for the module. Minimal inline comments added as well.<br>
<br>
Code reorganization in mpas_io_input.<br>
<br>
Bug fix in mpas_dmpar.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/multiple_blocks/src/framework/mpas_block_creator.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/framework/mpas_block_creator.F 2012-05-30 21:41:48 UTC (rev 1950)
+++ branches/omp_blocks/multiple_blocks/src/framework/mpas_block_creator.F 2012-05-31 19:28:35 UTC (rev 1951)
@@ -1,3 +1,19 @@
+!***********************************************************************
+!
+! mpas_block_creator
+!
+!> \brief This module is responsible for the intial creation and setup of the block data structures.
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This module provides routines for the creation of blocks, with both an
+!> arbitrary number of blocks per processor and an arbitrary number of halos for
+!> each block. The provided routines also setup the exchange lists for each
+!> block.
+!
+!-----------------------------------------------------------------------
+
module mpas_block_creator
use mpas_dmpar
@@ -10,15 +26,29 @@
contains
- subroutine mpas_block_creator_setup_blocks_and_0halo_cells(domain, indexToCellIDField, cellList, blockID, blockStart, blockCount)!{{{
+!***********************************************************************
+!
+! routine mpas_block_creator_setup_blocks_and_0halo_cells
+!
+!> \brief Initializes the list of blocks, and determines 0 halo cell indices.
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This routine sets up the linked list of blocks, and creates the
+!> indexToCellID field for the 0 halo. The information required to setup these
+!> structures is provided as input in cellList, blockID, blockStart, and
+!> blockCount.
+!
+!-----------------------------------------------------------------------
- type (domain_type), pointer :: domain
- type (block_type), pointer :: block
- type (field1dInteger), pointer :: indexToCellIDField
- integer, dimension(:), intent(in) :: cellList
- integer, dimension(:), intent(in) :: blockID
- integer, dimension(:), intent(in) :: blockStart
- integer, dimension(:), intent(in) :: blockCount
+ subroutine mpas_block_creator_setup_blocks_and_0halo_cells(domain, indexToCellID, cellList, blockID, blockStart, blockCount)!{{{
+ type (domain_type), pointer :: domain !< Input: Domain information
+ type (field1dInteger), pointer :: indexToCellID !< Input/Output: indexToCellID field
+ integer, dimension(:), intent(in) :: cellList !< Input: List of cell indices owned by this processor
+ integer, dimension(:), intent(in) :: blockID !< Input: List of block indices owned by this processor
+ integer, dimension(:), intent(in) :: blockStart !< Input: Indices of starting cell id in cellList for each block
+ integer, dimension(:), intent(in) :: blockCount !< Input: Number of cells from cellList owned by each block.
integer :: nHalos
type (block_type), pointer :: blockCursor
@@ -30,47 +60,77 @@
nBlocks = size(blockID)
nHalos = config_num_halos
- allocate(domain % blocklist)
- block => domain % blocklist
- nullify(block % prev)
- nullify(block % next)
-
- allocate(indexToCellIDField)
- nullify(indexToCellIDField % next)
+ if(nBlocks > 0) then
+ ! Setup first block
+ allocate(domain % blocklist)
+! block => domain % blocklist
+ nullify(domain % blocklist % prev)
+ nullify(domain % blocklist % next)
+
+ ! Setup first block field
+ allocate(indexToCellID)
+ nullify(indexToCellID % next)
+
+ ! Loop over blocks
+ blockCursor => domain % blocklist
+ fieldCursor => indexToCellID
+ do i = 1, nBlocks
+ ! Initialize block information
+ blockCursor % blockID = blockID(i)
+ blockCursor % localBlockID = i - 1
+ blockCursor % domain => domain
+
+ ! Link to block, and setup array size
+ fieldCursor % block => blockCursor
+ fieldCursor % dimSizes(1) = blockCount(i)
+ nullify(fieldCursor % ioinfo)
+
+ ! Initialize exchange lists
+ call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % sendList, nHalos)
+ call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % recvList, nHalos)
+ call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % copyList, nHalos)
+
+ ! Allocate array, and copy indices into array
+ allocate(fieldCursor % array(fieldCursor % dimSizes(1)))
+ fieldCursor % array(1:fieldCursor % dimSizes(1)) = cellList(blockStart(i)+1:blockStart(i)+blockCount(i))
+
+ ! Advance cursors, and create new blocks as needed
+ if(i < nBlocks) then
+ allocate(blockCursor % next)
+ allocate(fieldCursor % next)
+
+ blockCursor % next % prev => blockCursor
+
+ blockCursor => blockCursor % next
+ fieldCursor => fieldCursor % next
+ end if
+
+ ! Nullify next pointers
+ nullify(blockCursor % next)
+ nullify(fieldCursor % next)
+ end do
+ else
+ ! If no blocks, nullify block and field pointers
+ nullify(domain % blocklist)
+ nullify(indexToCellID)
+ end if
+ end subroutine mpas_block_creator_setup_blocks_and_0halo_cells!}}}
- blockCursor => block
- fieldCursor => indexToCellIDField
- do i = 1, nBlocks
- blockCursor % blockID = blockID(i)
- blockCursor % localBlockID = i - 1
- blockCursor % domain => domain
-
- fieldCursor % block => blockCursor
- fieldCursor % dimSizes(1) = blockCount(i)
- nullify(fieldCursor % ioinfo)
+!***********************************************************************
+!
+! routine mpas_block_creator_build_0halo_cell_fields
+!
+!> \brief Initializes 0 halo cell based fields requried to work out halos
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This routine uses the previously setup 0 halo cell field, and the blocks of
+!> data read in by other routhers to determine all of the connectivity for the 0
+!> halo cell fields on all blocks on a processor.
+!
+!-----------------------------------------------------------------------
- call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % sendList, nHalos)
- call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % recvList, nHalos)
- call mpas_dmpar_init_mulithalo_exchange_list(fieldCursor % copyList, nHalos)
-
- allocate(fieldCursor % array(fieldCursor % dimSizes(1)))
- fieldCursor % array(1:fieldCursor % dimSizes(1)) = cellList(blockStart(i)+1:blockStart(i)+blockCount(i))
-
- if(i < nBlocks) then
- allocate(blockCursor % next)
- allocate(fieldCursor % next)
-
- blockCursor % next % prev => blockCursor
-
- blockCursor => blockCursor % next
- fieldCursor => fieldCursor % next
- end if
-
- nullify(blockCursor % next)
- nullify(fieldCursor % next)
- end do
- end subroutine mpas_block_creator_setup_blocks_and_0halo_cells!}}}
-
subroutine mpas_block_creator_build_0halo_cell_fields(indexToCellIDBlock, nEdgesOnCellBlock, cellsOnCellBlock, verticesOnCellBlock, edgesOnCellBlock, indexToCellID_0Halo, nEdgesOnCell_0Halo, cellsOnCell_0Halo, verticesOnCell_0Halo, edgesOnCell_0Halo)!{{{
type(field1dInteger), pointer :: indexToCellIDBlock !< Input: Block of read in indexToCellID field
type(field1dInteger), pointer :: nEdgesOnCellBlock !< Input: Block of read in nEdgesOnCell field
@@ -92,27 +152,39 @@
integer :: nCellsInBlock, maxEdges, nHalos
integer :: i, iHalo
- nHalos = size(indexToCellID_0Halo % sendList % halos)
+ nHalos = config_num_halos
+ ! Only sending from halo layer 1 for setup
allocate(sendingHaloLayers(1))
sendingHaloLayers(1) = 1
maxEdges = cellsOnCellBlock % dimSizes(1)
+ ! Build exchange list from the block of read in data to each block's index fields.
call mpas_dmpar_get_exch_list(1, indexToCellIDBlock, indexToCellID_0Halo)
- allocate(nEdgesOnCell_0Halo)
- nullify(nEdgesOncell_0Halo % next)
+ ! Setup header fields if at least 1 block exists
+ if(associated(indexToCellID_0Halo)) then
+ allocate(nEdgesOnCell_0Halo)
+ nullify(nEdgesOncell_0Halo % next)
- allocate(cellsOnCell_0Halo)
- nullify(cellsOnCell_0Halo % next)
+ allocate(cellsOnCell_0Halo)
+ nullify(cellsOnCell_0Halo % next)
+
+ allocate(verticesOnCell_0Halo)
+ nullify(verticesOnCell_0Halo % next)
+
+ allocate(edgesOnCell_0Halo)
+ nullify(edgesOnCell_0Halo % next)
+ else
+ nullify(nEdgesOnCell_0Halo)
+ nullify(cellsOnCell_0Halo)
+ nullify(verticesOnCell_0Halo)
+ nullify(edgesOnCell_0Halo)
+ end if
- allocate(verticesOnCell_0Halo)
- nullify(verticesOnCell_0Halo % next)
- allocate(edgesOnCell_0Halo)
- nullify(edgesOnCell_0Halo % next)
-
+ ! Loop over blocks
indexCursor => indexToCellID_0Halo
nEdgesCursor => nEdgesOnCell_0Halo
cellsOnCellCursor => cellsOnCell_0Halo
@@ -121,16 +193,19 @@
do while(associated(indexCursor))
nCellsInBlock = indexCursor % dimSizes(1)
+ ! Link to block structure
nEdgesCursor % block => indexCursor % block
cellsOnCellCursor % block => indexCursor % block
verticesOnCellCursor % block => indexCursor % block
edgesOnCellCursor % block => indexCursor % block
+ ! Nullify ioinfo, since this data is not read in
nullify(nEdgesCursor % ioinfo)
nullify(cellsOnCellCursor % ioinfo)
nullify(verticesOnCellCursor % ioinfo)
nullify(edgesOnCellCursor % ioinfo)
+ ! Setup array sizes
nEdgesCursor % dimSizes(1) = nCellsInBlock
cellsOnCellCursor % dimSizes(1) = maxEdges
cellsOnCellCursor % dimSizes(2) = nCellsInBlock
@@ -139,6 +214,7 @@
edgesOnCellCursor % dimSizes(1) = maxEdges
edgesOnCellCursor % dimSizes(2) = nCellsInBlock
+ ! Link exchange lists
nEdgesCursor % sendList => indexCursor % sendList
nEdgesCursor % recvList => indexCursor % recvList
nEdgesCursor % copyList => indexCursor % copyList
@@ -152,11 +228,13 @@
edgesOnCellCursor % recvList => indexCursor % recvList
edgesOnCellCursor % copyList => indexCursor % copyList
+ ! Allocate arrays
allocate(nEdgesCursor % array(nEdgesCursor % dimSizes(1)))
allocate(cellsOnCellCursor % array(cellsOnCellCursor % dimSizes(1), cellsOnCellCursor % dimSizes(2)))
allocate(verticesOnCellCursor % array(verticesOnCellCursor % dimSizes(1), verticesOnCellCursor % dimSizes(2)))
allocate(edgesOnCellCursor % array(edgesOnCellCursor % dimSizes(1), edgesOnCellCursor % dimSizes(2)))
+ ! Create new blocks and advance cursors as needed
indexCursor => indexCursor % next
if(associated(indexCursor)) then
allocate(nEdgesCursor % next)
@@ -171,18 +249,36 @@
end if
+ ! Nullify next pointers
nullify(nEdgesCursor % next)
nullify(cellsOnCellCursor % next)
nullify(verticesOnCellCursor % next)
nullify(edgesOnCellCursor % next)
- end do
+ end do ! indexCursor loop over blocks
+ ! Communicate data from read in blocks to each block's fields
call mpas_dmpar_alltoall_field(nEdgesOnCellBlock, nEdgesOnCell_0Halo, sendingHaloLayers)
call mpas_dmpar_alltoall_field(cellsOnCellBlock, cellsOnCell_0Halo, sendingHaloLayers)
call mpas_dmpar_alltoall_field(verticesOnCellBlock, verticesOnCell_0Halo, sendingHaloLayers)
call mpas_dmpar_alltoall_field(edgesOnCellBlock, edgesOnCell_0Halo, sendingHaloLayers)
end subroutine mpas_block_creator_build_0halo_cell_fields!}}}
+!***********************************************************************
+!
+! routine mpas_block_creator_build_0_and_1halo_vertex_fields
+!
+!> \brief Initializes 0 and 1 halo vertex based fields requried to work out halos
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This routine uses the previously setup 0 halo cell fields, and the blocks of
+!> data read in by other routhers to determine which vertices are in a blocks
+!> 0 and 1 halo for all blocks on a processor.
+!> NOTE: This routine can be used on either vertices or edges
+!
+!-----------------------------------------------------------------------
+
subroutine mpas_block_creator_build_0_and_1halo_vertex_fields(indexToVertexIDBlock, cellsOnVertexBlock, indexToCellID_0Halo, nEdgesOnCell_0Halo, verticesOnCell_0Halo, indexToVertexID_0Halo, cellsOnVertex_0Halo, nVerticesSolve)!{{{
type (field1dInteger), pointer :: indexToVertexIDBlock !< Input: indexToVertexID read in field
type (field2dInteger), pointer :: cellsOnVertexBlock !< Input: cellsOnVertex read in field
@@ -208,51 +304,71 @@
integer :: haloStart
integer :: iBlock, i, j, k
+ ! Setup sendingHaloLayers
allocate(sendingHaloLayers(1))
sendingHaloLayers(1) = 1
+ ! Get dimension information
maxEdges = verticesOnCell_0Halo % dimSizes(1)
vertexDegree = cellsOnVertexBlock % dimSizes(1)
nHalos = config_num_halos
- allocate(cellsOnVertex_0Halo)
- allocate(indexToVertexID_0Halo)
+ ! Setup initial block for each field
+ if(associated(indexToCellCursor)) then
+ allocate(cellsOnVertex_0Halo)
+ allocate(indexToVertexID_0Halo)
+ nullify(cellsOnVertex_0Halo % next)
+ nullify(indexToVertexID_0Halo % next)
+ else
+ nullify(cellsOnVertex_0Halo)
+ nullify(indexToVertexID_0Halo)
+ end if
+
+ ! Loop over blocks
indexToCellCursor => indexToCellID_0Halo
verticesOnCellCursor => verticesOnCell_0Halo
nEdgesCursor => nEdgesOnCell_0Halo
indexToVertexCursor => indexToVertexID_0Halo
cellsOnVertexCursor => cellsOnVertex_0Halo
do while(associated(indexToCellCursor))
+ ! Determine number of cells in block
nCellsInBlock = indexToCellCursor % dimSizes(1)
+ ! Determine all vertices in block
call mpas_block_decomp_all_edges_in_block(maxEdges, nCellsInBlock, nEdgesCursor % array, verticesOnCellCursor % array, nVerticesLocal, localVertexList)
+ ! Setup indexToVertex block
indexToVertexCursor % block => indexToCellCursor % block
nullify(indexToVertexCursor % ioinfo)
indexToVertexCursor % dimSizes(1) = nVerticesLocal
allocate(indexToVertexCursor % array(indexToVertexCursor % dimSizes(1)))
indexToVertexCursor % array(:) = localVertexList(:)
+ ! Setup cellsOnVertex block
cellsOnVertexCursor % block => indexToCellCursor % block
nullify(cellsOnVertexCursor % ioinfo)
cellsOnVertexCursor % dimSizes(1) = vertexDegree
cellsOnVertexCursor % dimSizes(2) = nVerticesLocal
allocate(cellsOnVertexCursor % array(cellsOnVertexCursor % dimSizes(1), cellsOnVertexCursor % dimSizes(2)))
+ ! Setup exchange lists
call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexCursor % sendList, nHalos+1)
call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexCursor % recvList, nHalos+1)
call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexCursor % copyList, nHalos+1)
+ ! Link exchange lists
cellsOnVertexCursor % sendList => indexToVertexCursor % sendList
cellsOnVertexCursor % recvList => indexToVertexCursor % recvList
cellsOnVertexCursor % copyList => indexToVertexCursor % copyList
+ ! Remove localVertexList array
deallocate(localVertexList)
+
+ ! Advance cursors, and create new blocks if needed
indexToCellCursor => indexToCellCursor % next
verticesOnCellCursor => verticesOnCellCursor % next
nEdgescursor => nEdgesCursor % next
-
if(associated(indexToCellCursor)) then
allocate(indexToVertexCursor % next)
indexToVertexCursor => indexToVertexCursor % next
@@ -260,85 +376,32 @@
allocate(cellsOnVertexCursor % next)
cellsOnVertexCursor => cellsOnVertexCursor % next
end if
+
+ ! Nullify next pointers
nullify(indexToVertexCursor % next)
nullify(cellsOnVertexCursor % next)
- end do
+ end do ! indexToCursor loop over blocks
+ ! Build exchangel ists from read in blocks to owned blocks.
call mpas_dmpar_get_exch_list(1, indexToVertexIDBlock, indexToVertexID_0Halo)
-! write(6,*) 'CELLSONVERTEXBLOCK'
-! indexToVertexCursor => indexToVertexIDBlock
-! do while(associated(indexToVertexCursor))
-! write(6,*) 'sendLists on block',indexToVertexCursor % block % blockID
-! exchListPtr => indexToVertexCursor % sendList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'to', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! write(6,*) 'recvLists on block',indexToVertexCursor % block % blockID
-! exchListPtr => indexToVertexCursor % recvList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'from', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! write(6,*) 'copyLists on block',indexToVertexCursor % block % blockID
-! exchListPtr => indexToVertexCursor % copyList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'to', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! indexToVertexCursor => indexToVertexCursor % next
-! end do
-
-! write(6,*) 'CELLSONVERTEX_0HALO'
-! cellsOnVertexCursor => cellsOnVertex_0Halo
-! do while(associated(cellsOnVertexCursor))
-! write(6,*) 'sendLists on block',cellsOnVertexCursor % block % blockID
-! exchListPtr => cellsOnVertexCursor % sendList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'to', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! write(6,*) 'recvLists on block',cellsOnVertexCursor % block % blockID
-! exchListPtr => cellsOnVertexCursor % recvList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'from', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! write(6,*) 'copyLists on block',cellsOnVertexCursor % block % blockID
-! exchListPtr => cellsOnVertexCursor % copyList % halos(1) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'to', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) i, exchListPtr % srcList(i), exchListPtr % destList(i)
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! cellsOnVertexCursor => cellsOnVertexCursor % next
-! end do
-
+ ! Perform all to all to get owned block data
call mpas_dmpar_alltoall_field(cellsOnVertexBlock, cellsOnVertex_0Halo, sendingHaloLayers)
- allocate(haloIndices)
- allocate(offSetField)
- allocate(vertexLimitField)
- allocate(nVerticesSolve)
+ ! Setup first block's fields if there is at least 1 block.
+ if(associated(indexToVertexID_0Halo)) then
+ allocate(haloIndices)
+ allocate(offSetField)
+ allocate(vertexLimitField)
+ allocate(nVerticesSolve)
+ else
+ nullify(haloIndices)
+ nullify(offSetField)
+ nullify(vertexLimitField)
+ nullify(nVerticesSolve)
+ end if
+ ! Loop over blocks
indexToVertexCursor => indexToVertexID_0Halo
cellsOnVertexCursor => cellsOnVertex_0Halo
indexToCellCursor => indexToCellID_0Halo
@@ -347,37 +410,45 @@
vertexLimitCursor => vertexLimitField
nVerticesSolveCursor => nVerticesSolve
do while(associated(indexToVertexCursor))
+ ! Determine 0 and 1 halo vertices
call mpas_block_decomp_partitioned_edge_list(indexToCellCursor % dimSizes(1), indexToCellCursor % array, &
vertexDegree, indexToVertexCursor % dimSizes(1), cellsOnVertexCursor % array, &
indexToVertexCursor % array, haloStart)
+ ! Link blocks
haloCursor % block => indexToVertexCursor % block
offSetCursor % block => indexToVertexCursor % block
vertexLimitCursor % block => indexToVertexCursor % block
nVerticesSolveCursor % block => indexToVertexCursor % block
+ ! Nullify io info
nullify(haloCursor % ioinfo)
nullify(offSetCursor % ioinfo)
nullify(vertexLimitCursor % ioinfo)
nullify(nVerticesSolveCursor % ioinfo)
+ ! Setup haloIndices
haloCursor % dimSizes(1) = indexToVertexCursor % dimSizes(1) - (haloStart-1)
allocate(haloCursor % array(haloCursor % dimSizes(1)))
haloCursor % array(:) = indexToVertexCursor % array(haloStart:indexToVertexCursor % dimSizes(1))
+ ! Link exchange lists
haloCursor % sendList => indexToVertexCursor % sendList
haloCursor % recvList => indexToVertexCursor % recvList
haloCursor % copyList => indexToVertexCursor % copyList
+ ! Determine offSet and limit on 0 halo vertices for exchange list creation
offSetCursor % scalar = haloStart - 1
vertexLimitCursor % scalar = haloStart - 1
- nVerticesSolveCursor % dimSizes(1) = nHalos+2 ! 1 for 0 halo, 1 for nHalo+1
+ ! Setup nVerticesSolve
+ nVerticesSolveCursor % dimSizes(1) = nHalos+2
allocate(nVerticesSolveCursor % array(nVerticesSolve % dimSizes(1)))
nVerticesSolveCursor % array = -1
nVerticesSolveCursor % array(1) = haloStart - 1
nVerticesSolveCursor % array(2) = indexToVertexCursor % dimSizes(1)
+ ! Advance cursors, and create new blocks if needed
indexToVertexCursor => indexToVertexCursor % next
cellsOnVertexCursor => cellsOnVertexCursor % next
indexToCellCursor => indexToCellCursor % next
@@ -394,14 +465,18 @@
allocate(nVerticesSolveCursor % next)
nVerticesSolveCursor => nVerticesSolveCursor % next
end if
+
+ ! Nullify next pointers
nullify(haloCursor % next)
nullify(offSetCursor % next)
nullify(vertexLimitCursor % next)
nullify(nVerticesSolveCursor % next)
end do
+ ! Create exchange lists from 0 halo to 1 halo vertices
call mpas_dmpar_get_exch_list(1, indexToVertexID_0Halo, haloIndices, offSetField, vertexLimitField)
+ ! Deallocate fields that are not needed anymore.
call mpas_deallocate_field(haloIndices)
call mpas_deallocate_field(offSetField)
call mpas_deallocate_field(vertexLimitCursor)
@@ -409,6 +484,22 @@
end subroutine mpas_block_creator_build_0_and_1halo_vertex_fields!}}}
+!***********************************************************************
+!
+! routine mpas_block_creator_build_cell_halos
+!
+!> \brief Builds cell halos
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This routine uses the previously setup 0 halo cell fields to determine
+!> which cells fall in each halo layer for a block. During this process, each
+!> halo's exchange lists are created. This process is performed for all blocks on
+!> a processor.
+!
+!-----------------------------------------------------------------------
+
subroutine mpas_block_creator_build_cell_halos(indexToCellID, nEdgesOnCell, cellsOnCell, verticesOnCell, edgesOnCell, nCellsSolve)!{{{
type (field1dInteger), pointer :: indexToCellID !< Input/Output: indexToCellID field for all halos
type (field1dInteger), pointer :: nEdgesOnCell !< Input/Output: nEdgesOnCell field for all halos
@@ -442,29 +533,46 @@
nHalos = config_num_halos
dminfo => indexToCellID % block % domain % dminfo
allocate(sendingHaloLayers(1))
- allocate(nCellsSolve)
- allocate(cellLimitField)
- allocate(offSetField)
+ ! Setup header fields
+ if(associated(indexToCellID)) then
+ allocate(nCellsSolve)
+ allocate(cellLimitField)
+ allocate(offSetField)
+
+ nullify(nCellsSolve % next)
+ nullify(cellLimitField % next)
+ nullify(offSetField % next)
+ else
+ nullify(nCellsSolve)
+ nullify(cellLimitField)
+ nullify(offSetField)
+ end if
+
+ ! Loop over blocks
offSetCursor => offsetField
cellLimitCursor => cellLimitField
indexCursor => indexToCellID
nCellsSolveCursor => nCellsSolve
do while (associated(indexCursor))
+ ! Setup offset
offSetCursor % scalar = indexCursor % dimSizes(1)
offSetCursor % block => indexCursor % block
nullify(offSetCursor % ioinfo)
+ ! Setup nCellsSolve
nCellsSolveCursor % dimSizes(1) = nHalos+1
allocate(nCellsSolveCursor % array(nCellsSolveCursor % dimSizes(1)))
nCellsSolveCursor % array(1) = indexCursor % dimSizes(1)
nCellsSolveCursor % block => indexCursor % block
nullify(nCellsSolveCursor % ioinfo)
+ ! Setup owned cellLimit
cellLimitCursor % scalar = indexCursor % dimSizes(1)
cellLimitCursor % block => indexCursor % block
nullify(cellLimitCursor % ioinfo)
+ ! Advance cursors and create new blocks if needed
indexCursor => indexCursor % next
if(associated(indexCursor)) then
allocate(offSetCursor % next)
@@ -476,16 +584,26 @@
allocate(cellLimitCursor % next)
cellLimitCursor => cellLimitCursor % next
end if
+
+ ! Nullify next pointers
nullify(offSetCursor % next)
nullify(nCellssolveCursor % next)
nullify(cellLimitCursor % next)
end do
+ ! Loop over halos
do iHalo = 1, nHalos
+ ! Sending halo layer is the current halo
sendingHaloLayers(1) = iHalo
- allocate(haloIndices)
+ if(associated(indexToCellID)) then
+ allocate(haloIndices)
+ nullify(haloIndices % next)
+ else
+ nullify(haloIndices)
+ end if
+ ! Loop over blocks
indexCursor => indexToCellID
nEdgesCursor => nEdgesOnCell
cellsOnCellCursor => cellsOnCell
@@ -494,11 +612,14 @@
haloCursor => haloIndices
offSetCursor => offSetField
do while(associated(indexCursor))
+ ! Determine block dimensions
nCellsInBlock = indexCursor % dimSizes(1)
maxEdges = cellsOnCellCursor % dimSizes(1)
+ ! Setup offSet
offSetCursor % scalar = nCellsInBlock
+ ! Setup block graphs
allocate(blockGraphWithHalo)
allocate(blockGraph)
allocate(blockGraph % vertexID(nCellsInBlock))
@@ -510,14 +631,14 @@
blockGraph % maxDegree = maxEdges
blockGraph % ghostStart = nCellsInBlock + 1
-!write(6,*) 'max nedges on block', indexCursor % block % blockid, maxval(nEdgescursor % array)
-
blockGraph % vertexID(:) = indexCursor % array(:)
blockGraph % nAdjacent(:) = nEdgesCursor % array(:)
blockGraph % adjacencyList(:,:) = cellsOnCellCursor % array(:,:)
+ ! Determine all cell id's with the next halo added
call mpas_block_decomp_add_halo(dminfo, blockGraph, blockGraphWithHalo)
+ ! Setup haloIndices
haloCursor % dimSizes(1) = blockGraphWithHalo % nVerticesTotal - blockGraphWithHalo % nVertices
allocate(haloCursor % array(haloCursor % dimSizes(1)))
haloCursor % array(:) = -1
@@ -528,6 +649,7 @@
haloCursor % block => indexCursor % block
nullify(haloCursor % ioinfo)
+ ! Deallocate block graphs
deallocate(blockGraphWithHalo % vertexID)
deallocate(blockGraphWithHalo % nAdjacent)
deallocate(blockGraphWithHalo % adjacencyList)
@@ -538,6 +660,7 @@
deallocate(blockGraph % adjacencyList)
deallocate(blockGraph)
+ ! Advance cursors and create new block if needed
indexCursor => indexCursor % next
nEdgesCursor => nEdgesCursor % next
cellsOnCellCursor => cellsOnCellCursor % next
@@ -548,55 +671,14 @@
allocate(haloCursor % next)
haloCursor => haloCursor % next
end if
+ ! Nullify next pointer
nullify(haloCursor % next)
end do ! indexCursor loop over blocks
+ ! Create exchange lists for current halo layer
call mpas_dmpar_get_exch_list(iHalo, indexToCellID, haloIndices, offSetField, cellLimitField)
- ! dwj: 05/24/12 debugging
-! write(6,*) 'new send lists'
-! indexCursor => indexToCellID
-! do while(associated(indexCursor))
-! exchListPtr => indexCursor % sendList % halos(iHalo) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'send list from block ', indexCursor % block % localBlockID, ' to proc ', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) exchListPtr % srcList(i), exchListPtr % destList(i)!, indexCursor % array(exchListPtr % srcList(i))
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! indexCursor => indexCursor % next
-! end do
-
-! write(6,*) 'new recv lists'
-! indexCursor => indexToCellID
-! do while(associated(indexCursor))
-! exchListPtr => indexCursor % recvList % halos(iHalo) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'recv list to block ', indexCursor % block % localBlockID, ' from proc ', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) exchListPtr % srcList(i), exchListPtr % destList(i), indexCursor % array(exchListPtr % srcList(i))
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! indexCursor => indexCursor % next
-! end do
-
-
-! write(6,*) 'new copy lists'
-! indexCursor => indexToCellID
-! do while(associated(indexCursor))
-! exchListPtr => indexCursor % copyList % halos(iHalo) % exchList
-! do while(associated(exchListPtr))
-! write(6,*) 'copy list from block ', indexCursor % block % localBlockID, ' to ', exchListPtr % endPointID
-! do i = 1, exchListPtr % nList
-! write(6,*) exchListPtr % srcList(i), exchListPtr % destList(i), indexCursor % array(exchListPtr % srcList(i))
-! end do
-! exchListPtr => exchListPtr % next
-! end do
-! indexCursor => indexCursor % next
-! end do
-
+ ! Loop over blocks
indexCursor => indexToCellID
nEdgesCursor => nEdgesOnCell
cellsOnCellCursor => cellsOnCell
@@ -605,11 +687,14 @@
haloCursor => haloIndices
nCellsSolveCursor => nCellsSolve
do while(associated(indexCursor))
+ ! Determine block dimensions
nCellsInBlock = indexCursor % dimSizes(1)
nCellsInHalo = haloCursor % dimSizes(1)
+ ! Setup new layer's nCellsSolve
nCellsSolveCursor % array(iHalo+1) = nCellsInBlock + nCellsInHalo
+ ! Copy cell indices into indexToCellID field
field1dArrayHolder => indexCursor % array
indexCursor % dimSizes(1) = nCellsSolveCursor % array(iHalo+1)
allocate(indexCursor % array(indexCursor % dimSizes(1)))
@@ -618,6 +703,7 @@
indexCursor % array(nCellsInBlock+1:nCellsSolveCursor % array(iHalo+1)) = haloCursor % array(1:nCellsInHalo)
deallocate(field1dArrayHolder)
+ ! Allocate space in nEdgesOnCell
field1dArrayHolder => nEdgesCursor % array
nEdgesCursor % dimSizes(1) = nCellsSolveCursor % array(iHalo+1)
allocate(nEdgesCursor % array(nEdgesCursor % dimSizes(1)))
@@ -625,6 +711,7 @@
nEdgesCursor % array(1:nCellsInBlock) = field1dArrayHolder(:)
deallocate(field1dArrayHolder)
+ ! Allocate space in cellsOnCell
field2dArrayHolder => cellsOnCellCursor % array
cellsOnCellCursor % dimSizes(2) = nCellsSolveCursor % array(iHalo+1)
allocate(cellsOnCellCursor % array(cellsOnCellCursor % dimSizes(1), cellsOnCellCursor % dimSizes(2)))
@@ -632,6 +719,7 @@
cellsOnCellCursor % array(:,1:nCellsInBlock) = field2dArrayHolder(:,:)
deallocate(field2dArrayHolder)
+ ! Allocate space in verticesOnCell
field2dArrayHolder => verticesOnCellCursor % array
verticesOnCellCursor % dimSizes(2) = nCellsSolveCursor % array(iHalo+1)
allocate(verticesOnCellCursor % array(verticesOnCellCursor % dimSizes(1), verticesOnCellCursor % dimSizes(2)))
@@ -639,6 +727,7 @@
verticesOnCellCursor % array(:,1:nCellsInBlock) = field2dArrayHolder(:,:)
deallocate(field2dArrayHolder)
+ ! Allocate space in edgesOnCell
field2dArrayHolder => edgesOnCellCursor % array
edgesOnCellCursor % dimSizes(2) = nCellsSolveCursor % array(iHalo+1)
allocate(edgesOnCellCursor % array(edgesOnCellCursor % dimSizes(1), edgesOnCellCursor % dimSizes(2)))
@@ -655,49 +744,40 @@
nCellsSolveCursor => nCellsSolveCursor % next
end do
+ ! Perform allToAll communications
call mpas_dmpar_alltoall_field(indexToCellID, indexToCellID, sendingHaloLayers)
call mpas_dmpar_alltoall_field(nEdgesOnCell, nEdgesOncell, sendingHaloLayers)
call mpas_dmpar_alltoall_field(cellsOnCell, cellsOnCell, sendingHaloLayers)
call mpas_dmpar_alltoall_field(verticesOnCell, verticesOnCell, sendingHaloLayers)
call mpas_dmpar_alltoall_field(edgesOnCell, edgesOnCell, sendingHaloLayers)
- ! dwj: 05/24/12 debugging
-! indexCursor => indexToCellID
-! do while(associated(indexCursor))
-!write(6,*) 'indexToCellID with HALO', iHalo,' on block ', indexCursor % block % blockID
-! do i = 1, indexcursor % dimSizes(1)
-! write(6,*) i, indexCursor % array(i)
-! end do
-! indexCursor => indexCursor % next
-! end do
-!
-! nEdgesCursor => nEdgesOnCell
-! do while(associated(nEdgesCursor))
-!write(6,*) 'nEdgesOnCell with HALO', iHalo, ' on block', nEdgescursor % block % blockID
-! do i = 1, nEdgesCursor % dimSizes(1)
-! write(6,*) i, nEdgesCursor % array(i)
-! end do
-! nEdgesCursor => nEdgesCursor % next
-! end do
-!
-! cellsOnCellCursor => cellsOnCell
-! do while(associated(cellsOnCellCursor))
-!write(6,*) 'cellsOnCell with HALO', iHalo,' on block ', cellsOnCellCursor % block % blockID
-! do i = 1, cellsOnCellCursor % dimSizes(2)
-! write(6,*) i, cellsOnCellCursor % array(:,i)
-! end do
-! cellsOnCellCursor => cellsOnCellCursor % next
-! end do
-
-
+ ! Deallocate haloindices field
call mpas_deallocate_field(haloIndices)
end do ! iHalo loop over nHalos
+ ! Deallocate array and field.
deallocate(sendingHaloLayers)
call mpas_deallocate_field(offSetField)
end subroutine mpas_block_creator_build_cell_halos!}}}
+!***********************************************************************
+!
+! routine mpas_block_creator_build_vertex_halos
+!
+!> \brief Builds vertex halos
+!> \author Doug Jacobsen
+!> \date 05/31/12
+!> \version SVN:$Id$
+!> \details
+!> This routine uses the previously setup 0 and 1 vertex fields and 0 halo cell fields to determine
+!> which vertices fall in each halo layer for a block. During this process, each
+!> halo's exchange lists are created. This process is performed for all blocks on
+!> a processor.
+!> NOTE: This routine can be used on either vertices or edges
+!
+!-----------------------------------------------------------------------
+
subroutine mpas_block_creator_build_vertex_halos(indexToCellID, nEdgesOnCell, nCellsSolve, verticesOnCell, indexToVertexID, cellsOnVertex, nVerticesSolve)!{{{
type (field1dInteger), pointer :: indexToCellID !< Input: indexToCellID field for all halos
type (field1dInteger), pointer :: nEdgesOnCell !< Input: nEdgesOnCell field for all halos
@@ -731,11 +811,22 @@
! Allocate some needed arrays and fields
allocate(sendingHaloLayers(1))
- allocate(haloIndices)
- allocate(offSetField)
- allocate(vertexLimitField)
+ if(associated(indexToVertexID)) then
+ allocate(haloIndices)
+ allocate(offSetField)
+ allocate(vertexLimitField)
+ nullify(haloIndices % next)
+ nullify(offSetField % next)
+ nullify(vertexLimitField % next)
+ else
+ nullify(haloIndices)
+ nullify(offSetField)
+ nullify(vertexLimitField)
+ end if
+
! Determine number of blocks, and setup field lists
+ ! Loop over blocks
nBlocks = 0
indexToVertexCursor => indexToVertexID
haloCursor => haloIndices
@@ -745,21 +836,26 @@
do while(associated(indexToVertexCursor))
nBlocks = nBlocks + 1
+ ! Setup vertexLimit and offSet
vertexLimitCursor % scalar = nVerticesSolveCursor % array(1)
offSetCursor % scalar = nVerticesSolveCursor % array(2)
+ ! Link blocks
vertexLimitCursor % block => indexToVertexCursor % block
offSetCursor % block => indexToVertexCursor % block
haloCursor % block => indexToVertexCursor % block
+ ! Nullify ioinfo
nullify(vertexLimitCursor % ioinfo)
nullify(offSetCursor % ioinfo)
nullify(haloCursor % ioinfo)
+ ! Link exchange lists
haloCursor % sendList => indexToVertexCursor % sendList
haloCursor % recvList => indexToVertexCursor % recvList
haloCursor % copyList => indexToVertexCursor % copyList
+ ! Advance cursors and create new blocks if needed
indexToVertexCursor => indexToVertexCursor % next
nVerticesSolveCursor => nVerticesSolveCursor % next
if(associated(indexToVertexCursor)) then
@@ -772,6 +868,8 @@
allocate(vertexLimitCursor % next)
vertexLimitCursor =>vertexLimitCursor % next
end if
+
+ ! Nullify next pointers
nullify(haloCursor % next)
nullify(offSetCursor % next)
nullify(vertexLimitCursor % next)
@@ -800,6 +898,8 @@
! Append new unique vertex id's to indexToVertexID field.
do iHalo = 3, nHalos+2
sendingHaloLayers(1) = iHalo-1
+
+ ! Loop over blocks
indexToVertexCursor => indexToVertexID
nEdgesCursor => nEdgesOnCell
nCellsSolveCursor => nCellsSolve
@@ -812,7 +912,7 @@
nCellsInBlock = nCellsSolveCursor % array(iHalo-1)
offSetCursor % scalar = nVerticesSolveCursor % array(iHalo-1)
-! call mpas_block_decomp_all_edges_in_block(maxEdges, nCellsInBlock, nEdgesCursor % array(1:nCellsInBlock), verticesOnCellCursor % array(:,1:nCellsInBlock), nVerticesLocal, localVertexList)
+ ! Determine all vertices in block
call mpas_block_decomp_all_edges_in_block(maxEdges, nCellsInBlock, nEdgesCursor % array, verticesOnCellCursor % array, nVerticesLocal, localVertexList)
nVerticesSolveCursor % array(iHalo) = nVerticesLocal
@@ -821,6 +921,7 @@
allocate(haloCursor % array(haloCursor % dimSizes(1)))
+ ! Add all vertices into block, and figure out which are new vertices meaning they belong to the new halo layer
j = 1
do i = 1, nVerticesLocal
if(.not. mpas_hash_search(vertexList(iBlock), localVertexList(i))) then
@@ -832,6 +933,7 @@
deallocate(localVertexList)
+ ! Advance Cursors
indexToVertexCursor => indexToVertexCursor % next
nEdgesCursor => nEdgesCursor % next
nCellsSolveCursor => nCellsSolveCursor % next
@@ -841,13 +943,16 @@
offSetCursor => offSetCursor % next
end do
+ ! Build current layers exchange list
call mpas_dmpar_get_exch_list(iHalo-1, indexToVertexID, haloIndices, offSetField, vertexLimitField)
+ ! Loop over blocks
indexToVertexCursor => indexToVertexID
cellsOnVertexCursor => cellsOnVertex
nVerticesSolveCursor => nVerticesSolve
haloCursor => haloIndices
do while(associated(indexToVertexCursor))
+ ! Copy in new halo indices
array1dHolder => indexToVertexCursor % array
indexToVertexCursor % dimSizes(1) = nVerticesSolveCursor % array(iHalo)
allocate(indexToVertexCursor % array(indexToVertexCursor % dimSizes(1)))
@@ -855,20 +960,24 @@
indexToVertexCursor % array(nVerticesSolveCursor % array(iHalo-1)+1:nVerticesSolveCursor % array(iHalo)) = haloCursor % array(:)
deallocate(array1dHolder)
+ ! Allocate space in cellsOnVertex
array2dHolder => cellsOnVertexCursor % array
cellsOnVertexCursor % dimSizes(2) = nVerticesSolveCursor % array(iHalo)
allocate(cellsOnVertexCursor % array(cellsOnVertexCursor % dimSizes(1), cellsOnVertexCursor % dimSizes(2)))
cellsOnVertexCursor % array(:,1:nVerticesSolveCursor % array(iHalo-1)) = array2dHolder(:,:)
deallocate(array2dHolder)
+ ! Deallocate haloCursor array
deallocate(haloCursor % array)
+ ! Advance cursors
indexToVertexCursor => indexToVertexCursor % next
cellsOnVertexCursor => cellsOnVertexCursor % next
nVerticesSolveCursor => nVerticesSolveCursor % next
haloCursor => haloCursor % next
end do
+ ! Performe allToAll communication
call mpas_dmpar_alltoall_field(cellsOnVertex, cellsOnVertex, sendingHaloLayers)
end do
@@ -885,47 +994,272 @@
end subroutine mpas_block_creator_build_vertex_halos!}}}
-
!***********************************************************************
!
-! routine mpas_get_halo_cells_and_exchange_lists
+! routine mpas_block_creator_finalize_block_init
!
-!> \brief Determines cell indices for each halo layer, and builds exchange lists
+!> \brief Finalize block creation
!> \author Doug Jacobsen
-!> \date 04/30/12
+!> \date 05/31/12
!> \version SVN:$Id$
!> \details
-!> This routine builds the field indexToCellID_nHalos, which is an array of nHalos linked lists.
-!> Each index in indexToCellID_nHalos represnts a linked list of cells in a given halo.
-!> It creates the exchange lists for cells, and places them in the block structure.
-!> In order to call this routine, there are some assumptions made.
-!> The first assumption is that the 1 index of each array is setup correctly,
-!> ie block pointers are valid, dimSizes are valid, next pointers are valid, ets
-!> The second assumption is that the arrays in each field are allocated and full with their appropriate information.
-!> These assumptions lead to the conclusion that the 0 halo has to be properly setup prior to calling this routine.
+!> This routine finalizes the block initialization processor. It calls
+!> mpas_block_allocate to allocate space for all fields in a block. Then the 0
+!> halo indices for each element and the exchange lists are copied into the
+!> appropriate block. A halo update is required after this routien is called
+!> to make sure all data in a block is valid.
!
!-----------------------------------------------------------------------
+ subroutine mpas_block_creator_finalize_block_init(blocklist, maxEdges, maxEdges2, vertexDegree, nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, indexToCellID, indexToEdgeID, indexToVertexID)!{{{
+ type (block_type), pointer :: blocklist !< Input/Output: Linked List of blocks
+ integer :: maxEdges !< Input: maxEdges dimension
+ integer :: maxEdges2 !< Input: maxEdges2 dimension
+ integer :: vertexDegree !< Input: vertexDegree dimension
+ integer :: nVertLevels !< Input: nVertLEvels dimension
+ type (field1dInteger), pointer :: nCellsSolve !< Input: nCellsSolve field information
+ type (field1dInteger), pointer :: nEdgesSolve !< Input: nEdgesSolve field information
+ type (field1dInteger), pointer :: nVerticesSolve !< Input: nVerticesSolve field information
+ type (field1dInteger), pointer :: indexToCellID !< Input: indexToCellID field information
+ type (field1dInteger), pointer :: indexToEdgeID !< Input: indexToEdgeID field information
+ type (field1dInteger), pointer :: indexToVertexID !< Input: indexToVertexID field information
+
+ type (domain_type), pointer :: domain
+
+ type (block_type), pointer :: block_ptr
+ type (field1dInteger), pointer :: nCellsCursor, nEdgesCursor, nVerticesCursor
+ type (field1dInteger), pointer :: indexToCellCursor, indexToEdgeCursor, indexToVertexCursor
+
+ integer :: nCells, nVertices, nEdges, nHalos
+ integer :: nCellsSolve_0Halo, nVerticesSolve_0Halo, nEdgesSolve_0Halo
+ integer :: blockID, localBlockID
+
+ nHalos = config_num_halos
+ domain => blocklist % domain
+
+ ! Loop over blocks
+ block_ptr => blocklist
+ nCellsCursor => nCellsSolve
+ nEdgesCursor => nEdgesSolve
+ nVerticesCursor => nVerticesSolve
+ indexToCellCursor => indexToCellID
+ indexToEdgeCursor => indexToEdgeID
+ indexToVertexCursor => indexToVertexID
+ do while(associated(block_ptr))
+ ! Determine block dimensions
+ nCells = nCellsCursor % array(nHalos+1)
+ nEdges = nEdgesCursor % array(nHalos+2)
+ nVertices = nVerticesCursor % array(nHalos+2)
+
+ nCellsSolve_0Halo = nCellsCursor % array(1)
+ nEdgesSolve_0Halo = nEdgesCursor % array(1)
+ nVerticesSolve_0Halo = nVerticesCursor % array(1)
+
+ ! Determine block IDs
+ blockID = block_ptr % blockID
+ localBlockID = block_ptr % localBlockID
+
+ ! Allocate fields in block
+ call mpas_allocate_block(nHalos, block_ptr, domain, blockID, &
+#include "dim_dummy_args.inc"
+ )
+
+ ! Set block's local id
+ block_ptr % localBlockID = localBlockID
+
+ ! Set block's *Solve dimensions
+ block_ptr % mesh % nCellsSolve = nCellsSolve_0Halo
+ block_ptr % mesh % nEdgesSolve = nEdgesSolve_0Halo
+ block_ptr % mesh % nVerticesSolve = nVerticesSolve_0Halo
+
+ ! Set block's 0 halo indices
+ block_ptr % mesh % indexToCellID % array(1:nCellsSolve_0Halo) = indexToCellCursor % array(1:nCellsSolve_0Halo)
+ block_ptr % mesh % indexToEdgeID % array(1:nEdgesSolve_0Halo) = indexToEdgeCursor % array(1:nEdgesSolve_0Halo)
+ block_ptr % mesh % indexToVertexID % array(1:nVerticesSolve_0Halo) = indexToVertexCursor % array(1:nVerticesSolve_0Halo)
+
+ ! Set block's exchange lists and nullify unneeded exchange lists
+ block_ptr % parinfo % cellsToSend => indexToCellCursor % sendList
+ block_ptr % parinfo % cellsToRecv => indexToCellCursor % recvList
+ block_ptr % parinfo % cellsToCopy => indexToCellCursor % copyList
+ nullify(indexToCellCursor % sendList)
+ nullify(indexToCellCursor % recvList)
+ nullify(indexToCellCursor % copyList)
+
+ block_ptr % parinfo % edgesToSend => indexToEdgeCursor % sendList
+ block_ptr % parinfo % edgesToRecv => indexToEdgeCursor % recvList
+ block_ptr % parinfo % edgesToCopy => indexToEdgeCursor % copyList
+ nullify(indexToEdgeCursor % sendList)
+ nullify(indexToEdgeCursor % recvList)
+ nullify(indexToEdgeCursor % copyList)
+
+ block_ptr % parinfo % verticesToSend => indexToVertexCursor % sendList
+ block_ptr % parinfo % verticesToRecv => indexToVertexCursor % recvList
+ block_ptr % parinfo % verticesToCopy => indexToVertexCursor % copyList
+ nullify(indexToVertexCursor % sendList)
+ nullify(indexToVertexCursor % recvList)
+ nullify(indexToVertexCursor % copyList)
+
+ ! Advance cursors
+ block_ptr => block_ptr % next
+ nCellsCursor => nCellsCursor % next
+ nEdgesCursor => nEdgesCursor % next
+ nVerticesCursor => nVerticesCursor % next
+ indexToCellCursor => indexToCellCursor % next
+ indexToEdgeCursor => indexToEdgeCursor % next
+ indexToVertexCursor => indextoVertexcursor % next
+ end do
+
+ ! Link fields between blocks
+ block_ptr => blocklist
+ do while(associated(block_ptr))
+ call mpas_create_field_links(block_ptr)
+
+ block_ptr => block_ptr % next
+ end do
+ end subroutine mpas_block_creator_finalize_block_init!}}}
+
!***********************************************************************
!
-! routine mpas_get_vertex_ids_and_exchange_lists
+! routine mpas_block_creator_reindex_block_fields
!
-!> \brief Determines vertex indices for each halo layer, and builds exchange lists
+!> \brief Reindex mesh connectivity arrays
!> \author Doug Jacobsen
-!> \date 05/01/12
+!> \date 05/31/12
!> \version SVN:$Id$
!> \details
-!> This routine fills in the arrays for the indexToVertexID_0Halo, and indexToVertexID_nHalos
-!> indexToVertexID_0Halo represents all vertices in the 0 halo, while indexToVertexID_nHalos represnts
-!> the vertex id's for all vertices in all other halos. It is an array of linked lists where each
-!> index represents the linked list of vertex ids at that halo layer.
-!> It creates the exchange lists for vertices, and places them in the block structure.
-!> In order to call this routine, there are some assumptions made.
-!> The first assumption is that the 1 index of each array is setup correctly,
-!> ie block pointers are valid, dimSizes are valid, next pointers are valid, ets
-!> The second assumption is that the arrays in each field are allocated and full with their appropriate information.
-!> These assumptions lead to the conclusion that the 0 halo has to be properly setup prior to calling this routine.
+!> This routine re-indexes the connectivity arrays for the mesh data
+!> structure. Prior to this routine, all indices are given as global index (which
+!> can later be found in the indexTo* arrays). After this routine is called,
+!> indices are provided as local indices now (1:nCells+1 ... etc).
!
!-----------------------------------------------------------------------
+ subroutine mpas_block_creator_reindex_block_fields(blocklist)!{{{
+ type (block_type), pointer :: blocklist !< Input/Output: Linked list of blocks
+
+ type (block_type), pointer :: block_ptr
+
+ integer :: i, j, k
+ integer, dimension(:,:), pointer :: cellIDSorted, edgeIDSorted, vertexIDSorted
+
+ ! Loop over blocks
+ block_ptr => blocklist
+ do while(associated(block_ptr))
+ !
+ ! Rename vertices in cellsOnCell, edgesOnCell, etc. to local indices
+ !
+ allocate(cellIDSorted(2, block_ptr % mesh % nCells))
+ allocate(edgeIDSorted(2, block_ptr % mesh % nEdges))
+ allocate(vertexIDSorted(2, block_ptr % mesh % nVertices))
+
+ do i=1,block_ptr % mesh % nCells
+ cellIDSorted(1,i) = block_ptr % mesh % indexToCellID % array(i)
+ cellIDSorted(2,i) = i
+ end do
+ call quicksort(block_ptr % mesh % nCells, cellIDSorted)
+
+ do i=1,block_ptr % mesh % nEdges
+ edgeIDSorted(1,i) = block_ptr % mesh % indexToEdgeID % array(i)
+ edgeIDSorted(2,i) = i
+ end do
+ call quicksort(block_ptr % mesh % nEdges, edgeIDSorted)
+
+ do i=1,block_ptr % mesh % nVertices
+ vertexIDSorted(1,i) = block_ptr % mesh % indexToVertexID % array(i)
+ vertexIDSorted(2,i) = i
+ end do
+ call quicksort(block_ptr % mesh % nVertices, vertexIDSorted)
+
+
+ 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, &
+ block_ptr % mesh % cellsOnCell % array(j,i))
+ if (k <= block_ptr % mesh % nCells) then
+ block_ptr % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
+ else
+ block_ptr % mesh % cellsOnCell % array(j,i) = block_ptr % mesh % nCells + 1
+ end if
+
+ k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
+ block_ptr % mesh % edgesOnCell % array(j,i))
+ if (k <= block_ptr % mesh % nEdges) then
+ block_ptr % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
+ else
+ block_ptr % mesh % edgesOnCell % array(j,i) = block_ptr % mesh % nEdges + 1
+ end if
+
+ k = mpas_binary_search(vertexIDSorted, 2, 1, block_ptr % mesh % nVertices, &
+ block_ptr % mesh % verticesOnCell % array(j,i))
+ if (k <= block_ptr % mesh % nVertices) then
+ block_ptr % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
+ else
+ block_ptr % mesh % verticesOnCell % array(j,i) = block_ptr % mesh % nVertices + 1
+ end if
+ end do
+ end do
+
+ do i=1,block_ptr % mesh % nEdges
+ do j=1,2
+
+ k = mpas_binary_search(cellIDSorted, 2, 1, block_ptr % mesh % nCells, &
+ block_ptr % mesh % cellsOnEdge % array(j,i))
+ if (k <= block_ptr % mesh % nCells) then
+ block_ptr % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
+ else
+ block_ptr % mesh % cellsOnEdge % array(j,i) = block_ptr % mesh % nCells + 1
+ end if
+
+ k = mpas_binary_search(vertexIDSorted, 2, 1, block_ptr % mesh % nVertices, &
+ block_ptr % mesh % verticesOnEdge % array(j,i))
+ if (k <= block_ptr % mesh % nVertices) then
+ block_ptr % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
+ else
+ block_ptr % mesh % verticesOnEdge % array(j,i) = block_ptr % mesh % nVertices + 1
+ end if
+
+ end do
+
+ do j=1,block_ptr % mesh % nEdgesOnEdge % array(i)
+
+ k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
+ block_ptr % mesh % edgesOnEdge % array(j,i))
+ if (k <= block_ptr % mesh % nEdges) then
+ block_ptr % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
+ else
+ block_ptr % mesh % edgesOnEdge % array(j,i) = block_ptr % mesh % nEdges + 1
+ end if
+ end do
+ end do
+
+ do i=1,block_ptr % mesh % nVertices
+ do j=1,block_ptr % mesh % vertexDegree
+
+ k = mpas_binary_search(cellIDSorted, 2, 1, block_ptr % mesh % nCells, &
+ block_ptr % mesh % cellsOnVertex % array(j,i))
+ if (k <= block_ptr % mesh % nCells) then
+ block_ptr % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
+ else
+ block_ptr % mesh % cellsOnVertex % array(j,i) = block_ptr % mesh % nCells + 1
+ end if
+
+ k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
+ block_ptr % mesh % edgesOnVertex % array(j,i))
+ if (k <= block_ptr % mesh % nEdges) then
+ block_ptr % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
+ else
+ block_ptr % mesh % edgesOnVertex % array(j,i) = block_ptr % mesh % nEdges + 1
+ end if
+ end do
+ end do
+
+ deallocate(cellIDSorted)
+ deallocate(edgeIDSorted)
+ deallocate(vertexIDSorted)
+
+ block_ptr => block_ptr % next
+ end do
+
+ end subroutine mpas_block_creator_reindex_block_fields!}}}
+
end module mpas_block_creator
Modified: branches/omp_blocks/multiple_blocks/src/framework/mpas_dmpar.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/framework/mpas_dmpar.F 2012-05-30 21:41:48 UTC (rev 1950)
+++ branches/omp_blocks/multiple_blocks/src/framework/mpas_dmpar.F 2012-05-31 19:28:35 UTC (rev 1951)
@@ -4222,6 +4222,8 @@
integer :: nHalos
integer :: i
+ nHalos = size(exchList % halos)
+
do i = 1, nHalos
call mpas_dmpar_destroy_exchange_list(exchList % halos(i) % exchList)
end do
Modified: branches/omp_blocks/multiple_blocks/src/framework/mpas_io_input.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/framework/mpas_io_input.F 2012-05-30 21:41:48 UTC (rev 1950)
+++ branches/omp_blocks/multiple_blocks/src/framework/mpas_io_input.F 2012-05-31 19:28:35 UTC (rev 1951)
@@ -72,13 +72,13 @@
type (field2dInteger), pointer :: cellsOnEdgeField
type (field2dInteger), pointer :: cellsOnVertexField
-#ifdef HAVE_ZOLTAN
-#ifdef _MPI
+!#ifdef HAVE_ZOLTAN
+!#ifdef _MPI
type (field1dReal), pointer :: xCellField, yCellField, zCellField
type (field1dReal), pointer :: xEdgeField, yEdgeField, zEdgeField
type (field1dReal), pointer :: xVertexField, yVertexField, zVertexField
-#endif
-#endif
+!#endif
+!#endif
type (field1DChar) :: xtime
@@ -122,13 +122,13 @@
type (field1DInteger), pointer :: indexToEdgeID_tList
-#ifdef HAVE_ZOLTAN
-#ifdef _MPI
+!#ifdef HAVE_ZOLTAN
+!#ifdef _MPI
type (field1DReal), pointer :: xCell, yCell, zCell
type (field1DReal), pointer :: xEdge, yEdge, zEdge
type (field1DReal), pointer :: xVertex, yVertex, zVertex
-#endif
-#endif
+!#endif
+!#endif
integer, dimension(:,:), pointer :: cellIDSorted, vertexIDSorted, edgeIDSorted
integer, dimension(:), pointer :: local_cell_list, local_edge_list, local_vertex_list
@@ -170,35 +170,32 @@
nHalos = config_num_halos
-
if (config_do_restart) then
+ ! this get followed by set is to ensure that the time is in standard format
+ call mpas_set_time(curr_time=startTime, dateTimeString=config_start_time)
+ call mpas_get_time(curr_time=startTime, dateTimeString=timeStamp)
- ! this get followed by set is to ensure that the time is in standard format
- call mpas_set_time(curr_time=startTime, dateTimeString=config_start_time)
- call mpas_get_time(curr_time=startTime, dateTimeString=timeStamp)
+ call mpas_insert_string_suffix(trim(config_restart_name), timeStamp, filename)
- call mpas_insert_string_suffix(trim(config_restart_name), timeStamp, filename)
-
- input_obj % filename = trim(filename)
- input_obj % stream = STREAM_RESTART
+ input_obj % filename = trim(filename)
+ input_obj % stream = STREAM_RESTART
else
- input_obj % filename = trim(config_input_name)
- input_obj % stream = STREAM_INPUT
+ input_obj % filename = trim(config_input_name)
+ input_obj % stream = STREAM_INPUT
end if
inputHandle = MPAS_io_open(trim(input_obj % filename), MPAS_IO_READ, MPAS_IO_PNETCDF, ierr)
if (ierr /= MPAS_IO_NOERR) then
- write(0,*) ' '
- if (input_obj % stream == STREAM_RESTART) then
- write(0,*) 'Error opening restart file ''', trim(input_obj % filename), ''''
- else if (input_obj % stream == STREAM_INPUT) then
- write(0,*) 'Error opening input file ''', trim(input_obj % filename), ''''
- else if (input_obj % stream == STREAM_SFC) then
- write(0,*) 'Error opening sfc file ''', trim(input_obj % filename), ''''
- end if
- write(0,*) ' '
- call mpas_dmpar_abort(domain % dminfo)
+ write(0,*) ' '
+ if (input_obj % stream == STREAM_RESTART) then
+ write(0,*) 'Error opening restart file ''', trim(input_obj % filename), ''''
+ else if (input_obj % stream == STREAM_INPUT) then
+ write(0,*) 'Error opening input file ''', trim(input_obj % filename), ''''
+ else if (input_obj % stream == STREAM_SFC) then
+ write(0,*) 'Error opening sfc file ''', trim(input_obj % filename), ''''
+ end if
+ write(0,*) ' '
+ call mpas_dmpar_abort(domain % dminfo)
end if
-
!
! Read global number of cells/edges/vertices
@@ -226,413 +223,19 @@
readingBlock % domain => domain
readingBlock % blockID = domain % dminfo % my_proc_id
readingBlock % localBlockID = 0
-
-
+
!
! Allocate and read fields that we will need in order to ultimately work out
! which cells/edges/vertices are owned by each block, and which are ghost
!
- ! Global cell indices
- allocate(indexToCellIDField)
- allocate(indexToCellIDField % ioinfo)
- indexToCellIDField % ioinfo % fieldName = 'indexToCellID'
- indexToCellIDField % ioinfo % start(1) = readCellStart
- indexToCellIDField % ioinfo % count(1) = nReadCells
- allocate(indexToCellIDField % array(nReadCells))
- allocate(readIndices(nReadCells))
- do i=1,nReadCells
- readIndices(i) = i + readCellStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'indexToCellID', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'indexToCellID', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'indexToCellID', indexToCellIDField % array, ierr)
- indexToCellIDField % dimSizes(1) = nReadCells
- indexToCellIDField % block => readingBlock
-! call mpas_dmpar_init_multihalo_exchange_list(indexToCellIDField % sendList, nHalos)
-! call mpas_dmpar_init_mulithalo_exchange_list(indexToCellIDField % recvList, nHalos)
-! call mpas_dmpar_init_multihalo_exchange_list(indexToCellIDField % copyList, nHalos)
- allocate(indexToCellIDField % sendList)
- allocate(indexToCellIDField % recvList)
- allocate(indexToCellIDField % copyList)
- allocate(indexToCellIDField % sendList % halos(nHalos))
- allocate(indexToCellIDField % recvList % halos(nHalos))
- allocate(indexToCellIDField % copyList % halos(nHalos))
- do i = 1, nHalos
- nullify(indexToCellIDField % sendList % halos(i) % exchList)
- nullify(indexToCellIDField % recvList % halos(i) % exchList)
- nullify(indexToCellIDField % copyList % halos(i) % exchList)
- end do
- nullify(indexToCellIDField % next)
-
-#ifdef HAVE_ZOLTAN
-#ifdef _MPI
- ! Cell x-coordinates (in 3d Cartesian space)
- allocate(xCellField)
- allocate(xCellField % ioinfo)
- xCellField % ioinfo % fieldName = 'xCell'
- xCellField % ioinfo % start(1) = readCellStart
- xCellField % ioinfo % count(1) = nReadCells
- allocate(xCellField % array(nReadCells))
- call MPAS_io_inq_var(inputHandle, 'xCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'xCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'xCell', xCellField % array, ierr)
- xCellField % dimSizes(1) = nReadCells
- xCellField % block => readingBlock
- xCellField % sendList => indexToCellIDField % sendList
- xCellField % recvList => indexToCellIDField % recvList
- xCellField % copyList => indexToCellIDField % copyList
- nullify(xCellField % next)
+ call mpas_io_setup_cell_block_fields(inputHandle, nreadCells, readCellStart, readingBlock, maxEdges, indexTocellIDField, xCellField, &
+ yCellField, zCellField, nEdgesOnCellField, cellsOnCellField, edgesOnCellField, verticesOnCellField)
- ! Cell y-coordinates (in 3d Cartesian space)
- allocate(yCellField)
- allocate(yCellField % ioinfo)
- yCellField % ioinfo % fieldName = 'yCell'
- yCellField % ioinfo % start(1) = readCellStart
- yCellField % ioinfo % count(1) = nReadCells
- allocate(yCellField % array(nReadCells))
- call MPAS_io_inq_var(inputHandle, 'yCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'yCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'yCell', yCellField % array, ierr)
- yCellField % sendList => indexToCellIDField % sendList
- yCellField % recvList => indexToCellIDField % recvList
- yCellField % copyList => indexToCellIDField % copyList
- yCellField % dimSizes(1) = nReadCells
- yCellField % block => readingBlock
- nullify(yCellField % next)
+ call mpas_io_setup_edge_block_fields(inputHandle, nReadEdges, readEdgeStart, readingBlock, indexToEdgeIDField, xEdgeField, yEdgeField, zEdgeField, cellsOnEdgeField)
- ! Cell z-coordinates (in 3d Cartesian space)
- allocate(zCellField)
- allocate(zCellField % ioinfo)
- zCellField % ioinfo % fieldName = 'zCell'
- zCellField % ioinfo % start(1) = readCellStart
- zCellField % ioinfo % count(1) = nReadCells
- allocate(zCellField % array(nReadCells))
- call MPAS_io_inq_var(inputHandle, 'zCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'zCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'zCell', zCellField % array, ierr)
- zCellField % dimSizes(1) = nReadCells
- zCellField % block => readingBlock
- zCellField % sendList => indexToCellIDField % sendList
- zCellField % recvList => indexToCellIDField % recvList
- zCellField % copyList => indexToCellIDField % copyList
- nullify(zCellField % next)
-#endif
-#endif
- deallocate(readIndices)
-
-
- ! Global edge indices
- allocate(indexToEdgeIDField)
- allocate(indexToEdgeIDField % ioinfo)
- indexToEdgeIDField % ioinfo % fieldName = 'indexToEdgeID'
- indexToEdgeIDField % ioinfo % start(1) = readEdgeStart
- indexToEdgeIDField % ioinfo % count(1) = nReadEdges
- allocate(indexToEdgeIDField % array(nReadEdges))
- allocate(indexToEdgeIDField % array(nReadEdges))
- allocate(readIndices(nReadEdges))
- do i=1,nReadEdges
- readIndices(i) = i + readEdgeStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'indexToEdgeID', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'indexToEdgeID', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'indexToEdgeID', indexToEdgeIDField % array, ierr)
- indexToEdgeIDField % dimSizes(1) = nREadEdges
- indexToEdgeIDField % block => readingBlock
- allocate(indexToEdgeIDField % sendList)
- allocate(indexToEdgeIDField % recvList)
- allocate(indexToEdgeIDField % copyList)
- allocate(indexToEdgeIDField % sendList % halos(nHalos))
- allocate(indexToEdgeIDField % recvList % halos(nHalos))
- allocate(indexToEdgeIDField % copyList % halos(nHalos))
- do i = 1, nHalos
- nullify(indexToEdgeIDField % sendList % halos(i) % exchList)
- nullify(indexToEdgeIDField % recvList % halos(i) % exchList)
- nullify(indexToEdgeIDField % copyList % halos(i) % exchList)
- end do
- nullify(indexToEdgeIDField % next)
-
-#ifdef HAVE_ZOLTAN
-#ifdef _MPI
- ! Edge x-coordinates (in 3d Cartesian space)
- allocate(xEdgeField)
- allocate(xEdgeField % ioinfo)
- xEdgeField % ioinfo % fieldName = 'xEdge'
- xEdgeField % ioinfo % start(1) = readEdgeStart
- xEdgeField % ioinfo % count(1) = nReadEdges
- allocate(xEdgeField % array(nReadEdges))
- call MPAS_io_inq_var(inputHandle, 'xEdge', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'xEdge', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'xEdge', xEdgeField % array, ierr)
- xEdgeField % dimSizes(1) = nReadEdges
- xEdgeField % block => readingBlock
- xEdgeField % sendList => indexToEdgeIDField % sendList
- xEdgeField % recvList => indexToEdgeIDField % recvList
- xEdgeField % copyList => indexToEdgeIDField % copyList
- nullify(xEdgeField % next)
-
- ! Edge y-coordinates (in 3d Cartesian space)
- allocate(yEdgeField)
- allocate(yEdgeField % ioinfo)
- yEdgeField % ioinfo % fieldName = 'yEdge'
- yEdgeField % ioinfo % start(1) = readEdgeStart
- yEdgeField % ioinfo % count(1) = nReadEdges
- allocate(yEdgeField % array(nReadEdges))
- call MPAS_io_inq_var(inputHandle, 'yEdge', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'yEdge', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'yEdge', yEdgeField % array, ierr)
- yEdgeField % dimSizes(1) = nReadEdges
- yEdgeField % block => readingBlock
- yEdgeField % sendList => indexToEdgeIDField % sendList
- yEdgeField % recvList => indexToEdgeIDField % recvList
- yEdgeField % copyList => indexToEdgeIDField % copyList
- nullify(yEdgeField % next)
-
- ! Edge z-coordinates (in 3d Cartesian space)
- allocate(zEdgeField)
- allocate(zEdgeField % ioinfo)
- zEdgeField % ioinfo % fieldName = 'zEdge'
- zEdgeField % ioinfo % start(1) = readEdgeStart
- zEdgeField % ioinfo % count(1) = nReadEdges
- allocate(zEdgeField % array(nReadEdges))
- call MPAS_io_inq_var(inputHandle, 'zEdge', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'zEdge', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'zEdge', zEdgeField % array, ierr)
- zEdgeField % dimSizes(1) = nReadEdges
- zEdgeField % block => readingBlock
- zEdgeField % sendList => indexToEdgeIDField % sendList
- zEdgeField % recvList => indexToEdgeIDField % recvList
- zEdgeField % copyList => indexToEdgeIDField % copyList
- nullify(zEdgeField % next)
-#endif
-#endif
- deallocate(readIndices)
-
-
- ! Global vertex indices
- allocate(indexToVertexIDField)
- allocate(indexToVertexIDField % ioinfo)
- indexToVertexIDField % ioinfo % fieldName = 'indexToVertexID'
- indexToVertexIDField % ioinfo % start(1) = readVertexStart
- indexToVertexIDField % ioinfo % count(1) = nReadVertices
- allocate(indexToVertexIDField % array(nReadVertices))
- allocate(readIndices(nReadVertices))
- do i=1,nReadVertices
- readIndices(i) = i + readVertexStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'indexToVertexID', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'indexToVertexID', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'indexToVertexID', indexToVertexIDField % array, ierr)
- indexToVertexIDField % dimSizes(1) = nReadVertices
- indexToVertexIDField % block => readingBlock
- allocate(indexToVertexIDField % sendList)
- allocate(indexToVertexIDField % recvList)
- allocate(indexToVertexIDField % copyList)
- allocate(indexToVertexIDField % sendList % halos(nHalos))
- allocate(indexToVertexIDField % recvList % halos(nHalos))
- allocate(indexToVertexIDField % copyList % halos(nHalos))
- do i = 1, nHalos
- nullify(indexToVertexIDField % sendList % halos(i) % exchList)
- nullify(indexToVertexIDField % recvList % halos(i) % exchList)
- nullify(indexToVertexIDField % copyList % halos(i) % exchList)
- end do
- nullify(indexToVertexIDField % next)
-
-#ifdef HAVE_ZOLTAN
-#ifdef _MPI
- ! Vertex x-coordinates (in 3d Cartesian space)
- allocate(xVertexField)
- allocate(xVertexField % ioinfo)
- xVertexField % ioinfo % fieldName = 'xVertex'
- xVertexField % ioinfo % start(1) = readVertexStart
- xVertexField % ioinfo % count(1) = nReadVertices
- allocate(xVertexField % array(nReadVertices))
- call MPAS_io_inq_var(inputHandle, 'xVertex', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'xVertex', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'xVertex', xVertexField % array, ierr)
- xVertexField % dimSizes(1) = nReadVertices
- xVertexField % block => readingBlock
- xVertexField % sendList => indexToVertexIDField % sendList
- xVertexField % recvList => indexToVertexIDField % recvList
- xVertexField % copyList => indexToVertexIDField % copyList
- nullify(xVertexField % next)
-
- ! Vertex y-coordinates (in 3d Cartesian space)
- allocate(yVertexField)
- allocate(yVertexField % ioinfo)
- yVertexField % ioinfo % fieldName = 'yVertex'
- yVertexField % ioinfo % start(1) = readVertexStart
- yVertexField % ioinfo % count(1) = nReadVertices
- allocate(yVertexField % array(nReadVertices))
- call MPAS_io_inq_var(inputHandle, 'yVertex', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'yVertex', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'yVertex', yVertexField % array, ierr)
- yVertexField % dimSizes(1) = nReadVertices
- yVertexField % block => readingBlock
- yVertexField % sendList => indexToVertexIDField % sendList
- yVertexField % recvList => indexToVertexIDField % recvList
- yVertexField % copyList => indexToVertexIDField % copyList
- nullify(yVertexField % next)
-
- ! Vertex z-coordinates (in 3d Cartesian space)
- allocate(zVertexField)
- allocate(zVertexField % ioinfo)
- zVertexField % ioinfo % fieldName = 'zVertex'
- zVertexField % ioinfo % start(1) = readVertexStart
- zVertexField % ioinfo % count(1) = nReadVertices
- allocate(zVertexField % array(nReadVertices))
- call MPAS_io_inq_var(inputHandle, 'zVertex', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'zVertex', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'zVertex', zVertexField % array, ierr)
- zVertexField % dimSizes(1) = nReadVertices
- zVertexField % block => readingBlock
- zVertexField % sendList => indexToVertexIDField % sendList
- zVertexField % recvList => indexToVertexIDField % recvList
- zVertexField % copyList => indexToVertexIDField % copyList
- nullify(zVertexField % next)
-#endif
-#endif
- deallocate(readIndices)
-
- ! Number of cell/edges/vertices adjacent to each cell
- allocate(nEdgesOnCellField)
- allocate(nEdgesOnCellField % ioinfo)
- nEdgesOnCellField % ioinfo % fieldName = 'nEdgesOnCell'
- nEdgesOnCellField % ioinfo % start(1) = readCellStart
- nEdgesOnCellField % ioinfo % count(1) = nReadCells
- allocate(nEdgesOnCellField % array(nReadCells))
- allocate(readIndices(nReadCells))
- do i=1,nReadCells
- readIndices(i) = i + readCellStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'nEdgesOnCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'nEdgesOnCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'nEdgesOnCell', nEdgesOnCellField % array, ierr)
- nEdgesOnCellField % dimSizes(1) = nReadCells
- nEdgesOnCellField % block => readingBlock
- nEdgesOnCellField % sendList => indexToCellIDField % sendList
- nEdgesOnCellField % recvList => indexToCellIDField % recvList
- nEdgesOnCellField % copyList => indexToCellIDField % copyList
- nullify(nEdgesOnCellField % next)
-
- ! Global indices of cells adjacent to each cell
- allocate(cellsOnCellField)
- allocate(cellsOnCellField % ioinfo)
- cellsOnCellField % ioinfo % fieldName = 'cellsOnCell'
- cellsOnCellField % ioinfo % start(1) = 1
- cellsOnCellField % ioinfo % start(2) = readCellStart
- cellsOnCellField % ioinfo % count(1) = maxEdges
- cellsOnCellField % ioinfo % count(2) = nReadCells
- allocate(cellsOnCellField % array(maxEdges,nReadCells))
- call MPAS_io_inq_var(inputHandle, 'cellsOnCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'cellsOnCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'cellsOnCell', cellsOnCellField % array, ierr)
- cellsOnCellField % dimSizes(1) = maxEdges
- cellsOnCellField % dimSizes(2) = nReadCells
- cellsOnCellField % block => readingBlock
- cellsOnCellField % sendList => indexToCellIDField % sendList
- cellsOnCellField % recvList => indexToCellIDField % recvList
- cellsOnCellField % copyList => indexToCellIDField % copyList
- nullify(cellsOnCellField % next)
-
- ! Global indices of edges adjacent to each cell
- allocate(edgesOnCellField)
- allocate(edgesOnCellField % ioinfo)
- edgesOnCellField % ioinfo % fieldName = 'edgesOnCell'
- edgesOnCellField % ioinfo % start(1) = 1
- edgesOnCellField % ioinfo % start(2) = readCellStart
- edgesOnCellField % ioinfo % count(1) = maxEdges
- edgesOnCellField % ioinfo % count(2) = nReadCells
- allocate(edgesOnCellField % array(maxEdges,nReadCells))
- call MPAS_io_inq_var(inputHandle, 'edgesOnCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'edgesOnCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'edgesOnCell', edgesOnCellField % array, ierr)
- edgesOnCellField % dimSizes(1) = maxEdges
- edgesOnCellField % dimSizes(2) = nReadCells
- edgesOnCellField % block => readingBlock
- edgesOnCellField % sendList => indexToCellIDField % sendList
- edgesOnCellField % recvList => indexToCellIDField % recvList
- edgesOnCellField % copyList => indexToCellIDField % copyList
- nullify(edgesOnCellField % next)
-
- ! Global indices of vertices adjacent to each cell
- allocate(verticesOnCellField)
- allocate(verticesOnCellField % ioinfo)
- verticesOnCellField % ioinfo % fieldName = 'verticesOnCell'
- verticesOnCellField % ioinfo % start(1) = 1
- verticesOnCellField % ioinfo % start(2) = readCellStart
- verticesOnCellField % ioinfo % count(1) = maxEdges
- verticesOnCellField % ioinfo % count(2) = nReadCells
- allocate(verticesOnCellField % array(maxEdges,nReadCells))
- call MPAS_io_inq_var(inputHandle, 'verticesOnCell', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'verticesOnCell', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'verticesOnCell', verticesOnCellField % array, ierr)
- deallocate(readIndices)
- verticesOnCellField % dimSizes(1) = maxEdges
- verticesOnCellField % dimSizes(2) = nReadCells
- verticesOnCellField % block => readingBlock
- verticesOnCellField % sendList => indexToCellIDField % sendList
- verticesOnCellField % recvList => indexToCellIDField % recvList
- verticesOnCellField % copyList => indexToCellIDField % copyList
- nullify(verticesOnCellField % next)
-
-
-
- ! Global indices of cells adjacent to each edge
- ! used for determining which edges are owned by a block, where
- ! iEdge is owned iff cellsOnEdge(1,iEdge) is an owned cell
- allocate(cellsOnEdgeField)
- allocate(cellsOnEdgeField % ioinfo)
- cellsOnEdgeField % ioinfo % fieldName = 'cellsOnEdge'
- cellsOnEdgeField % ioinfo % start(1) = 1
- cellsOnEdgeField % ioinfo % start(2) = readEdgeStart
- cellsOnEdgeField % ioinfo % count(1) = 2
- cellsOnEdgeField % ioinfo % count(2) = nReadEdges
- allocate(cellsOnEdgeField % array(2,nReadEdges))
- allocate(readIndices(nReadEdges))
- do i=1,nReadEdges
- readIndices(i) = i + readEdgeStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'cellsOnEdge', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'cellsOnEdge', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'cellsOnEdge', cellsOnEdgeField % array, ierr)
- deallocate(readIndices)
- cellsOnEdgeField % dimSizes(1) = 2
- cellsOnEdgeField % dimSizes(2) = nReadEdges
- cellsOnEdgeField % block => readingBlock
- cellsOnEdgeField % sendList => indexToEdgeIDField % sendList
- cellsOnEdgeField % recvList => indexToEdgeIDField % recvList
- cellsOnEdgeField % copyList => indexToEdgeIDField % copyList
- nullify(cellsOnEdgeField % next)
-
-
- ! Global indices of cells adjacent to each vertex
- ! used for determining which vertices are owned by a block, where
- ! iVtx is owned iff cellsOnVertex(1,iVtx) is an owned cell
- allocate(cellsOnVertexField)
- allocate(cellsOnVertexField % ioinfo)
- cellsOnVertexField % ioinfo % fieldName = 'cellsOnVertex'
- cellsOnVertexField % ioinfo % start(1) = 1
- cellsOnVertexField % ioinfo % start(2) = readVertexStart
- cellsOnVertexField % ioinfo % count(1) = vertexDegree
- cellsOnVertexField % ioinfo % count(2) = nReadVertices
- allocate(cellsOnVertexField % array(vertexDegree,nReadVertices))
- allocate(readIndices(nReadVertices))
- do i=1,nReadVertices
- readIndices(i) = i + readVertexStart - 1
- end do
- call MPAS_io_inq_var(inputHandle, 'cellsOnVertex', ierr=ierr)
- call MPAS_io_set_var_indices(inputHandle, 'cellsOnVertex', readIndices, ierr=ierr)
- call mpas_io_get_var(inputHandle, 'cellsOnVertex', cellsOnVertexField % array, ierr)
- cellsOnVertexField % dimSizes(1) = vertexDegree
- cellsOnVertexField % dimSizes(2) = nReadVertices
- cellsOnVertexField % block => readingBlock
- cellsOnVertexField % sendList => indexToVertexIDField % sendList
- cellsOnVertexField % recvList => indexToVertexIDField % recvList
- cellsOnVertexField % copyList => indexToVertexIDField % copyList
- nullify(cellsOnVertexField % next)
- deallocate(readIndices)
-
+ call mpas_io_setup_vertex_block_fields(inputHandle, nReadVertices, readVertexStart, readingBlock, vertexDegree, indexToVertexIDField, &
+ xVertexField, yVertexField, zVertexField, cellsOnVertexField)
!
! Set up a graph derived data type describing the connectivity for the cells
! that were read by this process
@@ -651,7 +254,6 @@
partial_global_graph_info % vertexID(:) = indexToCellIDField % array(:)
partial_global_graph_info % nAdjacent(:) = nEdgesOnCellField % array(:)
partial_global_graph_info % adjacencyList(:,:) = cellsOnCellField % array(:,:)
-
! TODO: Ensure (by renaming or exchanging) that initial cell range on each proc is contiguous
! This situation may occur when reading a restart file with cells/edges/vertices written
@@ -689,6 +291,7 @@
write(6,*) 'Building vertex halos', nHalos
call mpas_block_creator_build_vertex_halos(indexToCellID_0Halo, nEdgesOnCell_0Halo, nCellsSolveField, edgesOnCell_0Halo, indexToEdgeID_0Halo, cellsOnEdge_0Halo, nEdgesSolveField)
+ ! DWJ DEBUGGING
! int1d_ptr => nCellsSolveField
! do while(associated(int1d_ptr))
! write(6,*) 'nCellsSolve on block', int1d_ptr % block % localBlockID
@@ -719,65 +322,17 @@
! int1d_ptr => int1d_ptr % next
! end do
- ! Allocate blocks, and copy indexTo arrays into blocks
- write(6,*) 'Allocate blocks, and copy indexTo arrays into blocks'
+ ! Allocate blocks, and copy indexTo arrays into blocks
+ write(6,*) 'Allocate blocks, and copy indexTo arrays into blocks'
+ call mpas_block_creator_finalize_block_init(domain % blocklist, maxEdges, maxEdges2, vertexDegree, nVertLevels, nCellsSolveField, nEdgesSolveField, nVerticesSolveField, indexToCellID_0Halo, indexToEdgeID_0Halo, indexToVertexID_0Halo)
+
+ write(6,*) 'initializing input object'
block_ptr => domain % blocklist
- int1d_ptr => nCellsSolveField
- int1d_ptr2 => nVerticesSolveField
- int1d_ptr3 => nEdgesSolveField
- indexToCellCursor => indexToCellID_0Halo
- indexToVertexCursor => indexToVertexID_0Halo
- indexToEdgeCursor => indexToEdgeID_0Halo
do while(associated(block_ptr))
- nCells = int1d_ptr % array(nHalos+1)
- nVertices = int1d_ptr2 % array(nHalos+2)
- nEdges = int1d_ptr3 % array(nHalos+2)
-
- call mpas_allocate_block(nHalos, block_ptr, domain , block_ptr % blockID, &
-#include "dim_dummy_args.inc"
- )
-
- block_ptr % mesh % indexToCellID % array(1:int1d_ptr % array(1)) = indexToCellCursor % array(1:int1d_ptr % array(1))
- block_ptr % mesh % indexToVertexID % array(1:int1d_ptr2 % array(1)) = indexToVertexCursor % array(1:int1d_ptr2 % array(1))
- block_ptr % mesh % indexToEdgeID % array(1:int1d_ptr3 % array(1)) = indexToEdgeCursor % array(1:int1d_ptr3 % array(1))
-
- block_ptr % parinfo % cellsToSend => indexToCellCursor % sendList
- block_ptr % parinfo % cellsToRecv => indexToCellCursor % recvList
- block_ptr % parinfo % cellsToCopy => indexToCellCursor % copyList
- nullify(indexToCellCursor % sendList)
- nullify(indexToCellCursor % recvList)
- nullify(indexToCellCursor % copyList)
-
- block_ptr % parinfo % verticesToSend => indexToVertexCursor % sendList
- block_ptr % parinfo % verticesToRecv => indexToVertexCursor % recvList
- block_ptr % parinfo % verticesToCopy => indexToVertexCursor % copyList
- nullify(indexToVertexCursor % sendList)
- nullify(indexToVertexCursor % recvList)
- nullify(indexToVertexCursor % copyList)
-
- block_ptr % parinfo % edgesToSend => indexToEdgeCursor % sendList
- block_ptr % parinfo % edgesToRecv => indexToEdgeCursor % recvList
- block_ptr % parinfo % edgesToCopy => indexToEdgeCursor % copyList
- nullify(indexToEdgeCursor % sendList)
- nullify(indexToEdgeCursor % recvList)
- nullify(indexToEdgeCursor % copyList)
-
- block_ptr % mesh % nCellsSolve = int1d_ptr % array(1)
- block_ptr % mesh % nVerticesSolve = int1d_ptr2 % array(1)
- block_ptr % mesh % nEdgesSolve = int1d_ptr3 % array(1)
-
+ call mpas_io_input_init(input_obj, block_ptr, domain % dminfo)
block_ptr => block_ptr % next
- int1d_ptr => int1d_ptr % next
- int1d_ptr2 => int1d_ptr2 % next
- int1d_ptr3 => int1d_ptr3 % next
- indexToCellCursor => indexToCellCursor % next
- indexToVertexCursor => indexToVertexCursor % next
- indexToEdgeCursor => indexToEdgeCursor % next
end do
- write(6,*) 'initializing input object'
- call mpas_io_input_init(input_obj, domain % blocklist, domain % dminfo)
-
write(6,*) 'getting file attributes'
call MPAS_readStreamAtt(input_obj % io_stream, 'sphere_radius', r_sphere_radius, ierr)
if (ierr /= MPAS_STREAM_NOERR) then
@@ -832,6 +387,23 @@
write(0,*) 'Restarting model from time ', timeStamp
end if
+ !DWJ DEBUGGING
+ write(6,*) 'BLOCK LOOP 1'
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ write(6,*) 'nEdges', block_ptr % mesh % nEdges
+ write(6,*) 'nEdgesSolve', block_ptr % mesh % nEdgesSolve
+ write(6,*) 'indexToCellID on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % indexToCellID % array(i)
+ end do
+ write(6,*) 'edgesOnCell on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % edgesOnCell % array(:, i)
+ end do
+ block_ptr => block_ptr % next
+ end do
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Do the actual work of reading all fields in from the input or restart file
! For each field:
@@ -850,16 +422,19 @@
call MPAS_io_close(inputHandle, ierr)
!DWJ DEBUGGING
-! write(6,*) 'nEdges', domain % blockList % mesh % nEdges
-! write(6,*) 'nEdgesSolve', domain % blockList % mesh % nEdgesSolve
-! write(6,*) 'indexToCellID 1'
-! do i = 1, domain % blockList % mesh % nCells
-! write(6,*) i, domain % blocklist % mesh % indexToCellID % array(i)
-! end do
-! write(6,*) 'edgesOnCell 1'
-! do i = 1, domain % blockList % mesh % nCells
-! write(6,*) i, domain % blocklist % mesh % edgesOnCell % array(:, i)
-! end do
+ write(6,*) 'BLOCK LOOP 2'
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ write(6,*) 'indexToCellID on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % indexToCellID % array(i)
+ end do
+ write(6,*) 'edgesOnCell on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % edgesOnCell % array(:, i)
+ end do
+ block_ptr => block_ptr % next
+ end do
!
! Exchange halos for all of the fields that were read from the input file
@@ -868,143 +443,45 @@
call mpas_exch_input_field_halos(domain, input_obj)
!DWJ DEBUGGING
-! write(6,*) 'indexToCellID 2'
-! do i = 1, domain % blockList % mesh % nCells
-! write(6,*) i, domain % blocklist % mesh % indexToCellID % array(i)
-! end do
-! write(6,*) 'edgesOnCell 2'
-! do i = 1, domain % blockList % mesh % nCells
-! write(6,*) i, domain % blocklist % mesh % edgesOnCell % array(:, i)
-! end do
-
- write(6,*) 're-index fields'
+ write(6,*) 'BLOCK LOOP 3'
block_ptr => domain % blocklist
do while(associated(block_ptr))
- !
- ! Rename vertices in cellsOnCell, edgesOnCell, etc. to local indices
- !
- allocate(cellIDSorted(2, block_ptr % mesh % nCells))
- allocate(edgeIDSorted(2, block_ptr % mesh % nEdges))
- allocate(vertexIDSorted(2, block_ptr % mesh % nVertices))
-
- do i=1,block_ptr % mesh % nCells
- cellIDSorted(1,i) = block_ptr % mesh % indexToCellID % array(i)
- cellIDSorted(2,i) = i
+ write(6,*) 'indexToCellID on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % indexToCellID % array(i)
end do
- call quicksort(block_ptr % mesh % nCells, cellIDSorted)
-
- do i=1,block_ptr % mesh % nEdges
- edgeIDSorted(1,i) = block_ptr % mesh % indexToEdgeID % array(i)
- edgeIDSorted(2,i) = i
+ write(6,*) 'edgesOnCell on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % edgesOnCell % array(:, i)
end do
- call quicksort(block_ptr % mesh % nEdges, edgeIDSorted)
-
- do i=1,block_ptr % mesh % nVertices
- vertexIDSorted(1,i) = block_ptr % mesh % indexToVertexID % array(i)
- vertexIDSorted(2,i) = i
- end do
- call quicksort(block_ptr % mesh % nVertices, vertexIDSorted)
-
-
- 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, &
- block_ptr % mesh % cellsOnCell % array(j,i))
- if (k <= block_ptr % mesh % nCells) then
- block_ptr % mesh % cellsOnCell % array(j,i) = cellIDSorted(2,k)
- else
- block_ptr % mesh % cellsOnCell % array(j,i) = block_ptr % mesh % nCells + 1
- end if
-
-! write(6,*) 'searching for edge', block_ptr % mesh % edgesOnCell % array(j,i)
- k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
- block_ptr % mesh % edgesOnCell % array(j,i))
- if (k <= block_ptr % mesh % nEdges) then
- block_ptr % mesh % edgesOnCell % array(j,i) = edgeIDSorted(2,k)
- else
-! write(6,*) 'Setting',block_ptr % mesh % edgesOnCell % array(j,i),'to',block_ptr % mesh % nEdges + 1
- block_ptr % mesh % edgesOnCell % array(j,i) = block_ptr % mesh % nEdges + 1
- end if
-
- k = mpas_binary_search(vertexIDSorted, 2, 1, block_ptr % mesh % nVertices, &
- block_ptr % mesh % verticesOnCell % array(j,i))
- if (k <= block_ptr % mesh % nVertices) then
- block_ptr % mesh % verticesOnCell % array(j,i) = vertexIDSorted(2,k)
- else
- block_ptr % mesh % verticesOnCell % array(j,i) = block_ptr % mesh % nVertices + 1
- end if
- end do
- end do
-
- do i=1,block_ptr % mesh % nEdges
- do j=1,2
-
- k = mpas_binary_search(cellIDSorted, 2, 1, block_ptr % mesh % nCells, &
- block_ptr % mesh % cellsOnEdge % array(j,i))
- if (k <= block_ptr % mesh % nCells) then
- block_ptr % mesh % cellsOnEdge % array(j,i) = cellIDSorted(2,k)
- else
- block_ptr % mesh % cellsOnEdge % array(j,i) = block_ptr % mesh % nCells + 1
- end if
-
- k = mpas_binary_search(vertexIDSorted, 2, 1, block_ptr % mesh % nVertices, &
- block_ptr % mesh % verticesOnEdge % array(j,i))
- if (k <= block_ptr % mesh % nVertices) then
- block_ptr % mesh % verticesOnEdge % array(j,i) = vertexIDSorted(2,k)
- else
- block_ptr % mesh % verticesOnEdge % array(j,i) = block_ptr % mesh % nVertices + 1
- end if
-
- end do
-
- do j=1,block_ptr % mesh % nEdgesOnEdge % array(i)
-
- k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
- block_ptr % mesh % edgesOnEdge % array(j,i))
- if (k <= block_ptr % mesh % nEdges) then
- block_ptr % mesh % edgesOnEdge % array(j,i) = edgeIDSorted(2,k)
- else
- block_ptr % mesh % edgesOnEdge % array(j,i) = block_ptr % mesh % nEdges + 1
- end if
-
- end do
- end do
-
- do i=1,block_ptr % mesh % nVertices
- do j=1,vertexDegree
-
- k = mpas_binary_search(cellIDSorted, 2, 1, block_ptr % mesh % nCells, &
- block_ptr % mesh % cellsOnVertex % array(j,i))
- if (k <= block_ptr % mesh % nCells) then
- block_ptr % mesh % cellsOnVertex % array(j,i) = cellIDSorted(2,k)
- else
- block_ptr % mesh % cellsOnVertex % array(j,i) = block_ptr % mesh % nCells + 1
- end if
-
- k = mpas_binary_search(edgeIDSorted, 2, 1, block_ptr % mesh % nEdges, &
- block_ptr % mesh % edgesOnVertex % array(j,i))
- if (k <= block_ptr % mesh % nEdges) then
- block_ptr % mesh % edgesOnVertex % array(j,i) = edgeIDSorted(2,k)
- else
- block_ptr % mesh % edgesOnVertex % array(j,i) = block_ptr % mesh % nEdges + 1
- end if
-
- end do
- end do
-
- deallocate(cellIDSorted)
- deallocate(edgeIDSorted)
- deallocate(vertexIDSorted)
-
block_ptr => block_ptr % next
end do
+ call mpas_block_creator_reindex_block_fields(domain % blocklist)
+
!DWJ DEBUGGING
-! write(6,*) 'edgesOnCell 3'
-! do i = 1, domain % blockList % mesh % nCells
-! write(6,*) i, domain % blocklist % mesh % edgesOnCell % array(:, i)
-! end do
+ write(6,*) 'BLOCK LOOP 4'
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ write(6,*) 'edgesOnCell on block', block_ptr % blockID
+ do i = 1, block_ptr % mesh % nCells
+ write(6,*) i, block_ptr % mesh % edgesOnCell % array(:, i)
+ end do
+ block_ptr => block_ptr % next
+ end do
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToCellIDField % sendList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToCellIDField % recvList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToCellIDField % copyList)
+
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToEdgeIDField % sendList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToEdgeIDField % recvList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToEdgeIDField % copyList)
+
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToVertexIDField % sendList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToVertexIDField % recvList)
+ call mpas_dmpar_destroy_mulithalo_exchange_list(indexToVertexIDField % copyList)
+
call mpas_deallocate_field(indexToCellIDField)
call mpas_deallocate_field(indexToEdgeIDField)
call mpas_deallocate_field(indexToVertexIDField)
@@ -1034,9 +511,6 @@
deallocate(block_start)
deallocate(block_count)
-! write(6,*) 'Stopping'
-! stop
-
#ifdef HAVE_ZOLTAN
#ifdef _MPI
! allocate(xCell(size(local_cell_list)))
@@ -1293,5 +767,434 @@
call MPAS_closeStream(input_obj % io_stream, nferr)
end subroutine mpas_io_input_finalize!}}}
+
+ subroutine mpas_io_setup_cell_block_fields(inputHandle, nReadCells, readCellStart, readingBlock, maxEdges, indexToCellID, xCell, yCell, zCell, nEdgesOnCell, cellsOnCell, edgesOnCell, verticesOnCell)!{{{
+ type (MPAS_IO_Handle_type) :: inputHandle
+ integer, intent(in) :: nReadCells
+ integer, intent(in) :: readCellStart
+ integer, intent(in) :: maxEdges
+ type (block_type), pointer :: readingBlock
+ type (field1dInteger), pointer :: indexToCellID
+ type (field1dReal), pointer :: xCell
+ type (field1dReal), pointer :: yCell
+ type (field1dReal), pointer :: zCell
+ type (field1dInteger), pointer :: nEdgesOnCell
+ type (field2dInteger), pointer :: cellsOnCell
+ type (field2dInteger), pointer :: edgesOnCell
+ type (field2dInteger), pointer :: verticesOnCell
+
+ integer :: i, nHalos
+ integer, dimension(:), pointer :: readIndices
+
+ nHalos = config_num_halos
+
+ !
+ ! Allocate and read fields that we will need in order to ultimately work out
+ ! which cells/edges/vertices are owned by each block, and which are ghost
+ !
+
+ ! Global cell indices
+ allocate(indexToCellID)
+ allocate(indexToCellID % ioinfo)
+ indexToCellID % ioinfo % fieldName = 'indexToCellID'
+ indexToCellID % ioinfo % start(1) = readCellStart
+ indexToCellID % ioinfo % count(1) = nReadCells
+ allocate(indexToCellID % array(nReadCells))
+ allocate(readIndices(nReadCells))
+ do i=1,nReadCells
+ readIndices(i) = i + readCellStart - 1
+ end do
+ call MPAS_io_inq_var(inputHandle, 'indexToCellID', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'indexToCellID', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'indexToCellID', indexToCellID % array, ierr)
+ indexToCellID % dimSizes(1) = nReadCells
+ indexToCellID % block => readingBlock
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToCellID % sendList, nHalos)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToCellID % recvList, nHalos)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToCellID % copyList, nHalos)
+ nullify(indexToCellID % next)
+
+#ifdef HAVE_ZOLTAN
+#ifdef _MPI
+ ! Cell x-coordinates (in 3d Cartesian space)
+ allocate(xCell)
+ allocate(xCell % ioinfo)
+ xCell % ioinfo % fieldName = 'xCell'
+ xCell % ioinfo % start(1) = readCellStart
+ xCell % ioinfo % count(1) = nReadCells
+ allocate(xCell % array(nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'xCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'xCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'xCell', xCell % array, ierr)
+ xCell % dimSizes(1) = nReadCells
+ xCell % block => readingBlock
+ xCell % sendList => indexToCellID % sendList
+ xCell % recvList => indexToCellID % recvList
+ xCell % copyList => indexToCellID % copyList
+ nullify(xCell % next)
+
+ ! Cell y-coordinates (in 3d Cartesian space)
+ allocate(yCell)
+ allocate(yCell % ioinfo)
+ yCell % ioinfo % fieldName = 'yCell'
+ yCell % ioinfo % start(1) = readCellStart
+ yCell % ioinfo % count(1) = nReadCells
+ allocate(yCell % array(nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'yCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'yCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'yCell', yCell % array, ierr)
+ yCell % sendList => indexToCellID % sendList
+ yCell % recvList => indexToCellID % recvList
+ yCell % copyList => indexToCellID % copyList
+ yCell % dimSizes(1) = nReadCells
+ yCell % block => readingBlock
+ nullify(yCell % next)
+
+ ! Cell z-coordinates (in 3d Cartesian space)
+ allocate(zCell)
+ allocate(zCell % ioinfo)
+ zCell % ioinfo % fieldName = 'zCell'
+ zCell % ioinfo % start(1) = readCellStart
+ zCell % ioinfo % count(1) = nReadCells
+ allocate(zCell % array(nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'zCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'zCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'zCell', zCell % array, ierr)
+ zCell % dimSizes(1) = nReadCells
+ zCell % block => readingBlock
+ zCell % sendList => indexToCellID % sendList
+ zCell % recvList => indexToCellID % recvList
+ zCell % copyList => indexToCellID % copyList
+ nullify(zCell % next)
+#endif
+#endif
+
+ ! Number of cell/edges/vertices adjacent to each cell
+ allocate(nEdgesOnCell)
+ allocate(nEdgesOnCell % ioinfo)
+ nEdgesOnCell % ioinfo % fieldName = 'nEdgesOnCell'
+ nEdgesOnCell % ioinfo % start(1) = readCellStart
+ nEdgesOnCell % ioinfo % count(1) = nReadCells
+ allocate(nEdgesOnCell % array(nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'nEdgesOnCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'nEdgesOnCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'nEdgesOnCell', nEdgesOnCell % array, ierr)
+ nEdgesOnCell % dimSizes(1) = nReadCells
+ nEdgesOnCell % block => readingBlock
+ nEdgesOnCell % sendList => indexToCellID % sendList
+ nEdgesOnCell % recvList => indexToCellID % recvList
+ nEdgesOnCell % copyList => indexToCellID % copyList
+ nullify(nEdgesOnCell % next)
+
+ ! Global indices of cells adjacent to each cell
+ allocate(cellsOnCell)
+ allocate(cellsOnCell % ioinfo)
+ cellsOnCell % ioinfo % fieldName = 'cellsOnCell'
+ cellsOnCell % ioinfo % start(1) = 1
+ cellsOnCell % ioinfo % start(2) = readCellStart
+ cellsOnCell % ioinfo % count(1) = maxEdges
+ cellsOnCell % ioinfo % count(2) = nReadCells
+ allocate(cellsOnCell % array(maxEdges,nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'cellsOnCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'cellsOnCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'cellsOnCell', cellsOnCell % array, ierr)
+ cellsOnCell % dimSizes(1) = maxEdges
+ cellsOnCell % dimSizes(2) = nReadCells
+ cellsOnCell % block => readingBlock
+ cellsOnCell % sendList => indexToCellID % sendList
+ cellsOnCell % recvList => indexToCellID % recvList
+ cellsOnCell % copyList => indexToCellID % copyList
+ nullify(cellsOnCell % next)
+
+ ! Global indices of edges adjacent to each cell
+ allocate(edgesOnCell)
+ allocate(edgesOnCell % ioinfo)
+ edgesOnCell % ioinfo % fieldName = 'edgesOnCell'
+ edgesOnCell % ioinfo % start(1) = 1
+ edgesOnCell % ioinfo % start(2) = readCellStart
+ edgesOnCell % ioinfo % count(1) = maxEdges
+ edgesOnCell % ioinfo % count(2) = nReadCells
+ allocate(edgesOnCell % array(maxEdges,nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'edgesOnCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'edgesOnCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'edgesOnCell', edgesOnCell % array, ierr)
+ edgesOnCell % dimSizes(1) = maxEdges
+ edgesOnCell % dimSizes(2) = nReadCells
+ edgesOnCell % block => readingBlock
+ edgesOnCell % sendList => indexToCellID % sendList
+ edgesOnCell % recvList => indexToCellID % recvList
+ edgesOnCell % copyList => indexToCellID % copyList
+ nullify(edgesOnCell % next)
+
+ ! Global indices of vertices adjacent to each cell
+ allocate(verticesOnCell)
+ allocate(verticesOnCell % ioinfo)
+ verticesOnCell % ioinfo % fieldName = 'verticesOnCell'
+ verticesOnCell % ioinfo % start(1) = 1
+ verticesOnCell % ioinfo % start(2) = readCellStart
+ verticesOnCell % ioinfo % count(1) = maxEdges
+ verticesOnCell % ioinfo % count(2) = nReadCells
+ allocate(verticesOnCell % array(maxEdges,nReadCells))
+ call MPAS_io_inq_var(inputHandle, 'verticesOnCell', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'verticesOnCell', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'verticesOnCell', verticesOnCell % array, ierr)
+ verticesOnCell % dimSizes(1) = maxEdges
+ verticesOnCell % dimSizes(2) = nReadCells
+ verticesOnCell % block => readingBlock
+ verticesOnCell % sendList => indexToCellID % sendList
+ verticesOnCell % recvList => indexToCellID % recvList
+ verticesOnCell % copyList => indexToCellID % copyList
+ nullify(verticesOnCell % next)
+
+ deallocate(readIndices)
+
+ end subroutine mpas_io_setup_cell_block_fields!}}}
+
+ subroutine mpas_io_setup_edge_block_fields(inputHandle, nReadEdges, readEdgeStart, readingBlock, indexToEdgeID, xEdge, yEdge, zEdge, cellsOnEdge)!{{{
+ type (MPAS_IO_Handle_type) :: inputHandle
+ integer, intent(in) :: nReadEdges
+ integer, intent(in) :: readEdgeStart
+ type (block_type), pointer :: readingBlock
+ type (field1dInteger), pointer :: indexToEdgeID
+ type (field1dReal), pointer :: xEdge
+ type (field1dReal), pointer :: yEdge
+ type (field1dReal), pointer :: zEdge
+ type (field2dInteger), pointer :: cellsOnEdge
+
+ integer :: i, nHalos
+ integer, dimension(:), pointer :: readIndices
+
+ nHalos = config_num_halos
+
+ !
+ ! Allocate and read fields that we will need in order to ultimately work out
+ ! which cells/edges/vertices are owned by each block, and which are ghost
+ !
+
+ allocate(readIndices(nReadEdges))
+
+ ! Global edge indices
+ allocate(indexToEdgeID)
+ allocate(indexToEdgeID % ioinfo)
+ indexToEdgeID % ioinfo % fieldName = 'indexToEdgeID'
+ indexToEdgeID % ioinfo % start(1) = readEdgeStart
+ indexToEdgeID % ioinfo % count(1) = nReadEdges
+ allocate(indexToEdgeID % array(nReadEdges))
+ allocate(indexToEdgeID % array(nReadEdges))
+ do i=1,nReadEdges
+ readIndices(i) = i + readEdgeStart - 1
+ end do
+ call MPAS_io_inq_var(inputHandle, 'indexToEdgeID', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'indexToEdgeID', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'indexToEdgeID', indexToEdgeID % array, ierr)
+ indexToEdgeID % dimSizes(1) = nREadEdges
+ indexToEdgeID % block => readingBlock
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToEdgeID % sendList, nHalos+1)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToEdgeID % recvList, nHalos+1)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToEdgeID % copyList, nHalos+1)
+ nullify(indexToEdgeID % next)
+
+#ifdef HAVE_ZOLTAN
+#ifdef _MPI
+ ! Edge x-coordinates (in 3d Cartesian space)
+ allocate(xEdge)
+ allocate(xEdge % ioinfo)
+ xEdge % ioinfo % fieldName = 'xEdge'
+ xEdge % ioinfo % start(1) = readEdgeStart
+ xEdge % ioinfo % count(1) = nReadEdges
+ allocate(xEdge % array(nReadEdges))
+ call MPAS_io_inq_var(inputHandle, 'xEdge', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'xEdge', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'xEdge', xEdge % array, ierr)
+ xEdge % dimSizes(1) = nReadEdges
+ xEdge % block => readingBlock
+ xEdge % sendList => indexToEdgeID % sendList
+ xEdge % recvList => indexToEdgeID % recvList
+ xEdge % copyList => indexToEdgeID % copyList
+ nullify(xEdge % next)
+
+ ! Edge y-coordinates (in 3d Cartesian space)
+ allocate(yEdge)
+ allocate(yEdge % ioinfo)
+ yEdge % ioinfo % fieldName = 'yEdge'
+ yEdge % ioinfo % start(1) = readEdgeStart
+ yEdge % ioinfo % count(1) = nReadEdges
+ allocate(yEdge % array(nReadEdges))
+ call MPAS_io_inq_var(inputHandle, 'yEdge', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'yEdge', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'yEdge', yEdge % array, ierr)
+ yEdge % dimSizes(1) = nReadEdges
+ yEdge % block => readingBlock
+ yEdge % sendList => indexToEdgeID % sendList
+ yEdge % recvList => indexToEdgeID % recvList
+ yEdge % copyList => indexToEdgeID % copyList
+ nullify(yEdge % next)
+
+ ! Edge z-coordinates (in 3d Cartesian space)
+ allocate(zEdge)
+ allocate(zEdge % ioinfo)
+ zEdge % ioinfo % fieldName = 'zEdge'
+ zEdge % ioinfo % start(1) = readEdgeStart
+ zEdge % ioinfo % count(1) = nReadEdges
+ allocate(zEdge % array(nReadEdges))
+ call MPAS_io_inq_var(inputHandle, 'zEdge', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'zEdge', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'zEdge', zEdge % array, ierr)
+ zEdge % dimSizes(1) = nReadEdges
+ zEdge % block => readingBlock
+ zEdge % sendList => indexToEdgeID % sendList
+ zEdge % recvList => indexToEdgeID % recvList
+ zEdge % copyList => indexToEdgeID % copyList
+ nullify(zEdge % next)
+#endif
+#endif
+
+
+ ! Global indices of cells adjacent to each edge
+ ! used for determining which edges are owned by a block, where
+ ! iEdge is owned iff cellsOnEdge(1,iEdge) is an owned cell
+ allocate(cellsOnEdge)
+ allocate(cellsOnEdge % ioinfo)
+ cellsOnEdge % ioinfo % fieldName = 'cellsOnEdge'
+ cellsOnEdge % ioinfo % start(1) = 1
+ cellsOnEdge % ioinfo % start(2) = readEdgeStart
+ cellsOnEdge % ioinfo % count(1) = 2
+ cellsOnEdge % ioinfo % count(2) = nReadEdges
+ allocate(cellsOnEdge % array(2,nReadEdges))
+ call MPAS_io_inq_var(inputHandle, 'cellsOnEdge', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'cellsOnEdge', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'cellsOnEdge', cellsOnEdge % array, ierr)
+ cellsOnEdge % dimSizes(1) = 2
+ cellsOnEdge % dimSizes(2) = nReadEdges
+ cellsOnEdge % block => readingBlock
+ cellsOnEdge % sendList => indexToEdgeID % sendList
+ cellsOnEdge % recvList => indexToEdgeID % recvList
+ cellsOnEdge % copyList => indexToEdgeID % copyList
+ nullify(cellsOnEdge % next)
+
+ deallocate(readIndices)
+
+ end subroutine mpas_io_setup_edge_block_fields!}}}
+
+ subroutine mpas_io_setup_vertex_block_fields(inputHandle, nReadVertices, readVertexStart, readingBlock, vertexDegree, indexToVertexID, xVertex, yVertex, zVertex, cellsOnVertex)!{{{
+ type (MPAS_IO_Handle_type) :: inputHandle
+ integer, intent(in) :: nReadVertices
+ integer, intent(in) :: readVertexStart
+ integer, intent(in) :: vertexDegree
+ type (block_type), pointer :: readingBlock
+ type (field1dInteger), pointer :: indexToVertexID
+ type (field1dReal), pointer :: xVertex
+ type (field1dReal), pointer :: yVertex
+ type (field1dReal), pointer :: zVertex
+ type (field2dInteger), pointer :: cellsOnVertex
+
+ integer :: i, nHalos
+ integer, dimension(:), pointer :: readIndices
+
+ nHalos = config_num_halos
+
+ ! Global vertex indices
+ allocate(indexToVertexID)
+ allocate(indexToVertexID % ioinfo)
+ indexToVertexID % ioinfo % fieldName = 'indexToVertexID'
+ indexToVertexID % ioinfo % start(1) = readVertexStart
+ indexToVertexID % ioinfo % count(1) = nReadVertices
+ allocate(indexToVertexID % array(nReadVertices))
+ allocate(readIndices(nReadVertices))
+ do i=1,nReadVertices
+ readIndices(i) = i + readVertexStart - 1
+ end do
+ call MPAS_io_inq_var(inputHandle, 'indexToVertexID', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'indexToVertexID', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'indexToVertexID', indexToVertexID % array, ierr)
+ indexToVertexID % dimSizes(1) = nReadVertices
+ indexToVertexID % block => readingBlock
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexID % sendList, nHalos+1)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexID % recvList, nHalos+1)
+ call mpas_dmpar_init_mulithalo_exchange_list(indexToVertexID % copyList, nHalos+1)
+ nullify(indexToVertexID % next)
+
+#ifdef HAVE_ZOLTAN
+#ifdef _MPI
+ ! Vertex x-coordinates (in 3d Cartesian space)
+ allocate(xVertex)
+ allocate(xVertex % ioinfo)
+ xVertex % ioinfo % fieldName = 'xVertex'
+ xVertex % ioinfo % start(1) = readVertexStart
+ xVertex % ioinfo % count(1) = nReadVertices
+ allocate(xVertex % array(nReadVertices))
+ call MPAS_io_inq_var(inputHandle, 'xVertex', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'xVertex', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'xVertex', xVertex % array, ierr)
+ xVertex % dimSizes(1) = nReadVertices
+ xVertex % block => readingBlock
+ xVertex % sendList => indexToVertexID % sendList
+ xVertex % recvList => indexToVertexID % recvList
+ xVertex % copyList => indexToVertexID % copyList
+ nullify(xVertex % next)
+
+ ! Vertex y-coordinates (in 3d Cartesian space)
+ allocate(yVertex)
+ allocate(yVertex % ioinfo)
+ yVertex % ioinfo % fieldName = 'yVertex'
+ yVertex % ioinfo % start(1) = readVertexStart
+ yVertex % ioinfo % count(1) = nReadVertices
+ allocate(yVertex % array(nReadVertices))
+ call MPAS_io_inq_var(inputHandle, 'yVertex', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'yVertex', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'yVertex', yVertex % array, ierr)
+ yVertex % dimSizes(1) = nReadVertices
+ yVertex % block => readingBlock
+ yVertex % sendList => indexToVertexID % sendList
+ yVertex % recvList => indexToVertexID % recvList
+ yVertex % copyList => indexToVertexID % copyList
+ nullify(yVertex % next)
+
+ ! Vertex z-coordinates (in 3d Cartesian space)
+ allocate(zVertex)
+ allocate(zVertex % ioinfo)
+ zVertex % ioinfo % fieldName = 'zVertex'
+ zVertex % ioinfo % start(1) = readVertexStart
+ zVertex % ioinfo % count(1) = nReadVertices
+ allocate(zVertex % array(nReadVertices))
+ call MPAS_io_inq_var(inputHandle, 'zVertex', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'zVertex', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'zVertex', zVertex % array, ierr)
+ zVertex % dimSizes(1) = nReadVertices
+ zVertex % block => readingBlock
+ zVertex % sendList => indexToVertexID % sendList
+ zVertex % recvList => indexToVertexID % recvList
+ zVertex % copyList => indexToVertexID % copyList
+ nullify(zVertex % next)
+#endif
+#endif
+
+
+ ! Global indices of cells adjacent to each vertex
+ ! used for determining which vertices are owned by a block, where
+ ! iVtx is owned iff cellsOnVertex(1,iVtx) is an owned cell
+ allocate(cellsOnVertex)
+ allocate(cellsOnVertex % ioinfo)
+ cellsOnVertex % ioinfo % fieldName = 'cellsOnVertex'
+ cellsOnVertex % ioinfo % start(1) = 1
+ cellsOnVertex % ioinfo % start(2) = readVertexStart
+ cellsOnVertex % ioinfo % count(1) = vertexDegree
+ cellsOnVertex % ioinfo % count(2) = nReadVertices
+ allocate(cellsOnVertex % array(vertexDegree,nReadVertices))
+ call MPAS_io_inq_var(inputHandle, 'cellsOnVertex', ierr=ierr)
+ call MPAS_io_set_var_indices(inputHandle, 'cellsOnVertex', readIndices, ierr=ierr)
+ call mpas_io_get_var(inputHandle, 'cellsOnVertex', cellsOnVertex % array, ierr)
+ cellsOnVertex % dimSizes(1) = vertexDegree
+ cellsOnVertex % dimSizes(2) = nReadVertices
+ cellsOnVertex % block => readingBlock
+ cellsOnVertex % sendList => indexToVertexID % sendList
+ cellsOnVertex % recvList => indexToVertexID % recvList
+ cellsOnVertex % copyList => indexToVertexID % copyList
+ nullify(cellsOnVertex % next)
+
+ deallocate(readIndices)
+
+ end subroutine mpas_io_setup_vertex_block_fields!}}}
+
end module mpas_io_input
</font>
</pre>