!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    io_chombo.f90 is part of AstroBEAR.
!
!    AstroBEAR is free software: you can redistribute it and/or modify	  
!    it under the terms of the GNU General Public License as published by 
!    the Free Software Foundation, either version 3 of the License, or    
!    (at your option) any later version.
!
!    AstroBEAR is distributed in the hope that it will be useful, 
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with AstroBEAR.  If not, see <http://www.gnu.org/licenses/>.
!
!#########################################################################
!> @defgroup IOChomboGroup IO Chombo Group
!! @ingroup IO
!! @brief Group containing modules used for outputting chombo formatted files

!> @file io_chombo.f90
!! @brief Main file for module IO Chombo

!> @defgroup IOChombo IO Chombo
!! @brief Provides a formaIOCommst-agnostic interface for AstroBEAR to read and write data files.
!! @ingroup IOChomboGroup

!> Provides a formaIOCommst-agnostic interface for AstroBEAR to read and write data files.
!! @ingroup IOChombo
!! @author Brandon D. Shroyer
!! @date 7-9-2010
MODULE IOChombo
#if defined _HDF5
  USE HDF5
  USE HDF5Declarations
  USE ChomboDeclarations
  USE GlobalDeclarations
  USE PhysicsDeclarations
  USE ProcessingDeclarations
  USE TreeDeclarations
  USE DataDeclarations
#if defined _MPI
  USE MessageDeclarations
  USE MpiTransmission
  USE MpiPacking
  USE CommunicationDeclarations
#endif
  USE EllipticDeclarations
  USE DistributionControl
  USE TimeStep, ONLY : olddt, overall_maxspeed,maxcfl
  IMPLICIT NONE
  PRIVATE
  INTEGER :: i_acc_childbox_data_size
  LOGICAL :: lChildDatasetsInitialized
  INTEGER, PARAMETER :: i_node=1, i_q=2, i_aux=3, i_child=4, i_all=5  !this order matters

  PUBLIC MakeChomboFile, ChomboRestartInit, ChomboRestartFinalize, ReadLevel !ChomboRestart


CONTAINS

  !---------------- Frame operations ------------------!

  !> Writes the current data frame to an HDF file using a modified Chombo specification.
  !! @param nframe the number of the sequence frame being written.
  !! @return time spent generating the chombo file.
  SUBROUTINE MakeChomboFile(nframe, text)
    USE ParticleControl, ONLY: lParticles, Particle_WriteData
    INTEGER, INTENT(IN) :: nframe
    CHARACTER(LEN=*) :: text
    CHARACTER(LEN=50) :: s_filename
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: level

    IF (MPI_ID == MASTER .OR. lParallelHDF5) THEN
       WRITE(s_filename, '(A10,A,I5.5,A4)') 'out/chombo', text, nframe, '.hdf'
       CALL CreateChomboHandle(trim(s_filename), chandle, CHOMBO_HANDLE_WRITE)
       CALL WriteGlobalAttributes(chandle, nframe)
       CALL WriteChomboGlobal(chandle)
    END IF
    DO level = -2, FinestLevel
       CALL WriteLevel(chandle, level)
       IF (MPI_ID == MASTER) write(*,*) 'level', level
    END DO
    IF (MPI_ID == MASTER .OR. lParallelHDF5) THEN
       CALL Particle_WriteData(chandle)
       CALL WritePointGravityObjects(chandle)
       CALL WriteOutflowObjects(chandle)
    END IF
    IF (MPI_ID == MASTER .OR. lParallelHDF5) CALL CloseChomboHandle(chandle)

  END SUBROUTINE MakeChomboFile


  !> Begins the process of restarting from the nframe-th chombo file.
  !! @param nframe The frame number to restart from.
  SUBROUTINE ChomboRestartInit(nframe)

    INTEGER :: nframe

    CHARACTER(LEN=I_FILENAME_LENGTH) :: s_filename
    !    TYPE(ChomboHandle), POINTER :: chandle
    ! Retrieve the nframe-th chombo file name.
    WRITE(s_filename, '(A10,I5.5,A4)') 'out/chombo', nframe, '.hdf'

    ! Create a chombo handle to manage the object.
    CALL CreateChomboHandle(s_filename, curr_handle, CHOMBO_HANDLE_READ)
    CALL ReadGlobalAttributes()
  END SUBROUTINE ChomboRestartInit

  SUBROUTINE ChomboRestartFinalize(nframe)
    INTEGER :: nframe
    ! close the chombo handle.
    CALL CloseChomboHandle(curr_handle)
  END SUBROUTINE ChomboRestartFinalize



  !-------------- Level I/O operations --------------!

  SUBROUTINE WriteLevel(chandle, level)
    INTEGER :: level
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER(8) :: datasize(5)
    INTEGER(8) :: mydatasize(5)
    INTEGER(8), ALLOCATABLE, DIMENSION(:,:) :: datasizes
    INTEGER :: DummyInt
#if defined _MPI
    INTEGER :: status(MPI_STATUS_SIZE)
#endif
    INTEGER :: iErr
    TYPE(NodeDefList), POINTER :: nodelist
    INTEGER :: proc
    IF (MPI_ID == Master .OR. lParallelHDF5) THEN
       IF (level >= -1) CALL IO_OpenCurrentLevel(chandle, level)
       CALL ClearChomboHandleOffsets(chandle)
    END IF
    NrInputVars=NrVars+NrDiagnosticVars
    mydatasize=SendChomboData_Precalc(level)
#if defined _MPI
    IF (lParallelHDF5) THEN
       CALL MPI_ALLReduce(mydatasize, datasize, 5, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, iErr)
    ELSE
       CALL MPI_Reduce(mydatasize, datasize, 5, MPI_INTEGER8, MPI_SUM, MASTER, MPI_COMM_WORLD, iErr)
    END IF
#else
    datasize=mydatasize
#endif
    IF (MPI_ID == Master .OR. lParallelHDF5) THEN
       IF (level >= 0) THEN
          CALL AddLevelAttributes(chandle, level)
          CALL AddDataAttributes(chandle)
          CALL Initialize_HDF5_Dataset_Double("data:datatype=0", chandle%level_id, datasize(i_q))
          IF (MaintainAuxArrays) CALL Initialize_HDF5_Dataset_Double("MHD_data", chandle%level_id, datasize(i_aux))
       END IF
       IF (level >= -1) CALL Initialize_HDF5_Dataset_Box("boxes", chandle%level_id, datasize(i_node))
       IF ((level < FinestLevel)) THEN
          IF (level >= -1) THEN
             CALL Initialize_HDF5_Dataset_Int("childbox_count", chandle%level_id, datasize(i_node))
             CALL Initialize_HDF5_Dataset_Box("child_boxes", chandle%level_id, datasize(i_child))
          ELSE
             CALL Add_HDF5_Attribute_Int("root_child_count", chandle%root_group_id, INT(datasize(i_child)))
             CALL Initialize_HDF5_Dataset_Box("root_child_boxes", chandle%root_group_id, datasize(i_child))
          END IF
       END IF
    END IF

    IF (lParallelHDF5) THEN
       CALL WriteNodes(chandle, level, mydatasize)
    ELSE
       ALLOCATE(datasizes(5,1:MPI_NP))
#if defined _MPI
       CALL MPI_Gather(mydatasize, 5, MPI_INTEGER8, datasizes, 5, MPI_INTEGER8, MASTER, MPI_COMM_WORLD, iErr)
#else
       datasizes(:,1)=mydatasize
#endif
       IF (MPI_ID == MASTER) THEN
          nodelist=>nodes(level)%p
          DO WHILE (ASSOCIATED(nodelist))
             CALL WriteNode(chandle, level, nodelist%self)
             nodelist=>nodelist%next
          END DO
          DO proc=1, MPI_NP-1
#if defined _MPI
             CALL MPI_SEND(DummyInt, 1, MPI_INTEGER, proc, TRANSMIT_IO_WORKER_GRIDS, MPI_COMM_WORLD, iErr)
             CALL RecvChomboData(level, chandle, proc, INT(datasizes(:, proc+1)))
#endif
          END DO
#if defined _MPI
       ELSE
          CALL MPI_RECV(DummyInt, 1, MPI_INTEGER, MASTER, TRANSMIT_IO_WORKER_GRIDS, MPI_COMM_WORLD, status, iErr)
          CALL SendChomboData(level, INT(mydatasize))
