<p><b>duda</b> 2012-03-22 18:24:20 -0600 (Thu, 22 Mar 2012)</p><p>BRANCH COMMIT<br>
<br>
Factor out reusable parts of MPAS_streamAddField to make it easier to implement<br>
this routine (using less code, too) for field types besides 1d integer fields.<br>
<br>
<br>
M src/framework/mpas_io_streams.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/io/src/framework/mpas_io_streams.F
===================================================================
--- branches/omp_blocks/io/src/framework/mpas_io_streams.F 2012-03-22 22:53:57 UTC (rev 1697)
+++ branches/omp_blocks/io/src/framework/mpas_io_streams.F 2012-03-23 00:24:20 UTC (rev 1698)
@@ -117,7 +117,8 @@
integer :: io_err
integer :: i
integer :: idim
- integer :: localDimSize, globalDimSize
+ integer :: totalDimSize, globalDimSize
+ logical :: isDecomposed
integer :: ndims
type (field1dInteger), pointer :: field_ptr
character (len=64), dimension(5) :: dimNames
@@ -142,18 +143,126 @@
write(0,*) '... field has ', ndims, ' dimensions'
- allocate(new_field_list_node)
- nullify(new_field_list_node % next)
+ !
+ ! Determine whether the field is decomposed, the indices that are owned by this task's blocks,
+ ! and the total number of outer-indices owned by this task
+ !
+ idim = ndims
+ totalDimSize = 0
+ field_ptr => field
+ if (trim(field % dimNames(idim)) == 'nCells') then
+write(0,*) '... outer dimension is nCells'
+ allocate(indices(0))
+ do while (associated(field_ptr))
+ call mergeArrays(indices, field_ptr % block % mesh % indexToCellID % array(1:field_ptr % block % mesh % nCellsSolve))
+ totalDimSize = totalDimSize + field_ptr % block % mesh % nCellsSolve
+ field_ptr => field_ptr % next
+ end do
+ call mpas_dmpar_sum_int(field % block % domain % dminfo, totalDimSize, globalDimSize)
+ isDecomposed = .true.
+ else if (trim(field % dimNames(idim)) == 'nEdges') then
+write(0,*) '... outer dimension is nEdges'
+ allocate(indices(0))
+ do while (associated(field_ptr))
+ call mergeArrays(indices, field_ptr % block % mesh % indexToEdgeID % array(1:field_ptr % block % mesh % nEdgesSolve))
+ totalDimSize = totalDimSize + field_ptr % block % mesh % nEdgesSolve
+ field_ptr => field_ptr % next
+ end do
+ call mpas_dmpar_sum_int(field % block % domain % dminfo, totalDimSize, globalDimSize)
+ isDecomposed = .true.
+ else if (trim(field % dimNames(idim)) == 'nVertices') then
+write(0,*) '... outer dimension is nVertices'
+ allocate(indices(0))
+ do while (associated(field_ptr))
+ call mergeArrays(indices, field_ptr % block % mesh % indexToVertexID % array(1:field_ptr % block % mesh % nVerticesSolve))
+ totalDimSize = totalDimSize + field_ptr % block % mesh % nVerticesSolve
+ field_ptr => field_ptr % next
+ end do
+ call mpas_dmpar_sum_int(field % block % domain % dminfo, totalDimSize, globalDimSize)
+ isDecomposed = .true.
+ else
+ isDecomposed = .false.
+ globalDimSize = field % dimSizes(idim)
+ totalDimSize = globalDimSize
+
+ if (field % block % domain % dminfo % my_proc_id == IO_NODE) then
+ ndims = 1
+ allocate(indices(field % dimSizes(ndims)))
+ do i=1,field % dimSizes(ndims)
+ indices(i) = i
+ end do
+ else
+ allocate(indices(0))
+ end if
+ end if
+
+
+ call MPAS_streamAddField_generic(stream, trim(field % fieldName), MPAS_IO_INT, field % dimNames, field % dimSizes, &
+ field % hasTimeDimension, isDecomposed, totalDimSize, globalDimSize, indices, ierr)
+
+ deallocate(indices)
+
+
+ !
+ ! Set field pointer and type in fieldList
+ !
+ new_field_list_node => stream % fieldList
+ do while (associated(new_field_list_node % next))
+ new_field_list_node => new_field_list_node % next
+ end do
new_field_list_node % field_type = FIELD_1D_INT
new_field_list_node % int1dField => field
+write(0,*) '... done adding field'
+
+ end subroutine MPAS_streamAddField_1dInteger
+
+
+ subroutine MPAS_streamAddField_generic(stream, fieldName, fieldType, dimNames, dimSizes, hasTimeDimension, isDecomposed, &
+ totalDimSize, globalDimSize, indices, ierr)
+
+ implicit none
+
+ type (MPAS_Stream_type), intent(inout) :: stream
+ character (len=*), intent(in) :: fieldName
+ integer, intent(in) :: fieldType
+ character (len=64), dimension(:), intent(in) :: dimNames
+ integer, dimension(:), intent(in) :: dimSizes
+ logical, intent(in) :: hasTimeDimension
+ logical, intent(in) :: isDecomposed
+ integer, intent(in) :: totalDimSize
+ integer, intent(in) :: globalDimSize
+ integer, dimension(:), intent(in) :: indices
+ integer, intent(out), optional :: ierr
+
+ integer :: io_err
+ integer :: i
+ integer :: idim
+ integer :: ndims
+ character (len=64), dimension(5) :: dimNamesLocal
+ character (len=64), dimension(:), pointer :: dimNamesInq
+ type (field_list_type), pointer :: field_list_cursor
+ type (field_list_type), pointer :: new_field_list_node
+
+ if (present(ierr)) ierr = MPAS_STREAM_NOERR
+
+ ndims = size(dimNames)
+
+write(0,*) '... field has ', ndims, ' dimensions'
+
+ allocate(new_field_list_node)
+ nullify(new_field_list_node % next)
+
if (stream % ioDirection == MPAS_IO_WRITE) then
write(0,*) '... defining field'
+ !
+ ! Define inner dimensions
+ !
do idim = 1, ndims-1
-write(0,*) '... defining dimension ', trim(field % dimNames(idim)), field % dimSizes(idim)
- write(dimNames(idim),'(a)') field % dimNames(idim)
- call MPAS_io_def_dim(stream % fileHandle, trim(field % dimNames(idim)), field % dimSizes(idim), io_err)
+write(0,*) '... defining dimension ', trim(dimNames(idim)), dimSizes(idim)
+ write(dimNamesLocal(idim),'(a)') dimNames(idim)
+ call MPAS_io_def_dim(stream % fileHandle, trim(dimNames(idim)), dimSizes(idim), io_err)
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR) then
if (present(ierr)) ierr = MPAS_IO_ERR
@@ -161,54 +270,32 @@
end if
end do
- ! Define outer-most dimension; handle known decomposed dimensions
+ !
+ ! Define outer-most dimension
+ !
idim = ndims
- localDimSize = 0
- field_ptr => field
- write(dimNames(idim),'(a)') field % dimNames(idim)
- if (trim(field % dimNames(idim)) == 'nCells') then
-write(0,*) '... outer dimension is nCells'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nCellsSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
- else if (trim(field % dimNames(idim)) == 'nEdges') then
-write(0,*) '... outer dimension is nEdges'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nEdgesSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
- else if (trim(field % dimNames(idim)) == 'nVertices') then
-write(0,*) '... outer dimension is nVertices'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nVerticesSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
+ write(dimNamesLocal(idim),'(a)') dimNames(idim)
+
+ if (isDecomposed) then
+ new_field_list_node % totalDimSize = totalDimSize
else
- globalDimSize = field % dimSizes(idim)
- new_field_list_node % isDecomposed = .false.
new_field_list_node % totalDimSize = globalDimSize
end if
-write(0,*) '... defining dimension ', trim(field % dimNames(idim)), globalDimSize
- call MPAS_io_def_dim(stream % fileHandle, trim(field % dimNames(idim)), globalDimSize, io_err)
+ new_field_list_node % isDecomposed = isDecomposed
+
+write(0,*) '... defining dimension ', trim(dimNames(idim)), globalDimSize
+ call MPAS_io_def_dim(stream % fileHandle, trim(dimNames(idim)), globalDimSize, io_err)
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR) then
if (present(ierr)) ierr = MPAS_IO_ERR
return
end if
+ !
! Define time dimension if necessary
- if (field % hasTimeDimension) then
+ !
+ if (hasTimeDimension) then
write(0,*) '... defining Time dimension '
call MPAS_io_def_dim(stream % fileHandle, 'Time', MPAS_IO_UNLIMITED_DIM, io_err)
call MPAS_io_err_mesg(io_err, .false.)
@@ -217,15 +304,15 @@
return
end if
ndims = ndims + 1
- write(dimNames(ndims),'(a)') 'Time'
+ write(dimNamesLocal(ndims),'(a)') 'Time'
end if
+ !
! Define variable to low-level interface
+ !
write(0,*) '... defining var to low-level interface with ndims ', ndims
-do i=1,ndims
- write(0,*) ' ... '//trim(dimNames(i))
-end do
- call MPAS_io_def_var(stream % fileHandle, trim(field % fieldName), MPAS_IO_INT, dimNames(1:ndims), io_err)
+
+ call MPAS_io_def_var(stream % fileHandle, trim(fieldName), fieldType, dimNamesLocal(1:ndims), io_err)
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR) then
if (present(ierr)) ierr = MPAS_IO_ERR
@@ -235,54 +322,27 @@
else if (stream % ioDirection == MPAS_IO_READ) then
write(0,*) '... inquiring about'
- call MPAS_io_inq_var(stream % fileHandle, trim(field % fieldName), dimnames=dimNamesInq, ierr=io_err)
+ call MPAS_io_inq_var(stream % fileHandle, trim(fieldName), dimnames=dimNamesInq, ierr=io_err)
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR) then
if (present(ierr)) ierr = MPAS_IO_ERR
return
end if
- do i=1,ndims
-write(0,*) 'Comparing '//trim(field % dimNames(i))//' '//trim(dimNamesInq(i))
- end do
+! Here, we should probably do a check to make sure the file agrees with what MPAS expects for the field
+! do i=1,ndims
+!write(0,*) 'Comparing '//trim(dimNames(i))//' '//trim(dimNamesInq(i))
+! end do
- ! Define outer-most dimension; handle known decomposed dimensions
- idim = ndims
- localDimSize = 0
- field_ptr => field
- if (trim(field % dimNames(idim)) == 'nCells') then
-write(0,*) '... outer dimension is nCells'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nCellsSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
- else if (trim(field % dimNames(idim)) == 'nEdges') then
-write(0,*) '... outer dimension is nEdges'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nEdgesSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
- else if (trim(field % dimNames(idim)) == 'nVertices') then
-write(0,*) '... outer dimension is nVertices'
- do while (associated(field_ptr))
- localDimSize = localDimSize + field_ptr % block % mesh % nVerticesSolve
- field_ptr => field_ptr % next
- end do
- call mpas_dmpar_sum_int(field % block % domain % dminfo, localDimSize, globalDimSize)
- new_field_list_node % isDecomposed = .true.
- new_field_list_node % totalDimSize = localDimSize
+ ! Set outer dimension sizes depending on whether the field is decomposed
+ if (isDecomposed) then
+ new_field_list_node % totalDimSize = totalDimSize
else
- globalDimSize = field % dimSizes(idim)
- new_field_list_node % isDecomposed = .false.
new_field_list_node % totalDimSize = globalDimSize
end if
+ new_field_list_node % isDecomposed = isDecomposed
+
deallocate(dimNamesInq)
end if
@@ -291,49 +351,7 @@
!
! Set variable indices
!
- if (trim(field % dimNames(1)) == 'nCells') then
- allocate(indices(0))
- field_ptr => field
- do while (associated(field_ptr))
- call mergeArrays(indices, field_ptr % block % mesh % indexToCellID % array(1:field_ptr % block % mesh % nCellsSolve))
- field_ptr => field_ptr % next
- end do
- call MPAS_io_set_var_indices(stream % fileHandle, trim(field % fieldName), indices, io_err)
-
- else if (trim(field % dimNames(1)) == 'nEdges') then
- allocate(indices(0))
- field_ptr => field
- do while (associated(field_ptr))
- call mergeArrays(indices, field_ptr % block % mesh % indexToEdgeID % array(1:field_ptr % block % mesh % nEdgesSolve))
- field_ptr => field_ptr % next
- end do
- call MPAS_io_set_var_indices(stream % fileHandle, trim(field % fieldName), indices, io_err)
-
- else if (trim(field % dimNames(1)) == 'nVertices') then
- allocate(indices(0))
- field_ptr => field
- do while (associated(field_ptr))
- call mergeArrays(indices, field_ptr % block % mesh % indexToVertexID % array(1:field_ptr % block % mesh % nVerticesSolve))
- field_ptr => field_ptr % next
- end do
- call MPAS_io_set_var_indices(stream % fileHandle, trim(field % fieldName), indices, io_err)
-
- else
-
- if (field % block % domain % dminfo % my_proc_id == IO_NODE) then
- ndims = 1
- allocate(indices(field % dimSizes(ndims)))
- do i=1,field % dimSizes(ndims)
- indices(i) = i
- end do
- else
- allocate(indices(0))
- end if
- call MPAS_io_set_var_indices(stream % fileHandle, trim(field % fieldName), indices, io_err)
- end if
-
- deallocate(indices)
-
+ call MPAS_io_set_var_indices(stream % fileHandle, trim(fieldName), indices, io_err)
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR) then
if (present(ierr)) ierr = MPAS_IO_ERR
@@ -356,11 +374,9 @@
field_list_cursor % next => new_field_list_node
end if
-write(0,*) '... done adding field'
+ end subroutine MPAS_streamAddField_generic
- end subroutine MPAS_streamAddField_1dInteger
-
subroutine MPAS_readStream(stream, frame, ierr)
implicit none
</font>
</pre>