<p><b>duda</b> 2012-03-22 09:57:46 -0600 (Thu, 22 Mar 2012)</p><p>BRANCH COMMIT<br>
<br>
Add functionality for adding 1d integer fields to a stream, and for writing a stream.<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-21 16:26:35 UTC (rev 1693)
+++ branches/omp_blocks/io/src/framework/mpas_io_streams.F 2012-03-22 15:57:46 UTC (rev 1694)
@@ -4,10 +4,31 @@
use mpas_grid_types
use mpas_io
+ type field_list_type
+ integer :: field_type
+ logical :: isDecomposed
+ integer :: totalDimSize ! Total size of outer dimension across all blocks for decomposed fields
+ type (field0dInteger), pointer :: int0dField => null()
+ type (field1dInteger), pointer :: int1dField => null()
+ type (field2dInteger), pointer :: int2dField => null()
+ type (field3dInteger), pointer :: int3dField => null()
+ type (field0dReal), pointer :: real0dField => null()
+ type (field1dReal), pointer :: real1dField => null()
+ type (field2dReal), pointer :: real2dField => null()
+ type (field3dReal), pointer :: real3dField => null()
+ type (field0dChar), pointer :: char0dField => null()
+ type (field1dChar), pointer :: char1dField => null()
+ type (field_list_type), pointer :: next => null()
+ end type field_list_type
+
type MPAS_Stream_type
logical :: isInitialized = .false.
+ integer :: ioFormat
+ integer :: ioDirection
+ integer :: framesPerFile
type (MPAS_IO_Handle_type) :: fileHandle
type (att_list_type), pointer :: attList => null()
+ type (field_list_type), pointer :: fieldList => null()
end type MPAS_Stream_type
@@ -31,12 +52,26 @@
module procedure MPAS_writeStreamAtt_text
end interface MPAS_writeStreamAtt
+
! Error codes
integer, parameter :: MPAS_STREAM_NOERR = 0, &
MPAS_STREAM_NOT_INITIALIZED = -1, &
MPAS_IO_ERR = -2
+ integer, parameter :: FIELD_0D_INT = 1, &
+ FIELD_1D_INT = 2, &
+ FIELD_2D_INT = 3, &
+ FIELD_3D_INT = 4, &
+ FIELD_0D_REAL = 5, &
+ FIELD_1D_REAL = 6, &
+ FIELD_2D_REAL = 7, &
+ FIELD_3D_REAL = 8, &
+ FIELD_0D_CHAR = 9, &
+ FIELD_1D_CHAR = 10
+ private mergeArrays
+
+
contains
@@ -57,8 +92,15 @@
stream % fileHandle = MPAS_io_open(fileName, ioDirection, ioFormat, io_err)
call MPAS_io_err_mesg(io_err, .false.)
- if (io_err /= MPAS_IO_NOERR .and. present(ierr)) ierr = MPAS_IO_ERR
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+ stream % ioDirection = ioDirection
+ stream % ioFormat = ioFormat
+ stream % framesPerFile = framesPerFile
+
stream % isInitialized = .true.
end subroutine MPAS_createStream
@@ -69,9 +111,20 @@
implicit none
type (MPAS_Stream_type), intent(inout) :: stream
- type (field1DInteger), intent(in) :: field
+ type (field1DInteger), intent(in), target :: field
integer, intent(out), optional :: ierr
+ integer :: io_err
+ integer :: i
+ integer :: idim
+ integer :: localDimSize, globalDimSize
+ integer :: ndims
+ type (field1dInteger), pointer :: field_ptr
+ character (len=64), dimension(5) :: dimNames
+ integer, dimension(:), pointer :: indices
+ type (field_list_type), pointer :: field_list_cursor
+ type (field_list_type), pointer :: new_field_list_node
+
if (present(ierr)) ierr = MPAS_STREAM_NOERR
!
@@ -82,6 +135,178 @@
return
end if
+write(0,*) '... Adding field '//trim(field % fieldName)//' to stream'
+
+ ndims = size(field % dimSizes)
+
+write(0,*) '... field has ', ndims, ' dimensions'
+
+ allocate(new_field_list_node)
+ nullify(new_field_list_node % next)
+ new_field_list_node % field_type = FIELD_1D_INT
+ new_field_list_node % int1dField => field
+
+ if (stream % ioDirection == MPAS_IO_WRITE) then
+write(0,*) '... defining field'
+
+ 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)
+ 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
+ end do
+
+ ! Define outer-most dimension; handle known decomposed dimensions
+ 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
+ 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)
+ 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
+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.)
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+ ndims = ndims + 1
+ write(dimNames(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_err_mesg(io_err, .false.)
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+
+ else if (stream % ioDirection == MPAS_IO_READ) then
+write(0,*) '... inquiring about'
+
+ end if
+
+
+ !
+ ! 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_err_mesg(io_err, .false.)
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+
+
+ !
+ ! Add field pointer to the list
+ !
+ if (.not. associated(stream % fieldList)) then
+write(0,*) 'Adding field to the head of the list'
+ stream % fieldList => new_field_list_node
+ else
+write(0,*) 'Adding field to the tail of the list'
+ field_list_cursor => stream % fieldList
+ do while (associated(field_list_cursor % next))
+ field_list_cursor => field_list_cursor % next
+ end do
+ field_list_cursor % next => new_field_list_node
+ end if
+
+write(0,*) '... done adding field'
+
end subroutine MPAS_streamAddField_1dInteger
@@ -114,6 +339,20 @@
integer, intent(in) :: frame
integer, intent(out), optional :: ierr
+ integer :: io_err
+ integer :: i
+ integer :: ownedSize
+ type (field1dInteger), pointer :: field_1dint_ptr
+ type (field_list_type), pointer :: field_cursor
+ integer :: int0d_temp
+ integer, dimension(:), pointer :: int1d_temp
+ integer, dimension(:,:), pointer :: int2d_temp
+ integer, dimension(:,:,:), pointer :: int3d_temp
+ real (kind=RKIND) :: real0d_temp
+ real (kind=RKIND), dimension(:), pointer :: real1d_temp
+ real (kind=RKIND), dimension(:,:), pointer :: real2d_temp
+ real (kind=RKIND), dimension(:,:,:), pointer :: real3d_temp
+
if (present(ierr)) ierr = MPAS_STREAM_NOERR
!
@@ -124,6 +363,69 @@
return
end if
+ call MPAS_io_set_frame(stream % fileHandle, frame, 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
+
+ field_cursor => stream % fieldList
+ do while (associated(field_cursor))
+ if (field_cursor % field_type == FIELD_0D_INT) then
+
+ else if (field_cursor % field_type == FIELD_1D_INT) then
+write(0,*) 'Writing out field '//trim(field_cursor % int1dField % fieldName)
+write(0,*) ' > is the field decomposed? ', field_cursor % isDecomposed
+write(0,*) ' > outer dimension size ', field_cursor % totalDimSize
+
+ allocate(int1d_temp(field_cursor % totalDimSize))
+
+ if (field_cursor % isDecomposed) then
+write(0,*) 'Gathering field from across all blocks'
+ ! Gather field from across multiple blocks
+ field_1dint_ptr => field_cursor % int1dField
+ i = 1
+ do while (associated(field_1dint_ptr))
+! Could this be something we do in streamAddField routines, since we already have such tests?
+ if (trim(field_1dint_ptr % dimNames(1)) == 'nCells') then
+ ownedSize = field_1dint_ptr % block % mesh % nCellsSolve
+ else if (trim(field_1dint_ptr % dimNames(1)) == 'nEdges') then
+ ownedSize = field_1dint_ptr % block % mesh % nEdgesSolve
+ else if (trim(field_1dint_ptr % dimNames(1)) == 'nVertices') then
+ ownedSize = field_1dint_ptr % block % mesh % nVerticesSolve
+ else
+ ownedSize = field_1dint_ptr % dimSizes(1)
+ end if
+write(0,*) ' > copying block ', i, i+ownedSize-1
+ int1d_temp(i:i+ownedSize-1) = field_1dint_ptr % array(1:ownedSize)
+ i = i + ownedSize
+ field_1dint_ptr => field_1dint_ptr % next
+ end do
+ else
+write(0,*) 'Copying field from first block'
+ int1d_temp(:) = field_cursor % int1dField % array(:)
+ end if
+
+write(0,*) 'MGD calling MPAS_io_put_var now...'
+ call MPAS_io_put_var(stream % fileHandle, field_cursor % int1dField % fieldName, int1d_temp, io_err)
+ call MPAS_io_err_mesg(io_err, .false.)
+ if (io_err /= MPAS_IO_NOERR .and. present(ierr)) ierr = MPAS_IO_ERR
+
+ deallocate(int1d_temp)
+
+ else if (field_cursor % field_type == FIELD_2D_INT) then
+ else if (field_cursor % field_type == FIELD_3D_INT) then
+ else if (field_cursor % field_type == FIELD_0D_REAL) then
+ else if (field_cursor % field_type == FIELD_1D_REAL) then
+ else if (field_cursor % field_type == FIELD_2D_REAL) then
+ else if (field_cursor % field_type == FIELD_3D_REAL) then
+ else if (field_cursor % field_type == FIELD_0D_CHAR) then
+ else if (field_cursor % field_type == FIELD_1D_CHAR) then
+ end if
+ field_cursor => field_cursor % next
+ end do
+
end subroutine MPAS_writeStream
@@ -415,6 +717,7 @@
integer, intent(out), optional :: ierr
integer :: io_err
+ type (field_list_type), pointer :: field_cursor
if (present(ierr)) ierr = MPAS_STREAM_NOERR
@@ -430,8 +733,43 @@
call MPAS_io_err_mesg(io_err, .false.)
if (io_err /= MPAS_IO_NOERR .and. present(ierr)) ierr = MPAS_IO_ERR
+write(0,*) 'Deallocating global attribute list'
+ call mpas_deallocate_attlist(stream % attList)
+
+write(0,*) 'Deallocating field list'
+ field_cursor => stream % fieldList
+ do while (associated(field_cursor))
+ stream % fieldList => stream % fieldList % next
+ deallocate(field_cursor)
+ field_cursor => stream % fieldList
+ end do
+
stream % isInitialized = .false.
end subroutine MPAS_closeStream
+
+ subroutine mergeArrays(array1, array2)
+
+ implicit none
+
+ integer, dimension(:), pointer :: array1
+ integer, dimension(:), intent(in) :: array2
+
+ integer :: n1, n2
+ integer, dimension(:), pointer :: newArray
+
+ n1 = size(array1)
+ n2 = size(array2)
+
+ allocate(newArray(n1+n2))
+
+ newArray(1:n1) = array1(:)
+ newArray(n1+1:n1+n2) = array2(:)
+
+ deallocate(array1)
+ array1 => newArray
+
+ end subroutine mergeArrays
+
end module mpas_io_streams
</font>
</pre>