#endif
       END IF

       DEALLOCATE(datasizes)
    END IF
 
    IF (MPI_ID == MASTER .OR. lParallelHDF5) THEN
       IF (level >= -1)  CALL IO_CloseCurrentLevel(chandle)
    END IF

  END SUBROUTINE WriteLevel

  SUBROUTINE ReadLevel(level)

    INTEGER :: level
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER, ALLOCATABLE, DIMENSION(:) :: displacements, sizes
    INTEGER :: iErr
    INTEGER :: mynodecount
    INTEGER, ALLOCATABLE, DIMENSION(:) :: nodecounts
    INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: patches, mypatches
    TYPE(NodeDefList), POINTER :: nodelist
    INTEGER :: i
    INTEGER :: box_count
    TYPE(NodeBox) :: box
    TYPE(NodeDef), POINTER :: node
    TYPE(InfoDef), POINTER :: info

    IF (MPI_ID == Master .OR. lParallelHDF5) THEN
       chandle=>curr_handle
       IF (level >= -1) CALL IO_OpenCurrentLevel(chandle, level)
       CALL ClearChomboHandleOffsets(chandle)
    END IF

    IF (lParallelHDF5) THEN
       CALL ReadNodes(chandle, level)
    ELSE
       !on master we need to gather all of the boxes on the current level, so we know where to send the data for each node
       !we will pack and send one node at a time so the file access will be sequential
       !Need to perform a gather
       mynodecount=nodecount(nodes(level)%p)    
       ALLOCATE(mypatches(3, 2, mynodecount))
       nodelist=>nodes(level)%p
       i=0
       DO WHILE (ASSOCIATED(nodelist))
          i=i+1
          mypatches(:,:,i)=nodelist%self%box%mGlobal
          nodelist=>nodelist%next
       END DO
       !IF (MPI_ID == MASTER) 
       ALLOCATE(nodecounts(MPI_NP))
#if defined _MPI
       CALL MPI_ALLGATHER(mynodecount, 1, MPI_INTEGER, nodecounts, 1, MPI_INTEGER, MPI_COMM_WORLD, iErr)
#else
       nodecounts(1)=mynodecount
#endif
       ALLOCATE(displacements(MPI_NP), sizes(MPI_NP))
       displacements(1)=0
       DO i=2,MPI_NP
          displacements(i)=displacements(i-1)+nodecounts(i-1)
       END DO
       !       ALLOCATE(patches(3, 2, sum(nodecounts)))
       ALLOCATE(patches(3,2,sum(nodecounts)))
#if defined _MPI
       CALL MPI_GATHERV(mypatches, size(mypatches), MPI_INTEGER, patches, 6*nodecounts, 6*displacements, MPI_INTEGER, MASTER, MPI_COMM_WORLD, iErr)
#else
       patches=mypatches
#endif

       IF (MPI_ID == Master) THEN
          IF (level >= -1) THEN
             box_count = IO_GetDatasetElementCount(chandle%level_id, "boxes")
          ELSE
             box_count = IO_GetDatasetElementCount(chandle%root_group_id, "root_child_boxes")
          END IF
          DO i = 1, box_count
             ! Obtain the next box and store it in the chombo handle structure.
             IF (level >= -1) THEN
                CALL IO_GetBoxFromChomboFile(chandle, "boxes", box%mGlobal)
             ELSE
                box%mGlobal=reshape((/1,1,1,GmX(:)/),(/3,2/))
                !             CALL IO_GetBoxFromChomboFile(chandle, "root_child_boxes", box%mGlobal)
             END IF
             box%mpi_ID=FindMatch(box%mGlobal,patches, displacements)-1
             IF (box%mpi_id == MASTER) THEN
                CALL StrictFindNode(level, box, node, 'readlevel')
                CALL ReadNode(chandle, level, node)
#if defined _MPI             
             ELSE
                ALLOCATE(node)
                CALL NullifyNodeFields(node)
                node%box=box
                IF (level >= 0) THEN
                   ALLOCATE(node%info)
                   CALL NullifyInfoFields(node%info)
                   info=>node%info
                   info%mGlobal=node%box%mGlobal
                   info%mX=info%mGlobal(:,2)-info%mGlobal(:,1)+1
                   ALLOCATE(Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrVars))
                   IF (MaintainAuxArrays) THEN
                      IF (nDim == 2) ALLOCATE(Info%aux(1:Info%mX(1)+1,1:Info%mX(2)+1,1:Info%mX(3),2))
                      IF (nDim == 3) ALLOCATE(Info%aux(1:Info%mX(1)+1,1:Info%mX(2)+1,1:Info%mX(3)+1,3))
                   END IF
                END IF
                CALL ReadNode(chandle, level, node)
                CALL SendNode(level, node)
                CALL DestroyNode(node)
#endif
             END IF
          END DO
#if defined _MPI
          DO i=1,MPI_NP-1
             CALL MPI_SEND(0, 1, MPI_INTEGER, i, TRANSMIT_FRAME_DATA, MPI_COMM_WORLD, iErr)
          END DO
#endif
#if defined _MPI
       ELSE
          CALL RecvNodes(level)
