Merge branch 'fhi-level' into 'master'

Fhi level

See merge request roy/pFoF_dev!6

Makefile

 include Make.inc
 
-all: pfof_snap pfof_cone conepartcreator conegravcreator conemapper
+all: pfof_snap.exe pfof_cone.exe amr2cell_hdf5.exe conepartcreator.exe conegravcreator.exe conemapper.exe
 
-pfof_snap: 
+pfof_snap.exe: 
 	cd pfof_snap/src ; make
 
-pfof_cone:
+pfof_cone.exe:
 	cd pfof_cone/src ; make
 
-conepartcreator:
+amr2cell_hdf5.exe:
+	cd amr2cell_hdf5/src ; make
+
+conepartcreator.exe:
 	cd tools/conepartcreator/src ; make
 
-conegravcreator:
+conegravcreator.exe:
 	cd tools/conegravcreator/src ; make
 
-conemapper:
+conemapper.exe:
 	cd tools/conemapper/src ; make
 
 clean:
 	cd pfof_snap/src ; make clean
 	cd pfof_cone/src ; make clean
+	cd amr2cell_hdf5/src ; make clean
 	cd tools/conepartcreator/src ; make clean
 	cd tools/conegravcreator/src ; make clean
 	cd tools/conemapper/src ; make clean

amr2cell_hdf5/src/Makefile

@@ -16,12 +16,12 @@ SRC = char_utils_m.f90			\
 	index_m.f90			\
 	ramses_info_m.f90 		\
 	amr2cell_hdf5_parameters_m.f90	\
+	ramses_hydro_descriptor_m.f90	\
 	cell_m.f90			\
 	cells_array_m.f90		\
 	cube_m.f90			\
 	cubes_array_m.f90		\
 	ramses_hydro_header_m.f90	\
-	ramses_hydro_descriptor_m.f90	\
 	ramses_hydro_data_m.f90		\
 	ramses_amr_header_m.f90 	\
 	ramses_amr_data_m.f90		\

amr2cell_hdf5/src/Makefile.read_cell_hdf5

+include ../../Make.inc
+
+EXE = read_cell_hdf5.exe
+
+COMMONDIR = ../../common/src/
+OBJDIR=../../obj/
+MODDIR=../../mod/
+VPATH= $(COMMONDIR):.
+
+SRC = char_utils_m.f90			\
+	error_handling_m.f90		\
+	constants_m.f90			\
+	mpi_process_m.f90		\
+	timer_m.f90 			\
+    mpi_utils_m.f90     \
+	shared_variables_m.f90 		\
+	index_m.f90			\
+	ramses_info_m.f90 		\
+	ramses_hydro_descriptor_m.f90	\
+	amr2cell_hdf5_parameters_m.f90	\
+	cell_m.f90			\
+	cells_array_m.f90		\
+	cube_m.f90			\
+	cubes_array_m.f90		\
+	ramses_hydro_header_m.f90	\
+	ramses_hydro_data_m.f90		\
+	ramses_amr_header_m.f90 	\
+	ramses_amr_data_m.f90		\
+	ramses_output_m.f90		\
+	ramses_output_to_cells_m.f90	\
+	read_cell_hdf5.f90
+
+OBJS_NOPREFIX=$(SRC:%.f90=%.o)
+OBJS=$(addprefix $(OBJDIR), $(OBJS_NOPREFIX))
+
+# variable containing if git is used or not
+ifeq ($(wildcard ../../.git/HEAD ../../.git/index),)
+	OBJ_GIT =
+	USE_GIT =
+else
+	USE_GIT = $(shell which git)
+	OBJ_GIT = ../../.git/HEAD ../../.git/index
+endif
+
+
+$(OBJDIR)%.o:%.f90
+	@echo "--------------------------------"
+	@echo building $*.o with parameters \
+	FC=$(FC), FCFLAGS=$(FCFLAGS)
+	$(FC) -c $(FCFLAGS) $< -o $@
+
+$(MODDIR)%.mod:%.f90
+	@echo "--------------------------------"
+	@echo building $*.mod with parameters \
+	FC=$(FC), FCFLAGS=$(FCFLAGS)
+	$(FC) -c $(FCFLAGS) $<
+
+
+all: directories release
+release: $(EXE)
+debug: FCFLAGS = $(DEBUGFLAGS)
+debug: $(EXE)
+
+$(EXE) : gitversion $(OBJS)
+	@echo "--------------------------------"
+	@echo $(ROOT_PFOF)
+	@echo building $(EXE) with parameters \
+	FC=$(FC), LDFLAGS=$(LDFLAGS)
+	$(FC) -g -o $(EXE) $(OBJS) $(LDFLAGS)
+
+cleanall :
+	@echo "--------------------------------"
+	@echo cleaning compiled objects, modules and executable
+	rm -f $(OBJDIR)*.o $(MODDIR)*.mod *~ $(EXE) $(COMMONDIR)/gitversion.h
+
+clean:
+	@echo "--------------------------------"
+	@echo cleaning executable
+	rm -f  $(EXE)
+
+
+# create file for git version
+gitversion: $(OBJ_GIT)
+ifneq ($(USE_GIT), )
+	@cd $(COMMONDIR) ;\
+	echo "#define GITVERSION \"$(shell git rev-parse HEAD)\"" > gitversion.h
+else
+	@cd $(COMMONDIR) ;\
+	echo "#define GITVERSION \"no_version_found\"" > gitversion.h
+endif
+
+directories: ${MODDIR} ${OBJDIR}
+
+${OBJDIR}:
+	${MKDIR_P} $@
+
+${MODDIR}:
+	${MKDIR_P} $@
+
+
+.PHONY: gitversion clean cleanall directories all