#endif
       END IF
       DEALLOCATE(mypatches)
       DEALLOCATE(nodecounts, displacements, patches)
    END IF
    IF (MPI_ID == MASTER .OR. lParallelHDF5) THEN
       IF (level >= -1)  CALL IO_CloseCurrentLevel(chandle)
    END IF
  END SUBROUTINE ReadLevel
  

  !---------------- Node I/O Operations -------------!

  !Node should have the following:
  !node%info%q
  !node%info%diagnostics
  !node%info%aux
  !node%info%mx
  !node%box%mglobal
  !node%children
  SUBROUTINE WriteNode(chandle, level, node)
    INTEGER :: level
    TYPE(NodeDef), POINTER :: node
    TYPE(InfoDef), POINTER :: info
    INTEGER :: child_count, nq, mx(3), i
    TYPE(NodeDefList), POINTER :: child
    TYPE(ChomboHandle), POINTER :: chandle
    info=>node%info
    child_count=NodeCount(node%children)
    IF (level > -2) THEN 
       CALL Write_Slab_To_Dataset_Box("boxes", chandle%level_id, reshape(node%box%mGlobal, (/6/)) - 1, chandle%box_offset)
       chandle%box_offset=chandle%box_offset+1
       IF (level < FinestLevel) THEN
          CALL Write_Slab_To_Dataset_Int("childbox_count", chandle%level_id,(/ child_count /), chandle%childbox_count_offset)
          chandle%childbox_count_offset = chandle%childbox_count_offset + 1
       END IF
    END IF
    IF (level < FinestLevel) THEN
       child => node%children
       DO WHILE (ASSOCIATED(child))
          IF (level > -2) THEN 
             CALL Write_Slab_To_Dataset_Box("child_boxes", chandle%level_id, reshape(child%self%box%mGlobal, (/6/)) - 1, chandle%childbox_offset)
          ELSE
             CALL Write_Slab_To_Dataset_Box("root_child_boxes", chandle%root_group_id, reshape(child%self%box%mGlobal, (/6/)) - 1, chandle%childbox_offset)
          END IF
          chandle%childbox_offset=chandle%childbox_offset+1             
          child => child%next
       END DO
    END IF
    IF (level >= 0) THEN
       nq=product(Info%mX)*NrVars
       CALL Write_Slab_To_Dataset_Double("data:datatype=0", chandle%level_id, reshape(Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrVars),(/nq/)), chandle%q_offset)
       chandle%q_offset=chandle%q_offset+nq
       IF (NrDiagnosticVars > 0) THEN
          nq=product(Info%mX)*NrDiagnosticVars
          CALL Write_Slab_To_Dataset_Double("data:datatype=0", chandle%level_id, reshape(Info%diagnostics(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrDiagnosticVars),(/nq/)), chandle%q_offset)
          chandle%q_offset=chandle%q_offset+nq
          DEALLOCATE(info%diagnostics)
       END IF
       IF (MaintainAuxArrays) THEN
          mx=Info%mX
          DO i=1,nDim
             mx(i)=mx(i)+1
             nq=product(mx)
             CALL Write_Slab_To_Dataset_Double("MHD_data", chandle%level_id, reshape(Info%aux(1:mX(1),1:mX(2),1:mX(3),i), (/nq/)), chandle%aux_offset)
             chandle%aux_offset=chandle%aux_offset+nq
             mx(i)=mx(i)-1
          END DO
       END IF
    END IF
  END SUBROUTINE WriteNode


  SUBROUTINE WriteNodes(chandle, level, datasize)
    INTEGER :: level
    TYPE(NodeDef), POINTER :: node
    TYPE(InfoDef), POINTER :: info
    INTEGER :: child_count, nq, mx(3), i
    TYPE(NodeDefList), POINTER :: child
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER(8), DIMENSION(5) :: datasize
    INTEGER(8), DIMENSION(5) :: offsets
    INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: boxes
    INTEGER, DIMENSION(:), ALLOCATABLE :: childboxcount
    INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: childboxes
    REAL(KIND=qPREC), DIMENSION(:), ALLOCATABLE :: qdata
    REAL(KIND=qPREC), DIMENSION(:), ALLOCATABLE :: auxdata
    TYPE(NodeDefList), POINTER :: nodelist
    INTEGER :: iErr

    ALLOCATE(boxes(3,2,0:datasize(i_node)-1))
    ALLOCATE(childboxcount(0:datasize(i_node)-1))
    ALLOCATE(childboxes(3,2,0:datasize(i_child)-1))
    ALLOCATE(qdata(0:datasize(i_q)-1))
    ALLOCATE(auxdata(0:datasize(i_aux)-1))
    nodelist=>nodes(level)%p
    DO WHILE (ASSOCIATED(nodelist))
       node=>nodelist%self
       info=>node%info
       child_count=NodeCount(node%children)
       IF (level > -2) THEN 
          boxes(:,:,chandle%box_offset)=node%box%mGlobal - 1
          chandle%box_offset=chandle%box_offset+1
          IF (level < FinestLevel) THEN
             childboxcount(chandle%childbox_count_offset)=child_count
             chandle%childbox_count_offset=chandle%childbox_count_offset + 1
          END IF
       END IF
       IF (level < FinestLevel) THEN
          child => node%children
          DO WHILE (ASSOCIATED(child))
             childboxes(:,:,chandle%childbox_offset)=child%self%box%mGlobal - 1
             chandle%childbox_offset=chandle%childbox_offset+1
             child => child%next
          END DO
       END IF
       IF (level >= 0) THEN
          nq=product(Info%mX)*NrVars
          qdata(chandle%q_offset:chandle%q_offset+nq-1)=reshape(Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrVars),(/nq/))
          chandle%q_offset=chandle%q_offset+nq
          IF (NrDiagnosticVars > 0) THEN
             nq=product(Info%mX)*NrDiagnosticVars
             qdata(chandle%q_offset:chandle%q_offset+nq-1)=reshape(Info%diagnostics(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrDiagnosticVars),(/nq/))
             chandle%q_offset=chandle%q_offset+nq
             DEALLOCATE(info%diagnostics)
          END IF
          IF (MaintainAuxArrays) THEN
             mx=Info%mX
             DO i=1,nDim
                mx(i)=mx(i)+1
                nq=product(mx)
                auxdata(chandle%aux_offset:chandle%aux_offset+nq-1)=reshape(Info%aux(1:mX(1),1:mX(2),1:mX(3),i), (/nq/))
                chandle%aux_offset=chandle%aux_offset+nq
                mx(i)=mx(i)-1
             END DO
          END IF
       END IF
       nodelist=>nodelist%next
    END DO
#if defined _MPI
    CALL MPI_SCAN(datasize, offsets, 5, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, iErr) !just do i_node, i_q, and i_aux    
#else
    offsets=datasize
#endif
    offsets=offsets-datasize !makes scan exclusive

    chandle%box_offset=offsets(i_node) !not necessary
    chandle%childbox_count_offset=offsets(i_node)
    chandle%q_offset=offsets(i_q)
    chandle%aux_offset=offsets(i_aux)
    chandle%childbox_offset=offsets(i_child)

    IF (level > -2) THEN
       CALL Write_Slab_To_Dataset_Box("boxes", chandle%level_id, reshape(boxes, (/size(boxes)/)), chandle%box_offset)
       IF (level < FinestLevel) THEN
          CALL Write_Slab_To_Dataset_Int("childbox_count", chandle%level_id,childboxcount, chandle%childbox_count_offset)
       END IF
    END IF
    IF (level < FinestLevel) THEN
       IF (level > -2) THEN 
          CALL Write_Slab_To_Dataset_Box("child_boxes", chandle%level_id, reshape(childboxes, (/size(childboxes)/)), chandle%childbox_offset)
       ELSE
          CALL Write_Slab_To_Dataset_Box("root_child_boxes", chandle%root_group_id, reshape(childboxes, (/size(childboxes)/)), chandle%childbox_offset)
       END IF
    END IF
    IF (level >= 0) THEN
       CALL Write_Slab_To_Dataset_Double("data:datatype=0", chandle%level_id, qdata, chandle%q_offset)
       IF (MaintainAuxArrays) THEN
          CALL Write_Slab_To_Dataset_Double("MHD_data", chandle%level_id, auxdata, chandle%aux_offset)
       END IF
    END IF

    DEALLOCATE(boxes)
    DEALLOCATE(childboxcount)
    DEALLOCATE(childboxes)
    DEALLOCATE(qdata)
    DEALLOCATE(auxdata)

  END SUBROUTINE WriteNodes


  SUBROUTINE ReadNodes(chandle, level)
    INTEGER :: level
    TYPE(NodeDef), POINTER :: node
    TYPE(InfoDef), POINTER :: info
    INTEGER :: nq, mx(3), i, j, nd
    TYPE(NodeDef), POINTER :: child
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER, DIMENSION(:), ALLOCATABLE :: childboxcount
    INTEGER, DIMENSION(:), ALLOCATABLE :: childboxes
    INTEGER, DIMENSION(:), ALLOCATABLE :: boxes
    REAL(KIND=qPREC), DIMENSION(:), POINTER :: qdata
    REAL(KIND=qPREC), DIMENSION(:), POINTER :: auxdata
    TYPE(NodeDefList), POINTER :: nodelist, lastnode, newnodelist
    INTEGER(8), DIMENSION(5) :: datasize, offsets
    INTEGER(8) :: nchildren
    INTEGER(8) :: offset
    INTEGER :: iErr
    TYPE(NodeBox) :: child_box
    INTEGER, DIMENSION(:,:,:), POINTER :: childgrids
    TYPE(NodeBox) :: box


    datasize=SendChomboData_Precalc(level)  
#if defined _MPI
    CALL MPI_SCAN(datasize, offsets, 5, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, iErr) !just do i_node, i_q, and i_aux    
#else
    offsets=datasize
#endif
    offsets=offsets-datasize !makes scan exclusive

    chandle%box_offset=offsets(i_node) !not necessary
    chandle%childbox_count_offset=offsets(i_node)
    chandle%q_offset=offsets(i_q)
    chandle%aux_offset=offsets(i_aux)
    ALLOCATE(childboxcount(0:datasize(i_node)-1))
    IF (level == -2) THEN
       IF (datasize(i_node) > 0)  childboxcount = 1
    ELSEIF (level < RestartLevel) THEN
       CALL Read_Slab_From_Dataset_Int("childbox_count", chandle%level_id, childboxcount, chandle%childbox_count_offset)

    END IF

    IF (level < RestartLevel) THEN
       !       datasize(i_child)=sum(childboxcount)
       nchildren=sum(childboxcount)
       datasize(i_child)=nchildren
#if defined _MPI
       CALL MPI_SCAN(nchildren, offset, 1, MPI_INTEGER8, MPI_SUM, MPI_COMM_WORLD, iErr)
#else
       offset=nchildren
#endif
       offset=offset-nchildren
       offsets(i_child)=offset
       chandle%childbox_offset=offsets(i_child)
    END IF
!    write(*,'(2I4,A,5I8,A,5I8)') level, MPI_ID, ':', offsets, ':', datasize
    ALLOCATE(boxes(0:3*2*datasize(i_node)-1))
    ALLOCATE(childboxes(0:3*2*datasize(i_child)-1))
    ALLOCATE(qdata(0:datasize(i_q)-1))
    ALLOCATE(auxdata(0:datasize(i_aux)-1))
    IF (level > -2) THEN
       CALL Read_Slab_From_Dataset_Box("boxes", chandle%level_id, boxes, chandle%box_offset)
    ELSE
       if (datasize(i_node) == 1) THEN
          boxes=reshape(nodes(-2)%p%self%box%mGlobal, (/6/))-1
       end if
    END IF
    IF (level < RestartLevel) THEN
       !IF (size(childboxes) > 0) THEN
       IF (level > -2) THEN 
          CALL Read_Slab_From_Dataset_Box("child_boxes", chandle%level_id, childboxes, chandle%childbox_offset)
       ELSE
          CALL Read_Slab_From_Dataset_Box("root_child_boxes", chandle%root_group_id, childboxes, chandle%childbox_offset)
       END IF
       !END IF
    END IF
    IF (level >= 0) THEN
       CALL Read_Slab_From_Dataset_Double("data:datatype=0", chandle%level_id, qdata, chandle%q_offset)
       IF (MaintainAuxArrays) THEN
          CALL Read_Slab_From_Dataset_Double("MHD_data", chandle%level_id, auxdata, chandle%aux_offset)
       END IF
    END IF





    chandle%box_offset=0
    chandle%childbox_count_offset=0
    chandle%childbox_offset=0
    chandle%q_offset=0
    chandle%aux_offset=0


    !asdf
    ! Need to sort nodes based on order in file - so that children will be sorted as well.



    if (.false.) then
       nullify(newnodelist, lastnode)
       nodelist=>Nodes(level)%p
       DO j=1,datasize(i_node)
          !       node=>nodelist%self
          box%mGlobal=reshape(boxes(6*chandle%box_offset:6*(chandle%box_offset+1)-1)+1, (/3,2/))
          chandle%box_offset=chandle%box_offset+1
          box%MPI_ID=MPI_ID
          NULLIFY(node)
          CALL StrictFindNode(level, box, node, 'ReadNodes')
          CALL AddNodeToList(node, lastnode, newnodelist)
          !       nodelist=>nodelist%next
       END DO
       CALL ClearNodeList(nodes(level)%p)
       nodes(level)%p=>newnodelist
       LastLocalNode(level)%p=>lastnode
    end if

    nodelist=>nodes(level)%p
    DO WHILE (ASSOCIATED(nodelist))
       node=>nodelist%self
       IF (level >= 0) info=>node%info
       !       child_count=NodeCount(node%children)
       !       IF (level > -2) THEN 
       !          boxes(:,:,chandle%box_offset)=node%box%mGlobal - 1
       !          chandle%box_offset=chandle%box_offset+1
       !          IF (level < FinestLevel) THEN
       !             childboxcount(chandle%childbox_count_offset)=child_count
       !             chandle%childbox_count_offset=chandle%childbox_count_offset + 1
       !          END IF
       !       END IF
       IF (level < RestartLevel) THEN
          DO i=1,childboxcount(chandle%childbox_count_offset)
             child_box%mGlobal=reshape(childboxes(6*chandle%childbox_offset:6*(chandle%childbox_offset+1)-1)+1, (/3,2/))
             chandle%childbox_offset=chandle%childbox_offset+1
             child_Box%MPI_ID = MPI_ID
             IF (MPI_NP == 1) THEN
                NULLIFY(child)
                CALL AddNode(level+1, child_box, child)
                CALL AddParent(child, node)
             ELSE
                ALLOCATE(child)
                CALL NullifyNodeFields(child)
                child%box=child_box
             END IF
             CALL AddChild(node, child)