amr2cell_hdf5/src/amr2cell_hdf5.f90

@@ -56,7 +56,7 @@ program amr2cell_hdf5
   type(ramses_output_t) :: ramses_output
 
   call mpi_init(mpierr)
-  call Hdf5_init()
+  call hdf5_init()
 
   call timer%Init()
 
@@ -66,10 +66,10 @@ program amr2cell_hdf5
 
   call amr2cell_hdf5_parameters%Read()
 
-  call amr2cell_hdf5_parameters%Check()
-
   call ramses_info%Read(trim(amr2cell_hdf5_parameters%data_dir)//trim(amr2cell_hdf5_parameters%info_filename), mpi_process)
 
+  call amr2cell_hdf5_parameters%Check(ramses_info)
+
   call ramses_output%Init(amr2cell_hdf5_parameters)
 
   call ramses_output%Read()
@@ -78,7 +78,7 @@ program amr2cell_hdf5
 
   call ramses_output%Deallocate()
 
-  call cells_array%Clean(amr2cell_hdf5_parameters%cube_level)
+  call cells_array%Clean(amr2cell_hdf5_parameters)
 
   call cubes_array%Init(amr2cell_hdf5_parameters, ramses_info)
 
@@ -90,14 +90,14 @@ program amr2cell_hdf5
 
   call cubes_array%Associate(cells_array, global_cells_per_cube)
 
-  call cubes_array%Write()
+  call cubes_array%Write(amr2cell_hdf5_parameters, ramses_output%hydro_descriptor)
 
   call timer%Finalize()
   if ( mpi_process%rank == 0 ) then
      call timer%Print(OUTPUT_UNIT)
   end if
 
-  call Hdf5_finalize()
+  call hdf5_finalize()
   call mpi_finalize(mpierr)
 
 end program amr2cell_hdf5

amr2cell_hdf5/src/amr2cell_hdf5_parameters_m.f90

@@ -51,20 +51,26 @@ module amr2cell_hdf5_parameters_m
      character(len=LEN_STRING) :: info_filename
      character(len=LEN_STRING) :: cube_filename
      integer :: cube_level
+     integer :: cell_level_min
+     integer :: cell_level_max
    contains
      procedure :: Check => Check_amr2cell_hdf5_parameters
      procedure :: Print => Print_amr2cell_hdf5_parameters
      procedure :: Read => Read_amr2cell_hdf5_parameters
+     procedure :: Read_hdf5 => Read_hdf5_amr2cell_hdf5_parameters
+     procedure :: Write_hdf5 => Write_hdf5_amr2cell_hdf5_parameters
   end type amr2cell_hdf5_parameters_t
 
 contains
   !----------------------------------------------------------------------------------------------------------------------------------
-  subroutine Check_amr2cell_hdf5_parameters(this)
+  subroutine Check_amr2cell_hdf5_parameters(this, ramses_info)
 
     use char_utils_m, only : Int_to_char5
     use constants_m, only : FILENAME_LEN
+    use ramses_info_m, only : ramses_info_t
 
-    class(amr2cell_hdf5_parameters_t), intent(in) :: this
+    class(amr2cell_hdf5_parameters_t), intent(inout) :: this
+    type(ramses_info_t), intent(in) :: ramses_info
 
     character(len=FILENAME_LEN) :: filename
     integer :: ncubes
@@ -85,6 +91,17 @@ contains
        call Algorithmic_error('Check_amr2cell_hdf5_parameters', 'Number of processes used should be <= number of cubes written.', 1000, mpi_process%rank)
     end if
 
+    ! check cell_level_min and cell_level_max
+    if( this%cell_level_min == 0 .or. this%cell_level_min < ramses_info%levelmin ) then
+      this%cell_level_min = ramses_info%levelmin
+      write(OUTPUT_UNIT,'(a,i0)') 'cell_level_min was undefined or < Ramses level_min ; its value has been adjusted to Ramses level_min = ', ramses_info%levelmin
+    end if
+
+    if( this%cell_level_max == 0 .or. this%cell_level_max > ramses_info%levelmax ) then
+      this%cell_level_max = ramses_info%levelmax
+      write(OUTPUT_UNIT,'(a,i0)') 'cell_level_max was undefined or > Ramses level_max ; its value has been adjusted to Ramses level_max = ', ramses_info%levelmax
+    end if
+
     ! test the existance of the files
     inquire(directory=this%data_dir, exist=data_dir_exists )
     if(.not. data_dir_exists) then
@@ -136,7 +153,7 @@ contains
     integer(kind=MPI_ADDRESS_KIND), allocatable, dimension(:) :: disp
     integer(kind=MPI_ADDRESS_KIND) :: base
     integer(kind=4) :: mpierr
-    integer(kind=4), parameter :: TYPESIZE = 7
+    integer(kind=4), parameter :: TYPESIZE = 9
     integer(kind=4) :: alloc_stat
     character(len=ERR_MSG_LEN) :: error_message
 
@@ -165,13 +182,15 @@ contains
     call mpi_get_address(parameters%info_filename, disp(5), mpierr)
     call mpi_get_address(parameters%cube_filename, disp(6), mpierr)
     call mpi_get_address(parameters%cube_level, disp(7), mpierr)
+    call mpi_get_address(parameters%cell_level_min, disp(8), mpierr)
+    call mpi_get_address(parameters%cell_level_max, disp(9), mpierr)
 
     base = disp(1)
     disp(:) = disp(:) - base
     blocklen(1:6) = LEN_STRING
-    blocklen(7:7) = 1
+    blocklen(7:9) = 1
     type(1:6) = MPI_CHARACTER
-    type(7:7) = MPI_INTEGER
+    type(7:9) = MPI_INTEGER
 
     call mpi_type_create_struct(typesize, blocklen, disp, type, mpi_type_a2cp, mpierr)
     call mpi_type_commit(mpi_type_a2cp, mpierr)
@@ -205,6 +224,8 @@ contains
     write(unit_number,'(a)') 'Output parameters:'    
     write(unit_number,'(a,a)') ' - cube filename: ', trim(this%cube_filename)
     write(unit_number,'(a,i0)') ' - cube level: ', this%cube_level
+    write(unit_number,'(a,i0)') ' - minimum cell level: ', this%cell_level_min
+    write(unit_number,'(a,i0)') ' - maximum cell level: ', this%cell_level_max
     write(unit_number,'(a,a)') '---------------------------------------------------------------------------------'
 
 #ifdef DEBUG
@@ -221,6 +242,8 @@ contains
     character(len=LEN_STRING) :: data_directory, amr_filename, hydro_filename, info_filename, &
          hydro_descriptor_filename
     character(len=LEN_STRING) :: cube_filename
+    integer :: cell_level_max = 0
+    integer :: cell_level_min = 0
     integer :: cube_level
     character(ERR_MSG_LEN) :: error_message
     integer :: ioerr
@@ -233,7 +256,7 @@ contains
 
     namelist / input_parameters / data_directory, amr_filename, hydro_filename, info_filename, &
          hydro_descriptor_filename
-    namelist / output_parameters / cube_level, cube_filename
+    namelist / output_parameters / cube_level, cell_level_min, cell_level_max, cube_filename
 
 #ifdef DEBUG
     write(ERROR_UNIT,'(a,i0)') 'Read_amr2cell_hdf5_parameters begins on process', mpi_process%rank
@@ -275,7 +298,9 @@ contains
        this%hydro_descriptor_filename = hydro_descriptor_filename
        this%info_filename = info_filename
        this%cube_filename = cube_filename
-       this%cube_level = cube_level       
+       this%cube_level = cube_level
+       this%cell_level_min = cell_level_min
+       this%cell_level_max = cell_level_max
     end if
     call timer%Inc_inp()
 
@@ -287,11 +312,114 @@ contains
     call timer%Inc_out()
 
 #ifdef DEBUG
-    ! check parameters on process 1:
     write(ERROR_UNIT,'(a,i0)') 'Read_amr2cell_hdf5_parameter ends on process', mpi_process%rank
 #endif
 
   end subroutine Read_amr2cell_hdf5_parameters
 
+!----------------------------------------------------------------------------------------------------------------------------------
+  subroutine Read_hdf5_amr2cell_hdf5_parameters(this, hdf5_id, mpi_process)
+
+    use fortran_hdf5_constants_m, only : H5_STR_LEN
+    use fortran_hdf5_manage_groups_m, only : Hdf5_close_group, Hdf5_open_group
+    use fortran_hdf5_read_attribute_m, only : Hdf5_read_attr
+    use hdf5, only : HID_T
+    use mpi_process_m, only : mpi_process_t
+
+    integer(HID_T), intent(in) :: hdf5_id
+    type(mpi_process_t), intent(in) :: mpi_process
+    class(amr2cell_hdf5_parameters_t), intent(out) :: this
+
+    character(len=H5_STR_LEN) :: groupname
+    character(len=H5_STR_LEN) :: aname
+    integer(kind=HID_T) :: gr_id
+   
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_amr2cell_hdf5_parameters begins on process', mpi_process%rank
+#endif
+
+    groupname = 'amr2cell_hdf5_parameters'
+    call Hdf5_open_group(hdf5_id,groupname,gr_id)
+
+    aname = 'data_dir'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%data_dir)
+    aname = 'amr_filename'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%amr_filename)
+    aname = 'hydro_filename'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%hydro_filename)
+    aname = 'hydro_descriptor_filename'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%hydro_descriptor_filename)
+    aname = 'info_filename'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%info_filename)
+    aname = 'cube_filename'
+    call Hdf5_read_attr(gr_id,aname,LEN_STRING,this%cube_filename)
+    aname = 'cube_level'
+    call Hdf5_read_attr(gr_id,aname,this%cube_level)
+    aname = 'cell_level_min'
+    call Hdf5_read_attr(gr_id,aname,this%cell_level_min)
+    aname = 'cell_level_max'
+    call Hdf5_read_attr(gr_id,aname,this%cell_level_max)
+
+    call Hdf5_close_group(gr_id)
+
+
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_amr2cell_hdf5_parameters ends on process', mpi_process%rank
+#endif
+
+  end subroutine Read_hdf5_amr2cell_hdf5_parameters
+
+  !----------------------------------------------------------------------------------------------------------------------------------
+  subroutine Write_hdf5_amr2cell_hdf5_parameters(this, hdf5_id, mpi_process)
+
+    use fortran_hdf5_constants_m, only : H5_STR_LEN
+    use fortran_hdf5_manage_groups_m, only : Hdf5_close_group, Hdf5_create_group
+    use fortran_hdf5_write_attribute_m, only : Hdf5_write_attr
+    use hdf5, only : HID_T
+    use mpi_process_m, only : mpi_process_t
+
+    integer(HID_T), intent(in) :: hdf5_id
+    type(mpi_process_t), intent(in) :: mpi_process
+    class(amr2cell_hdf5_parameters_t), intent(in) :: this
+
+    character(len=H5_STR_LEN) :: groupname
+    character(len=H5_STR_LEN) :: aname
+    integer(kind=HID_T) :: gr_id
+   
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Write_hdf5_amr2cell_hdf5_parameters begins on process', mpi_process%rank
+#endif
+
+    groupname = 'amr2cell_hdf5_parameters'
+    call Hdf5_create_group(hdf5_id,groupname,gr_id)
+
+    aname = 'data_dir'
+    call Hdf5_write_attr(gr_id,aname,this%data_dir)
+    aname = 'amr_filename'
+    call Hdf5_write_attr(gr_id,aname,this%amr_filename)
+    aname = 'hydro_filename'
+    call Hdf5_write_attr(gr_id,aname,this%hydro_filename)
+    aname = 'hydro_descriptor_filename'
+    call Hdf5_write_attr(gr_id,aname,this%hydro_descriptor_filename)
+    aname = 'info_filename'
+    call Hdf5_write_attr(gr_id,aname,this%info_filename)
+    aname = 'cube_filename'
+    call Hdf5_write_attr(gr_id,aname,this%cube_filename)
+    aname = 'cube_level'
+    call Hdf5_write_attr(gr_id,aname,this%cube_level)
+    aname = 'cell_level_min'
+    call Hdf5_write_attr(gr_id,aname,this%cell_level_min)
+    aname = 'cell_level_max'
+    call Hdf5_write_attr(gr_id,aname,this%cell_level_max)
+
+    call Hdf5_close_group(gr_id)
+
+
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Write_hdf5_amr2cell_hdf5_parameters ends on process', mpi_process%rank
+#endif
+
+  end subroutine Write_hdf5_amr2cell_hdf5_parameters
+
 end module amr2cell_hdf5_parameters_m
 