!             write(*,'(I4, A, 6I6)') MPI_ID, 'adding child ', child%box
          END DO
          IF ( MPI_NP == 1 .AND. level >= 0) THEN
             nullify(childgrids)
             CALL BuildChildGrids(node, childgrids)
             CALL AllocChildFixups(node%info, childgrids)
             IF (ASSOCIATED(childgrids)) DEALLOCATE(childgrids)
          END IF
          chandle%childbox_count_offset=chandle%childbox_count_offset+1
       END IF
       IF (level >= 0) THEN
          nq=product(Info%mX)*NrVars
          nd=product(Info%mX)*NrInputVars
          Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1:NrVars)=reshape(qdata(chandle%q_offset:chandle%q_offset+nq-1), (/Info%mX(1), Info%mX(2), Info%mX(3), NrVars/))
          chandle%q_offset=chandle%q_offset+nd
          IF (MaintainAuxArrays) THEN
             mx=Info%mX
             DO i=1,nDim
                mx(i)=mx(i)+1
                nq=product(mx)
                Info%aux(1:mX(1),1:mX(2),1:mX(3),i) = reshape(auxdata(chandle%aux_offset:chandle%aux_offset+nq-1), mx)
                chandle%aux_offset=chandle%aux_offset+nq
                mx(i)=mx(i)-1
             END DO
          END IF
       END IF
       nodelist=>nodelist%next
    END DO

    DEALLOCATE(childboxcount)
    DEALLOCATE(boxes)
    DEALLOCATE(childboxes)
    DEALLOCATE(qdata)
    DEALLOCATE(auxdata)


    ! A.)  Use the order in the chombo file to distribute the grids
    ! Each processor has their list of boxes for level n
    ! 1.)  Numerical integration is done on the boxes per proc using MPI_SCAN
    ! 2.)  Each processor reads in their childounts.
    ! 3.)  Numerical integration is done to get processor offsets for childboxes, qdata, and auxdata using MPI_SCAN
    ! 4.)  childboxes, qdata, and auxdata is read in in parallel
    ! 5.)  childboxes are partitioned and distributed in strict processor order (using MPI_SCAN)?
  END SUBROUTINE ReadNodes



  SUBROUTINE ReadNode(chandle, level, node)
    INTEGER :: level
    TYPE(NodeDef), POINTER :: node
    INTEGER :: child_count
    TYPE(NodeDef), POINTER :: child
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: i
    TYPE(NodeBox) :: child_box
    INTEGER, DIMENSION(:,:,:), POINTER :: childgrids
    IF (level > -2) THEN 
       IF (level < RestartLevel) THEN
          CALL IO_GetChildBoxCountFromChomboFile(chandle, child_count)
       END IF
    ELSE
       child_count=1 !assume just 1 domain
    END IF
    IF (level < RestartLevel) THEN
       DO i=1, child_count
          IF (level > -2) THEN 
             CALL IO_GetBoxFromChomboFile(chandle, "child_boxes", child_box%mGlobal)
          ELSE
             CALL IO_GetBoxFromChomboFile(chandle, "root_child_boxes", child_box%mGlobal)
          END IF
          child_box%MPI_ID = MPI_ID
          IF (MPI_NP == 1) THEN
             NULLIFY(child)
             CALL AddNode(level+1, child_box, child)
             CALL AddParent(child, node)
          ELSE
             ALLOCATE(child)
             CALL NullifyNodeFields(child)
             child%box=child_box
          END IF
          CALL AddChild(node, child)
       END DO
       IF (MPI_NP == 1 .AND. level >= 0) THEN
          nullify(childgrids)
          CALL BuildChildGrids(node, childgrids)
          CALL AllocChildFixups(node%info, childgrids)
          IF (ASSOCIATED(childgrids)) DEALLOCATE(childgrids)
       END IF
    END IF
    IF (level >= 0) THEN
       CALL IO_ReadQFromChomboFile(chandle, node, node%info) !also writes aux
    END IF
  END SUBROUTINE ReadNode



  !----------------- Box IO Routines -------------!


  !> Extract a grid box from an active chombo file.
  !! @param chandle An active chombo handle.
  !! @param box_set_name String describing box set
  !! @param box A 3x2 integer array to be filled with a grid box.
  SUBROUTINE IO_GetBoxFromChomboFile(chandle, box_set_name, box)

    TYPE(ChomboHandle), POINTER :: chandle
    CHARACTER(LEN=*) :: box_set_name
    INTEGER, DIMENSION(3,2) :: box

    INTEGER, DIMENSION(6) :: box_data
    INTEGER(HID_T) :: source_group_id

    ! The only major box dataset that is not associated with some kind of level ID
    ! (including the ones in domain_level) is "root_child_boxes", which holds the children
    ! of the level -2 node.
    IF (box_set_name /= "root_child_boxes") THEN
       source_group_id = chandle%level_id
    ELSE
       source_group_id = chandle%root_group_id
    END IF

    ! Read box data in from chombo file and increment the appropriate offset.
    IF (box_set_name == "boxes") THEN
       CALL Read_Slab_From_Dataset_Box(box_set_name, source_group_id, box_data, chandle%box_offset)
       chandle%box_offset = chandle%box_offset + 1
    ELSE
       CALL Read_Slab_From_Dataset_Box(box_set_name, source_group_id, box_data, chandle%childbox_offset)
       chandle%childbox_offset = chandle%childbox_offset + 1
    END IF

    ! Reorder box data to match output array.
    box(1:3,1) = box_data(1:3) + 1
    box(1:3,2) = box_data(4:6) + 1

  END SUBROUTINE IO_GetBoxFromChomboFile

  !------------------ ChildCount IO Routines ----------------------!



  !> Extract a single child-count (integer element) for a grid from a chombo file.
  !! @param chandle An active chombo handle.
  !! @param child_count An integer to hold the child count.
  SUBROUTINE IO_GetChildBoxCountFromChomboFile(chandle, child_count)

    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: child_count
    INTEGER, DIMENSION(1):: childcount_array


    ! Reads cell-centered data in from chombo file.
    CALL Read_Slab_From_Dataset_Int("childbox_count", chandle%level_id, childcount_array, chandle%childbox_count_offset)
    child_count = childcount_array(1)

    ! Increments chombo handle's cell-centered data offset.
    chandle%childbox_count_offset = chandle%childbox_count_offset + 1

  END SUBROUTINE IO_GetChildBoxCountFromChomboFile



  !------------------ q and aux routines for I/O -----------------------!

  SUBROUTINE IO_ReadQFromChomboFile(chandle, node, info)
    TYPE(ChomboHandle), POINTER :: chandle
    TYPE(NodeDef), POINTER :: node
    TYPE(InfoDef), POINTER :: info
    INTEGER :: mx, my, mz, nax, nay, naz, nq, nd, ncells
    REAL(KIND=qPREC), DIMENSION(:), POINTER :: dbla_chombo_qvar, dbla_chombo_aux

    ! Obtain the grid dimensions from the coordinate box.
    mx = node%box%mGlobal(1,2) - node%box%mGlobal(1,1) + 1
    my = node%box%mGlobal(2,2) - node%box%mGlobal(2,1) + 1
    mz = node%box%mGlobal(3,2) - node%box%mGlobal(3,1) + 1
    ncells=mx*my*mz
    nq=ncells*NrVars
    nd=ncells*NrInputVars

    ! Allocate memory for the current node's grid.
    ALLOCATE(dbla_chombo_qvar(nq))

    ! Read the data array to the Chombo file using the supplied offset.
    CALL Read_Slab_From_Dataset_Double("data:datatype=0", chandle%level_id, dbla_chombo_qvar, chandle%q_offset)
    Info%q(1:mx, 1:my, 1:mz, 1:NrVars)=reshape(dbla_chombo_qvar, (/mx,my,mz,NrVars/))


    ! Adjust the offset by the size of the grid just written.
    chandle%q_offset = chandle%q_offset +nd
    DEALLOCATE(dbla_chombo_qvar)

    ! If MHD is running, then do the same thing for the Info%aux array.
    IF (MaintainAuxArrays) THEN 
       nax=(mx+1)*my*mz
       nay=(my+1)*mx*mz
       naz=(mz+1)*mx*my
       IF (nDim <= 2) naz=0
       ! Allocate memory for the current node's aux grid.
       ALLOCATE(dbla_chombo_aux(nax+nay+naz))                  

       ! Write the data array to the Chombo file using the supplied offset.
       CALL Read_Slab_From_Dataset_Double("MHD_data", chandle%level_id, dbla_chombo_aux, chandle%aux_offset)
       ! Adjust the offset by the size of the grid just written.
       chandle%aux_offset = chandle%aux_offset + nax+nay+naz

       Info%aux(1:mx+1,1:my,1:mz,1)=RESHAPE(dbla_chombo_aux(1:nax), (/mx+1,my,mz/))
       Info%aux(1:mx,1:my+1,1:mz,2)=RESHAPE(dbla_chombo_aux(nax+1:nax+nay), (/mx,my+1,mz/))
       IF (nDim >= 3) Info%aux(1:mx,1:my,1:mz+1,3)=RESHAPE(dbla_chombo_aux(nax+nay+1:nax+nay+naz), (/mx,my,mz+1/))

       ! Release array memory so that it can be used by the next grid.
       DEALLOCATE(dbla_chombo_aux)
    END IF         ! (if MaintainAuxArrays).
  END SUBROUTINE IO_ReadQFromChomboFile




  !------------------ Communication Routines for output ------------------!

#if defined _MPI
  SUBROUTINE RecvChomboData(n, chandle, proc, datasize)
    INTEGER :: n
    INTEGER :: proc
    TYPE(ChomboHandle), POINTER :: chandle
    TYPE(PackedMessage), POINTER :: message
    TYPE(InfoDef), POINTER :: info
    TYPE(NodeDef), POINTER :: node
    INTEGER, DIMENSION(:) :: datasize
    INTEGER :: nchildren, i, mglobal(3,2)
    TYPE(NodeDef), POINTER :: child
    CALL CreatePackedMessage(n, proc, TRANSMIT_FRAME_DATA, STAGE_RECV, message, datasize(i_all))
    CALL ReceivePackedMessage(message)
    DO WHILE(GetNextBox(message, mGlobal))
       ALLOCATE(node)
       CALL NullifyNodeFields(node)
       node%box%mglobal=mglobal
       ALLOCATE(node%info)
       CALL NullifyInfoFields(node%info)
       info=>node%info
       info%mx=node%box%mGlobal(:,2)-node%box%mGlobal(:,1)+1
       IF (n >= 0) THEN
          ALLOCATE(Info%q(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3), 1:NrVars))
          CALL UnPackData(message, info%q)
          IF (NrDiagnosticVars > 0) THEN
             ALLOCATE(info%diagnostics(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3), 1:NrDiagnosticVars))
             CALL UnPackData(message, info%diagnostics)
          END IF
          IF (MaintainAuxArrays) THEN
             IF (nDim == 2) THEN 
                ALLOCATE(Info%aux(1:Info%mX(1)+1, 1:Info%mX(2)+1, 1:Info%mX(3),2))
             ELSEIF (nDim == 3) THEN
                ALLOCATE(Info%aux(1:Info%mX(1)+1, 1:Info%mX(2)+1, 1:Info%mX(3)+1,3))
             END IF
             CALL UnPackData(message, Info%aux(1:Info%mX(1)+1, 1:Info%mX(2), 1:Info%mX(3), 1))
             CALL UnPackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2)+1, 1:Info%mX(3), 2))
             IF (nDim == 3) CALL UnPackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3)+1, 3))
          END IF

       END IF
       CALL UnPackData(message, nChildren)
       DO i=1,nchildren
          ALLOCATE(child)
          CALL NullifyNodeFields(child)
          CALL AddChild(node, child)
          CALL UNPackData(message, child%box%mglobal)
       END DO
       CALL WriteNode(chandle, n,node)
       !       CALL DestroyInfo(node%info)
       CALL DestroyNode(node)
    END DO
    CALL DestroyPackedMessage(message)
  END SUBROUTINE RecvChomboData

  SUBROUTINE SendChomboData(n, datasize)
    INTEGER :: n
    TYPE(PackedMessage), POINTER :: message
    TYPE(NodeDefList), POINTER :: nodelist, childlist
    TYPE(InfoDef), POINTER :: info
    INTEGER, DIMENSION(:) :: datasize
    TYPE(NodeDef), POINTER :: node
    CALL CreatePackedMessage(n, Master, TRANSMIT_FRAME_DATA, STAGE_SEND, message, datasize(i_all))
    nodelist=>Nodes(n)%p
    DO WHILE (ASSOCIATED(nodelist))
       node=>nodelist%self
       info=>node%info
       CALL PackData(message, node%box%mglobal)
       IF (n >= 0) THEN
          CALL PackData(message, info%q(1:Info%mX(1), 1:Info%mX(2), 1:Info%mx(3), 1:NrVars))
          IF (NrDiagnosticVars > 0) THEN
             CALL PackData(message, info%diagnostics(1:Info%mX(1), 1:Info%mX(2), 1:Info%mx(3), 1:NrDiagnosticVars))
             DEALLOCATE(Info%diagnostics) !this needs to be done on workers as part of IO process.
          END IF
          IF (MaintainAuxArrays) THEN
             CALL PackData(message, Info%aux(1:Info%mX(1)+1, 1:Info%mX(2), 1:Info%mX(3), 1))
             CALL PackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2)+1, 1:Info%mX(3), 2))
             IF (nDim == 3) CALL PackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3)+1, 3))
          END IF
       END IF
       CALL PackData(message, NodeCount(node%children))
       childlist=>node%children
       DO WHILE (ASSOCIATED(childlist))
          CALL PackData(message, childlist%self%box%mglobal)
          childlist=>childlist%next
       END DO
       nodelist=>nodelist%next
    END DO
    CALL PackData(message, TERMINATIONBOX)
    CALL ClosePackedMessage(message)
    CALL DestroyPackedMessage(message)
  END SUBROUTINE SendChomboData
#endif

  FUNCTION SendChomboData_Precalc(n)
    INTEGER, DIMENSION(5) :: SendChomboData_Precalc
    INTEGER :: n
    INTEGER :: datasize(5)
    TYPE(InfoDef), POINTER :: info
    TYPE(NodeDefList), POINTER :: nodelist, childlist
    TYPE(NodeDef), POINTER :: node
    datasize=0
    nodelist=>Nodes(n)%p
    DO WHILE (ASSOCIATED(nodelist))
       node=>nodelist%self
       info=>node%info
       datasize(i_node)=datasize(i_node)+1
       IF (n >= 0) THEN
          datasize(i_q)=datasize(i_q)+product(info%mx)*NrInputVars
          IF (MaintainAuxArrays) THEN
             datasize(i_aux)=datasize(i_aux)+(Info%mX(1)*Info%mX(2)*Info%mX(3))*nDim
             datasize(i_aux)=datasize(i_aux)+(Info%mX(1)+Info%mX(2))*Info%mX(3)
             IF (nDim >=3) datasize(i_aux)=datasize(i_aux)+(Info%mX(1)*Info%mX(2))

!             datasize(i_aux)=datasize(i_aux)+(Info%mX(1)+1)*(Info%mX(2)*Info%mX(3))
!             datasize(i_aux)=datasize(i_aux)+(Info%mX(2)+1)*(Info%mX(1)*Info%mX(3))
!             IF (nDim >=3) datasize(i_aux)=datasize(i_aux)+(Info%mX(3)+1)*(Info%mX(1)*Info%mX(2))
          END IF
       END IF
       childlist=>node%children
       DO WHILE (ASSOCIATED(childlist))
          datasize(i_child)=datasize(i_child)+1
          childlist=>childlist%next
       END DO
       nodelist=>nodelist%next
    END DO
    datasize(i_all)=6*PACK_INTEGER_SIZE + & !for termination box
         datasize(i_node)*6*PACK_INTEGER_SIZE + & !for mglobals
         datasize(i_node)*1*PACK_INTEGER_SIZE + & !for nchildren
         datasize(i_child)*6*PACK_INTEGER_SIZE + & !for child mglobals
         datasize(i_q)*PACK_DOUBLE_SIZE + & !for grid data
         datasize(i_aux)*PACK_DOUBLE_SIZE !for aux data
    SendChomboData_Precalc=datasize
  END FUNCTION SendChomboData_Precalc

  !-------------- Communication Routines for input -------------------!

  FUNCTION SendNode_Precalc(n, node)
    INTEGER :: SendNode_Precalc
    INTEGER :: n
    INTEGER :: datasize
    TYPE(InfoDef), POINTER :: info
    TYPE(NodeDefList), POINTER :: childlist
    TYPE(NodeDef), POINTER :: node
    datasize=0
    info=>node%info
    datasize=datasize+6*PACK_INTEGER_SIZE
 
    IF (n >= 0) THEN
       datasize=datasize+product(info%mx)*NrVars*PACK_DOUBLE_SIZE
       IF (MaintainAuxArrays) THEN
          datasize=datasize+(Info%mX(1)*Info%mX(2)*Info%mX(3))*nDim*PACK_DOUBLE_SIZE
          datasize=datasize+(Info%mX(1)+Info%mX(2))*Info%mX(3)*PACK_DOUBLE_SIZE
          IF (nDim >=3) datasize=datasize+(Info%mX(1)*Info%mX(2))*PACK_DOUBLE_SIZE
       END IF
    END IF
    IF (n < MaxLevel) THEN
       datasize=datasize+PACK_INTEGER_SIZE
       childlist=>node%children
       DO WHILE (ASSOCIATED(childlist))
          datasize=datasize+6*PACK_INTEGER_SIZE
          childlist=>childlist%next
       END DO
    END IF
    SendNode_Precalc=datasize
  END FUNCTION SendNode_Precalc


#if defined _MPI
  SUBROUTINE SendNode(level, node)
    INTEGER :: level
    TYPE(NodeDef), POINTER :: node
    TYPE(PackedMessage), POINTER :: message
    TYPE(NodeDefList), POINTER :: childlist
    TYPE(InfoDef), POINTER :: info
    INTEGER :: datasize
    INTEGER :: iErr
    datasize=INT(SendNode_Precalc(level, node))
    CALL MPI_SEND(datasize, 1, MPI_INTEGER, node%box%MPI_ID, TRANSMIT_FRAME_DATA, MPI_COMM_WORLD, iErr)
    CALL CreatePackedMessage(level, node%box%MPI_ID, TRANSMIT_FRAME_DATA, STAGE_SEND, message, datasize)
    info=>node%info
    CALL PackData(message, node%box%mglobal)
    IF (level >= 0) THEN
       CALL PackData(message, info%q(1:Info%mX(1), 1:Info%mX(2), 1:Info%mx(3), 1:NrVars))
       IF (MaintainAuxArrays) THEN
          CALL PackData(message, Info%aux(1:Info%mX(1)+1, 1:Info%mX(2), 1:Info%mX(3), 1))
          CALL PackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2)+1, 1:Info%mX(3), 2))
          IF (nDim == 3) CALL PackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3)+1, 3))
       END IF
    END IF
    IF (level < MaxLevel) THEN
       CALL PackData(message, NodeCount(node%children))
       childlist=>node%children
       DO WHILE (ASSOCIATED(childlist))
          CALL PackData(message, childlist%self%box%mglobal)
          childlist=>childlist%next
       END DO
    END IF
    CALL ClosePackedMessage(message)
    CALL DestroyPackedMessage(message)
  END SUBROUTINE SendNode

  SUBROUTINE RecvNodes(level)
    INTEGER :: level
    TYPE(PackedMessage), POINTER :: message
    TYPE(InfoDef), POINTER :: info
    TYPE(NodeDef), POINTER :: node
    INTEGER :: datasize
    INTEGER :: nchildren, i
    TYPE(NodeDef), POINTER :: child
    INTEGER :: status(MPI_STATUS_SIZE)
    INTEGER :: iErr
    TYPE(NodeBox) :: Box
    DO
       CALL MPI_RECV(datasize, 1, MPI_INTEGER, Master, TRANSMIT_FRAME_DATA, MPI_COMM_WORLD, status, iErr)
       IF (datasize == 0) EXIT
       CALL CreatePackedMessage(level, Master, TRANSMIT_FRAME_DATA, STAGE_RECV, message, INT(datasize))
       CALL ReceivePackedMessage(message)
       CALL UnPackData(message, box%mglobal)
       box%MPI_ID=MPI_ID
       CALL StrictFindNode(level, box, node,'recvnodes')
       IF (level >= 0) THEN
          info=>node%info
          CALL UnPackData(message, info%q(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3), 1:NrVars))
          IF (MaintainAuxArrays) THEN
             CALL UnPackData(message, Info%aux(1:Info%mX(1)+1, 1:Info%mX(2), 1:Info%mX(3), 1))
             CALL UnPackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2)+1, 1:Info%mX(3), 2))
             IF (nDim == 3) CALL UnPackData(message, Info%aux(1:Info%mX(1), 1:Info%mX(2), 1:Info%mX(3)+1, 3))
          END IF
       END IF
       IF (level < MaxLevel) THEN
          CALL UnPackData(message, nChildren)
          DO i=1,nchildren
             ALLOCATE(child)
             CALL NullifyNodeFields(child)
             CALL AddChild(node, child)
             CALL UNPackData(message, child%box%mglobal)
          END DO
       END IF
       CALL DestroyPackedMessage(message)
    END DO
  END SUBROUTINE RecvNodes