amr2cell_hdf5/src/cell_m.f90

@@ -40,8 +40,8 @@ module cell_m
   public :: cell_t, Init_cell_mpi_type
 
   type cell_t(n_non_thermal_pressures, n_passive_scalars)
-     integer, len :: n_non_thermal_pressures
-     integer, len :: n_passive_scalars
+     integer(kind=4), len :: n_non_thermal_pressures
+     integer(kind=4), len :: n_passive_scalars
      integer(kind=4) :: level
      integer(kind=4) :: process_id
      integer(kind=4) :: cube_id
@@ -56,7 +56,8 @@ module cell_m
      procedure :: Is_eq
      procedure :: Is_sup
      procedure :: Print => Print_cell
-     procedure :: Write => Write_hdf5_cell
+     procedure :: Read_hdf5 => Read_hdf5_cell
+     procedure :: Write_hdf5 => Write_hdf5_cell
   end type cell_t
 
 contains
@@ -250,6 +251,56 @@ contains
 
   end subroutine Print_cell
 
+  !------------------------------------------------------------------------------------------------------------------------------------
+  subroutine Read_hdf5_cell(this, hdf5_id)
+
+   use char_utils_m, only : Int_to_char12, Int_to_char5
+   use fortran_hdf5_constants_m
+   use fortran_hdf5_manage_groups_m
+   use fortran_hdf5_read_attribute_m
+   use fortran_hdf5_read_data_m
+   use hdf5
+
+   class(cell_t(*,*)),intent(inout) :: this
+   integer(kind=hid_t),intent(in) :: hdf5_id
+
+   character(len=H5_STR_LEN) :: attr_name
+   character(len=H5_STR_LEN) :: groupname
+   integer(kind=hid_t) :: cell_gr_id
+   integer :: idim
+
+#ifdef DEBUG
+   write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_cell begins on process', mpi_process%rank
+#endif
+   call timer%Set_ref()
+
+   attr_name = 'identity'
+   call Hdf5_Read_data(hdf5_id, attr_name, this%identity)
+   attr_name = 'level'
+   call Hdf5_Read_attr(hdf5_id, attr_name, this%level)
+   attr_name = 'position'
+   call Hdf5_Read_attr(hdf5_id, attr_name, 3, this%position)
+   attr_name = 'velocity field'
+   call Hdf5_Read_attr(hdf5_id, attr_name, 3, this%velocity_field)
+   attr_name = 'density'
+   call Hdf5_Read_attr(hdf5_id, attr_name, this%density)
+   attr_name = 'thermal pressure'
+   call Hdf5_Read_attr(hdf5_id, attr_name, this%thermal_pressure)
+   if ( this%n_non_thermal_pressures /= 0 ) then
+      attr_name = 'non thermal pressures'
+!      call Hdf5_Read_attr(hdf5_id, attr_name, this%n_non_thermal_pressures, this%non_thermal_pressures)
+   end if
+   if( this%n_passive_scalars /= 0 ) then 
+      attr_name = 'passive scalars'
+      call Hdf5_Read_attr(hdf5_id, attr_name, this%n_passive_scalars, this%passive_scalars)
+   end if
+
+   call timer%Inc_inp()
+#ifdef DEBUG
+   write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_cell ends on process', mpi_process%rank
+#endif
+
+ end subroutine Read_hdf5_cell
 
   !------------------------------------------------------------------------------------------------------------------------------------
   subroutine Write_hdf5_cell(this, hdf5_id)

amr2cell_hdf5/src/cells_array_m.f90

@@ -49,13 +49,15 @@ module cells_array_m
      procedure Deallocate => Deallocate_cells_array
      procedure Distribute => Distribute_cells_array
      procedure Clean => Clean_unused_cells 
+     procedure Read_hdf5 => Read_hdf5_cells_array
      procedure, private :: Remove_duplicate_cells
      procedure, private :: Remove_large_cells
+     procedure, private :: Remove_small_cells
      procedure Sort => Sort_cells_array
      procedure, private :: Sort_heapify_cells
      procedure, private :: Sort_sink_cells
      procedure, private :: Sort_swap_cells
-     procedure :: Write => Write_hdf5_cells_array
+     procedure :: Write_hdf5 => Write_hdf5_cells_array
   end type cells_array_t
 
 contains
@@ -89,14 +91,17 @@ contains
   end subroutine Allocate_cells_array
 
   !------------------------------------------------------------------------------------------------------------------------------------
-  subroutine Clean_unused_cells(this, cube_level)
-
+  subroutine Clean_unused_cells(this, amr2cell_hdf5_parameters)
+    
+    use amr2cell_hdf5_parameters_m, only : amr2cell_hdf5_parameters_t
+    
+    type(amr2cell_hdf5_parameters_t), intent(in) :: amr2cell_hdf5_parameters
     class(cells_array_t),intent(inout) :: this