#endif


  !-------------- Global Attribute Routines ------------------!

  SUBROUTINE WriteGlobalAttributes(chandle, nframe)
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: nframe
    CHARACTER(LEN = I_MAX_CNAME_LENGTH) :: s_component_tag
    CHARACTER(LEN = I_MAX_CNAME_LENGTH), DIMENSION(0:I_MAX_COMPONENTS) :: sa_component_names
    TYPE(DiagnosticListDef), POINTER :: Diagnostic
    INTEGER :: NrCons, i
    INTEGER (HSIZE_T):: i_index        ! Used for any incidental loops along the way

    ! time attribute (uses Info's current time).
    CALL Add_HDF5_Attribute_Double("time", chandle%root_group_id, levels(0)%tnow)
    CALL Add_HDF5_Attribute_Double("olddt", chandle%root_group_id, olddt)
    CALL Add_HDF5_Attribute_Double("overall_maxspeed", chandle%root_group_id, overall_maxspeed)
    CALL Add_HDF5_Attribute_Double("maxcfl", chandle%root_group_id, maxcfl)

    ! iteration attribute (uses Info's current frame number).
    CALL Add_HDF5_Attribute_Int("iteration", chandle%root_group_id, nFrame)

    ! num_components attribute (use Info's number of output variables).
    CALL Add_HDF5_Attribute_Int("num_components", chandle%root_group_id, NrVars+NrDiagnosticVars)

    ! node centering (fixes contour maps in visit)
    CALL Add_HDF5_Attribute_Int("data_centering", chandle%root_group_id, 0) ! cell centered

!!! Storing the domain boundaries in the Chombo file.
    CALL Add_HDF5_Attribute_FloatVector("lower_bound", chandle%root_group_id, GxBounds(:,1))
    CALL Add_HDF5_Attribute_FloatVector("upper_bound", chandle%root_group_id, GxBounds(:,2))

    ! num_levels attribute (Chombo format requires that this be 1 greater the maximum
    ! level--essentially so that it acts as a level-index.
    CALL Add_HDF5_Attribute_Int("num_levels", chandle%root_group_id, FinestLevel + 1)

    ! max_level attribute (use Info's current FinestLevel value).
    CALL Add_HDF5_Attribute_Int("max_level", chandle%root_group_id, FinestLevel)

!!!This pro_lo attribute causes pseudocolor slices to not work with multiple root patches in visit
    !CALL Add_HDF5_Attribute_FloatVector("prob_lo", chandle%root_group_id, GxBounds(:,1))

!!!!!    The scaling factors below are not part of the Chombo HDF specification, but they're
!!!!!    useful enough that we're including them anyway.

    CALL Add_HDF5_Attribute_Double("n_density_scale", chandle%root_group_id, nScale)
    CALL Add_HDF5_Attribute_Double("rho_scale", chandle%root_group_id, rScale)
    CALL Add_HDF5_Attribute_Double("length_scale", chandle%root_group_id, lScale)
    CALL Add_HDF5_Attribute_Double("velocity_scale", chandle%root_group_id, velScale)
    CALL Add_HDF5_Attribute_Double("pressure_scale", chandle%root_group_id, pScale)
    CALL Add_HDF5_Attribute_Double("temperature_scale", chandle%root_group_id, TempScale)
    IF (ivx /= 0) CALL Add_HDF5_Attribute_Double("px_scale", chandle%root_group_id, rScale * velScale)
    IF (ivy /= 0) CALL Add_HDF5_Attribute_Double("py_scale", chandle%root_group_id, rScale * velScale)
    IF (ivz /= 0) CALL Add_HDF5_Attribute_Double("pz_scale", chandle%root_group_id, rScale * velScale)
    IF (iE .ne. 0) CALL Add_HDF5_Attribute_Double("E_scale", chandle%root_group_id, pScale)
    IF (iBx .ne. 0) CALL Add_HDF5_Attribute_Double("Bx_scale", chandle%root_group_id, pScale)
    IF (iBy .ne. 0) CALL Add_HDF5_Attribute_Double("By_scale", chandle%root_group_id, pScale)
    IF (iBz .ne. 0) CALL Add_HDF5_Attribute_Double("Bz_scale", chandle%root_group_id, pScale)

    !Uses vrar_index to keep track of what variables are in what slots
    NrCons = NrHydroVars - NrTracerVars

    ! Create tags for conservative variables
    DO i=1, NrCons
       sa_component_names(i-1)= trim(FieldName(i))
    END DO

    ! Creates tags for tracer variables.
    DO i = 1, NrTracerVars 
       sa_component_names(nTracerLo-1+i-1) = trim(TracerName(i))
    END DO

    ! Creates tags for Elliptic variables.
    DO i = 1, NrEllipticVars 
       sa_component_names(nEllipticLo-1+i-1) = trim(EllipticName(i))       
    END DO

    ! Creates tags for Diagnostic variables
    Diagnostic=>FirstDiagnostic
    i=1
    DO WHILE (ASSOCIATED(Diagnostic))
       sa_component_names(NrHydroVars+NrEllipticVars+i-1) = trim(Diagnostic%Field%Name)
       i=i+1
       Diagnostic=>Diagnostic%next
    END DO


    ! This part is separated out so that changing the assembly of the components
    ! is easy, but the addition of attributes may eventually have to be done
    ! inline in order to reduce overhead.
    !
    ! Construct new component_n attributes for each sa_component_name(n).
    DO i_index = 0, NrVars + NrDiagnosticVars - 1

       ! Create component tag (attribute name).
       IF (i_index <= 9) THEN
          WRITE(s_component_tag,'(A10,I1.1)') 'component_', i_index
       ELSE
          WRITE(s_component_tag,'(A10,I2.2)') 'component_', i_index
       END IF

       CALL Add_HDF5_Attribute_String(s_component_tag, chandle%root_group_id, TRIM(sa_component_names(i_index)))
    END DO

  END SUBROUTINE WriteGlobalAttributes

  SUBROUTINE ReadGlobalAttributes()

    current_time = Read_HDF5_Attribute_Double("time", curr_handle%root_group_id)
    olddt = Read_HDF5_Attribute_Double("olddt", curr_handle%root_group_id)
    overall_maxspeed = Read_HDF5_Attribute_Double("overall_maxspeed", curr_handle%root_group_id)
    maxcfl=Read_HDF5_Attribute_Double("maxcfl", curr_handle%root_group_id)
    restartlevel=Read_HDF5_Attribute_Int("max_level", curr_handle%root_group_id)
    levels(:)%tnow = current_time
    levels(:)%dt = 0.d0
    levels(:)%step = 1
    RestartLevel=min(MaxLevel,restartlevel)
    NrInputVars = READ_HDF5_Attribute_Int("num_components", curr_handle%root_group_id)
    IF (NrInputVars < NrVars) THEN
       IF (MPI_ID == 0) PRINT*, 'Error - not enough variables in chombo file fro restart'
       STOP
    END IF
  END SUBROUTINE ReadGlobalAttributes


  SUBROUTINE WriteChomboGlobal(chandle)
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: i_Err
    INTEGER(HID_T) :: hid_chombo_global_group_id

    ! Create Chombo_global group.
    CALL h5gcreate_f(chandle%file_handle, "Chombo_global", hid_chombo_global_group_id, i_err)

    ! Add testReal attribute to Chombo_global group (value is a literal).
    CALL Add_HDF5_Attribute_Double("testReal", hid_chombo_global_group_id, DBL_TEST_REAL_DAT) 

    ! Add SpaceDim attribute to Chombo_global group.  We can just use Info%nDim

    IF (nDim == 1) THEN !Chombo reader for visit requires nDim to be at least 2.
       CALL Add_HDF5_Attribute_Int("SpaceDim", hid_chombo_global_group_id, 2)
    ELSE
       CALL Add_HDF5_Attribute_Int("SpaceDim", hid_chombo_global_group_id, nDim)
    END IF
    ! Close the Chombo_global group.
    CALL h5gclose_f(hid_chombo_global_group_id, i_err)

  END SUBROUTINE WriteChomboGlobal


  !------------- Level Attribute Routines -----------------!

  SUBROUTINE AddLevelAttributes(chandle, i_level)
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: i_level
    INTEGER, DIMENSION(I_CHOMBO_BOX_SIZE) :: ia_box_global        ! array for holding box data.

    ! Create the dx attribute
    CALL Add_HDF5_Attribute_Double("dx", chandle%level_id, levels(i_level)%dx)
    ! Create the anisotropic attribute (optional chomobo attrib for anisotropic spacing)

    CALL Add_HDF5_Attribute_FloatVector("anisotropic", chandle%level_id, (/1d0,1d0,1d0/))

    ! Create the ref_ratio attribute.
    CALL Add_HDF5_Attribute_Int("ref_ratio", chandle%level_id, levels(i_level)%CoarsenRatio)

    ! Create prob_domain attribute using a box object.
    ia_box_global=0
    ia_box_global(3+1:3+nDim) = Gmx(1:nDim) * product(levels(0:i_level-1)%coarsenratio)-1
    CALL Add_HDF5_Attribute_Box("prob_domain", hid_box_id, chandle%level_id, ia_box_global)
  END SUBROUTINE AddLevelAttributes

  SUBROUTINE AddDataAttributes(chandle)
    TYPE(ChomboHandle), POINTER :: chandle
    INTEGER :: i_err
    INTEGER :: hid_data_attributes_group_id
    ! Create the data attributes group.
    CALL h5gcreate_f(chandle%level_id, "data_attributes", hid_data_attributes_group_id, i_err)

    ! Adds the output_ghost attribute using the literal intvector [0, 0, 0].
    CALL Add_HDF5_Attribute_IntVector("outputGhost", hid_intvect_id, hid_data_attributes_group_id, (/0,0,0/))

    ! Adds the comps attribute (the same as num_components)
    CALL Add_HDF5_Attribute_Int("comps", hid_data_attributes_group_id, NrVars+NrDiagnosticVars)

    ! Adds the objectType attribute.
    CALL Add_HDF5_Attribute_String("objectType", hid_data_attributes_group_id, "FArrayBox")

    ! Close data attributes group.
    CALL h5gclose_f(hid_data_attributes_group_id, i_err)
  END SUBROUTINE AddDataAttributes


  !-------------- Miscellaneous I/O routines -----------------!


  FUNCTION FindMatch(mGlobal, patches, displacements)
    INTEGER :: mGlobal(:,:), patches(:,:,:), displacements(:), FindMatch
    INTEGER :: i, j, n
    n=size(displacements)
    FindMatch=-1
    DO i=1,size(patches, 3)
       IF (ALL(patches(:,:,i)==mGlobal)) THEN
          FindMatch=n
          DO j=1, n
             IF (displacements(j) >= i) THEN
                FindMatch=j-1
                EXIT
             END IF
          END DO
          RETURN
       END IF
    END DO
    IF (FindMatch==-1) THEN
       write(*,*) 'Error did not find a match for'
       write(*,'(6I6)') mGlobal
       write(*,'(6I6)') patches
       write(*,'(100I5)') displacements
       STOP
    END IF
  END FUNCTION FindMatch


  SUBROUTINE ClearNode(proc_id, info, node, nodelist)
    INTEGER :: proc_id
    TYPE(InfoDef), POINTER :: info
    TYPE(NodeDef), POINTER :: node
    TYPE(NodeDefList), POINTER :: nodelist

    ! Destroy the structures created from remote data.
    IF (proc_id > Master) THEN
       ! Deallocate the cell-centered data array.
       DEALLOCATE(Info%q)
       NULLIFY(Info%q)
       ! Deallocate the aux array imported from the remote processor
       IF (MaintainAuxArrays) THEN
          DEALLOCATE(Info%aux)
          NULLIFY(Info%aux)
       END IF
       DEALLOCATE(Info)
       NULLIFY(Info)
       IF (ASSOCIATED(node%children))  CALL DestroyNodeList(node%children)
       DEALLOCATE(node)
       NULLIFY(node)
    ELSE
       nodelist => nodelist%next
    END IF
  END SUBROUTINE ClearNode



  !> Creates a group for particle data within the Chombo file and populates it with the contents of the SinkParticles list.
  !! @param chandle A chombo file handle.
  SUBROUTINE WritePointGravityObjects(chandle)

    USE PointGravitySrc, ONLY: FirstPointGravityObj, PointGravityDef, PointGravity_CountObjects, PointGravity_InitChomboDatasets, PointGravity_WriteObjectToChombo

    TYPE(ChomboHandle), POINTER :: chandle

    TYPE(PointGravityDef), POINTER :: gravity_object
    INTEGER :: i_err
    INTEGER :: obj_count


    ! Open a new particles group within the root level
    CALL h5gcreate_f(chandle%root_group_id, "point_gravity_objects", chandle%source_group_id, i_err)

    chandle%source_offset = 0

    ! Retrieve the number of point gravity objects from the object list.
    obj_count = PointGravity_CountObjects()
    ! Store the number of objects.
    CALL Add_HDF5_Attribute_Int("num_objects", chandle%source_group_id, obj_count)

    ! If there are any gravity objects, then create datasets to store their data.  Otherwise,
    ! don't create the datasets (since there's no need to take up file space).
    IF (obj_count > 0) THEN

       CALL PointGravity_InitChomboDatasets(chandle, obj_count)

       gravity_object => FirstPointGravityObj

       ! Loop over the list of point gravity objects and write their data to the new component datasets.
       DO WHILE (ASSOCIATED(gravity_object))
          IF (gravity_object%lStoreInChombo) THEN
             CALL PointGravity_WriteObjectToChombo(chandle, gravity_object)
          END IF
          gravity_object => gravity_object%next

       END DO

    END IF

    ! Close the source term group.
    CALL h5gclose_f(chandle%source_group_id, i_err)

  END SUBROUTINE WritePointGravityObjects

  !> Creates a group for outflow data within the Chombo file and populates it with the contents of the Outflow list.
  !! @param chandle A chombo file handle.
  SUBROUTINE WriteOutflowObjects(chandle)

    USE Outflows, ONLY: OutflowDef, &
         Outflow_InitChomboDatasets, Outflow_WriteObjectToChombo
    USE ObjectControl

    TYPE(ChomboHandle), POINTER :: chandle
    TYPE(ObjectDef), POINTER :: object
    INTEGER :: i_err
    INTEGER :: obj_count


    ! Open a new outflows group within the root level
    CALL h5gcreate_f(chandle%root_group_id, "outflow_objects", chandle%source_group_id, i_err)

    chandle%source_offset = 0

    ! Retrieve the number of point gravity objects from the object list.
    obj_count = ObjectsCount(outflowobj)

    ! Store the number of objects.
    CALL Add_HDF5_Attribute_Int("num_objects", chandle%source_group_id, obj_count)

    ! If there are any outflow objects, then create datasets to store their data.  Otherwise,
    ! don't create the datasets (since there's no need to take up file space).
    IF (obj_count > 0) THEN

       CALL Outflow_InitChomboDatasets(chandle, obj_count)

       Object => ListHead

       ! Loop over the list of outflow objects and write their data to the new component datasets.
       DO WHILE(ASSOCIATED(Object))
          IF (Object%type == outflowobj) THEN
             pOutflow = transfer(Object%storage,pOutflow)
             IF (pOutflow%ptr%lStoreInChombo) THEN
                CALL Outflow_WriteObjectToChombo(chandle, pOutflow%ptr)
             END IF

          END IF
          Object => Object%next
       END DO

    END IF

    ! Close the source term group.
    CALL h5gclose_f(chandle%source_group_id, i_err)

  END SUBROUTINE WriteOutflowObjects
#endif

END MODULE IOChombo