-    integer, intent(in) :: cube_level
 
     call this%Sort()
     call this%Remove_duplicate_cells()
-    call this%Remove_large_cells(cube_level)
+    call this%Remove_large_cells(amr2cell_hdf5_parameters%cell_level_min)
+    call this%Remove_small_cells(amr2cell_hdf5_parameters%cell_level_max)
 
   end subroutine Clean_unused_cells
 
@@ -322,7 +327,51 @@ contains
   end subroutine Distribute_cells_array
 
   !------------------------------------------------------------------------------------------------------------------------------------
+  subroutine Read_hdf5_cells_array(this, hdf5_id, ramses_hydro_descriptor)
+
+    use char_utils_m, only : int_to_char12
+    use fortran_hdf5_constants_m
+    use fortran_hdf5_manage_groups_m
+    use fortran_hdf5_read_attribute_m
+    use fortran_hdf5_read_data_m
+    use hdf5
+    use ramses_hydro_descriptor_m, only : ramses_hydro_descriptor_t
+
+    class(cells_array_t),intent(out) :: this
+    integer(kind=hid_t),intent(in) :: hdf5_id
+    type(ramses_hydro_descriptor_t), intent(in) :: ramses_hydro_descriptor
+
+    character(len=H5_STR_LEN) :: attr_name
+    character(len=H5_STR_LEN) :: groupname
+    integer(kind=hid_t) :: cell_gr_id
+    integer :: icell
+    integer(IDKIND) :: ncells
+    integer(kind=4) :: nps
+    integer(kind=4) :: nthp
+
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_cells_array begins on process', mpi_process%rank
+#endif
+    ! read number of cells in this array
+    attr_name = 'ncells'
+    call Hdf5_read_data(hdf5_id,attr_name, ncells)
 
+    call this%Allocate(ramses_hydro_descriptor%n_non_thermal_pressures,ramses_hydro_descriptor%n_passive_scalars,ncells)
+
+    do icell = 1, this%ncells
+       groupname = 'cell_'//int_to_char12(icell)
+       call Hdf5_open_group(hdf5_id, groupname, cell_gr_id)
+       call this%elts(icell)%Read_hdf5(cell_gr_id)
+       call Hdf5_close_group(cell_gr_id)
+    end do
+
+#ifdef DEBUG
+    write(ERROR_UNIT,'(a,i0)') 'Read_hdf5_cells_array ends on process', mpi_process%rank
+#endif
+
+  end subroutine Read_hdf5_cells_array
+
+  !------------------------------------------------------------------------------------------------------------------------------------
   subroutine Remove_duplicate_cells(this)
 
     class(cells_array_t), intent(inout) :: this
@@ -383,10 +432,10 @@ contains
   end subroutine Remove_duplicate_cells