!*****************
! File Description:  
!   Library to manage HDF5 file (working with SparkCode)
! Creator: Hoang-Trong Minh Tuan 
! 
! Rev: 0.1 
! Update: 
!    x.x: 
!*******************/

MODULE hdf5_utility
  use precisioninfo
  USE hdf5
  IMPLICIT NONE
  ! 0 or 1
#define _ENABLE_CACHE_ 0

#define _ENABLE_SZIP_ 0


  INTEGER, PARAMETER :: MAX_NUM_HANDLE = 100, MAX_NUM_GROUPS = 4
  ! keep the handle to the HDF5 objects
  TYPE HDFFILE !regular data
     INTEGER(HID_T) :: file_ID
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_dspace ! (data) space
     INTEGER ::  num_dspace
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_dset ! dataset
     INTEGER :: num_dset
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_plist !property list
     INTEGER :: num_plist
  END TYPE HDFFILE
  
  TYPE HDFFILE2 !to deal with spatial data
     INTEGER(HID_T) :: file_ID
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_dspace ! (data) space
     INTEGER ::  num_dspace
     !first = vlayer
     !second= hlayer
     INTEGER(HID_T), DIMENSION(2,MAX_NUM_HANDLE,MAX_NUM_GROUPS) :: arr_dset ! dataset (1:2,dataset_id, group_id)
     INTEGER, DIMENSION(MAX_NUM_GROUPS) :: num_dset !number of dset in each group
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_plist !property list
     INTEGER :: num_plist
     INTEGER(HID_T), DIMENSION(MAX_NUM_GROUPS) :: arr_group ! group
     INTEGER :: num_group
     REAL(KIND=dp), DIMENSION(MAX_NUM_HANDLE*MAX_NUM_GROUPS) :: time_data
  END TYPE HDFFILE2

  TYPE HDFFILE3 !spatial whole-cell data
     INTEGER(HID_T) :: file_ID
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_dspace ! (data) space
     INTEGER ::  num_dspace
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_dset ! dataset
     INTEGER :: num_dset
     INTEGER(HID_T), DIMENSION(MAX_NUM_HANDLE) :: arr_plist !property list
     INTEGER :: num_plist
     REAL(KIND=dp), DIMENSION(MAX_NUM_HANDLE) :: time_data
  END TYPE HDFFILE3

  ! file size of each family member
  INTEGER(HSIZE_T), PARAMETER :: SP_FAMILY_SIZE = 100*1024*1024 !(100 MB)

CONTAINS

  !! Init HDF5
  FUNCTION hdf_init() RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    !local
    INTEGER :: hdf_error

    CALL h5open_f(hdf_error)
    res = hdf_error
  END FUNCTION  hdf_init

  !! Close HDF5
  FUNCTION hdf_close() RESULT(res)
    IMPLICIT NONE
    INTEGER :: res
    !local
    INTEGER :: hdf_error

    CALL h5close_f(hdf_error)
    res = hdf_error
  END FUNCTION  hdf_close


  !! Open file and setting up its structure
  ! No compression is supported
  !!NOTE: Currently works with CHUNKS_PER_FILE = 1 only
  ! Will be resolved when hyperslab is added to hdf_spark_write_simple()
  !!!
  SUBROUTINE hdf_open_simple(filename, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, numfields)
    ! NSFU = number of SFU 
    ! CHUNK_SIZE = the number of record to keep before flushing out (A.K.A. a chunk)
    ! CHUNKS_PER_FILE = number of chunk in a file
    ! [OUT]file_struct = structure that hold all information about the HDF file
    ! [OUT]numfields = number of data per record 
    IMPLICIT NONE
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: NSFU
    INTEGER, INTENT(IN) :: CHUNK_SIZE, CHUNKS_PER_FILE
    TYPE(HDFFILE), INTENT(INOUT) :: file_struct
    INTEGER, INTENT(IN) :: numfields
    ! constant
    INTEGER, PARAMETER :: RANK = 2 ! 2D-dataspace
    ! local
    INTEGER (HID_T) :: file_ID, space_id, dset_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER :: hdf_error
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0

    !create the HDF5 file

    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
    !setting up dataspace + dataset
    ! (1) time
    ! (2) V_m
    ! (3) Ca_myo
    ! (4) Ca_nsr
    ! (5) isfu(1:NSFU)
    ! (6) Ca_ds(1:NSFU)
    ! (7) Ca_jsr(1:NSFU)
    ! (8) pR_nj(1:mR)  --- leak
    !     (pR_nj(+)pL_nj) mR*mL  --- ECC
    ! NOTE: grows in columns, not in rows
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    current_dims = (/ numfields, CHUNK_SIZE*CHUNKS_PER_FILE /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  1
    file_struct%arr_dspace(file_struct%num_dspace) = space_id

    CALL h5dcreate_f(file_ID, "/LEAK", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    file_struct%num_dset = 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
  END SUBROUTINE hdf_open_simple

  SUBROUTINE hdf_open_simple_leak(filename, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, numfields)
    ! NSFU = number of SFU 
    ! CHUNK_SIZE = the number of record to keep before flushing out (A.K.A. a chunk)
    ! CHUNKS_PER_FILE = number of chunk in a file
    ! [OUT]file_struct = structure that hold all information about the HDF file
    ! [OUT]numfields = number of data per record 
    IMPLICIT NONE
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: NSFU
    INTEGER, INTENT(IN) :: CHUNK_SIZE, CHUNKS_PER_FILE
    TYPE(HDFFILE), INTENT(INOUT) :: file_struct
    INTEGER, INTENT(IN) :: numfields
    ! constant
    INTEGER, PARAMETER :: RANK = 2 ! 2D-dataspace
    ! local
    INTEGER (HID_T) :: file_ID, space_id, dset_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER :: hdf_error
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0

    !create the HDF5 file

    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
    !setting up dataspace + dataset
    ! (1) time
    ! (2) V_m
    ! (3) Ca_myo
    ! (4) Ca_nsr
    ! (5) isfu(1:NSFU)
    ! (6) Ca_ds(1:NSFU)
    ! (7) Ca_jsr(1:NSFU)
    ! (8) pR_nj(1:mR)  --- leak
    !     (pR_nj(+)pL_nj) mR*mL  --- ECC
    ! NOTE: grows in columns, not in rows
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    current_dims = (/ numfields, CHUNK_SIZE*CHUNKS_PER_FILE /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  1
    file_struct%arr_dspace(file_struct%num_dspace) = space_id

    CALL h5dcreate_f(file_ID, "/LEAK", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    file_struct%num_dset = 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
  END SUBROUTINE hdf_open_simple_leak

  SUBROUTINE hdf_open_simple_ecc(filename, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, numfields)
    ! NSFU = number of SFU 
    ! CHUNK_SIZE = the number of record to keep before flushing out (A.K.A. a chunk)
    ! CHUNKS_PER_FILE = number of chunk in a file
    ! [OUT]file_struct = structure that hold all information about the HDF file
    ! [OUT]numfields = number of data per record 
    IMPLICIT NONE
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: NSFU
    INTEGER, INTENT(IN) :: CHUNK_SIZE, CHUNKS_PER_FILE
    TYPE(HDFFILE), INTENT(INOUT) :: file_struct
    INTEGER, INTENT(IN) :: numfields
    ! constant
    INTEGER, PARAMETER :: RANK = 2 ! 2D-dataspace
    ! local
    INTEGER (HID_T) :: file_ID, space_id, dset_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER :: hdf_error
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0

    !create the HDF5 file

    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
    !setting up dataspace + dataset
    ! (1) time
    ! (2) V_m
    ! (3) Ca_myo
    ! (4) Ca_nsr
    ! (5) isfu(1:NSFU)
    ! (6) Ca_ds(1:NSFU)
    ! (7) Ca_jsr(1:NSFU)
    ! (8) pR_nj(1:mR)  --- leak
    !     (pR_nj(+)pL_nj) mR*mL  --- ECC
    ! NOTE: grows in columns, not in rows
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    current_dims = (/ numfields, CHUNK_SIZE*CHUNKS_PER_FILE /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  1
    file_struct%arr_dspace(file_struct%num_dspace) = space_id

    CALL h5dcreate_f(file_ID, "/ECC", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    file_struct%num_dset = 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
  END SUBROUTINE hdf_open_simple_ecc

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! Open file and setting up its structure
  !! use file families (not supported by MatLab)
  SUBROUTINE hdf_spark_open_filefamily(NSFU, BUFFER_SIZE, file_struct, prealloc)
    ! NSFU = number of SFU 
    ! BUFFER_SIZE = the number of record to keep before flushing out (whole package is called a chunk)
    ! file_struct = structure that hold all information about the HDF file
    ! prealloc = amount of data (in multiple of BUFFERSIZE) to be allocated in the beginning
    ! (default: CONST_PREALLOC chunks)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NSFU
    INTEGER, INTENT(IN) :: BUFFER_SIZE
    TYPE(HDFFILE), INTENT(INOUT) :: file_struct
    INTEGER, OPTIONAL :: prealloc
    ! constant
    INTEGER, PARAMETER :: RANK = 2 ! 2D-dataspace
    INTEGER, PARAMETER :: CONST_PREALLOC = 100
    ! local
    CHARACTER(LEN=30) :: filename
    INTEGER (HID_T) :: file_ID, space_id, dset_id,  space_id_2
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims, max_dims, chunk_dims
    INTEGER :: numfields, hdf_error
    INTEGER :: allocfirst
    INTEGER(HID_T) :: faccess_pl_id, dcreation_pl_id, daccess_pl_id, &
         dcreation_pl_id_2
    INTEGER :: iN, pixels_per_block, option_mask
    LOGICAL  :: avail
    INTEGER :: filter_info, filter_info_both
    INTEGER :: ch_row, ch_col
    INTEGER(SIZE_T) :: rdcc_nbytes, rdcc_nelmts
    INTEGER :: mdc_nelmts
    REAL :: rdcc_w0


    IF (PRESENT(prealloc)) THEN
       allocfirst = prealloc
    ELSE
       allocfirst = CONST_PREALLOC
    END IF
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0

#if _ENABLE_SZIP_ == 1    
    !check SZIP
    CALL h5zfilter_avail_f(H5Z_FILTER_SZIP_F, avail, hdf_error)
    IF (.NOT. avail) THEN
       WRITE(*, *) "SZIP filter not available"
       STOP
    END IF
    CALL h5zget_filter_info_f(H5Z_FILTER_SZIP_F, filter_info, hdf_error)
    filter_info_both = IOR(H5Z_FILTER_ENCODE_ENABLED_F, H5Z_FILTER_DECODE_ENABLED_F)
    IF (filter_info .NE. filter_info_both) THEN
       WRITE (*,*) "SZIP filter not available for encoding and decoding"
       STOP
    ENDIF
#endif    
    !create the HDF5 file (file family)
    filename = "sprk_%d_.h5"
    !+ access-property list
    CALL H5Pcreate_f(H5P_FILE_ACCESS_F, faccess_pl_id, hdf_error) !default file-access property
    CALL H5Pset_fapl_family_f(faccess_pl_id, SP_FAMILY_SIZE, &
         H5P_DEFAULT_F, hdf_error)
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = faccess_pl_id
    !+ set chunk cache
#if _ENABLE_CACHE_ == 1    
    mdc_nelmts = 0 ! not important
    rdcc_nelmts = 521 ! # objects in chunk cache
    rdcc_nbytes = 20*1024*1024  ! cache size (in 20MB)
    rdcc_w0 = 0.5
    CALL h5pset_cache_f(faccess_pl_id, mdc_nelmts, rdcc_nelmts, rdcc_nbytes, rdcc_w0, hdf_error)
#endif    
    CALL h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F, faccess_pl_id)
    file_struct%file_ID = file_ID
    !setting up dataspace + dataset
    ! (1) time
    ! (2) V_m
    ! (3) Ca_myo
    ! (4) Ca_nsr
    ! (5) isfu(1:NSFU)
    ! (6) Ca_ds(1:NSFU)
    ! (7) Ca_jsr(1:NSFU)
    ! NOTE: grows in columns, not in rows
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !(1) 
    numfields = 1
    current_dims = (/ numfields, allocfirst*BUFFER_SIZE /) !size at
    !creation time
    max_dims = (/ numfields, H5S_UNLIMITED_F /)
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error, max_dims)
    file_struct%num_dspace = file_struct%num_dspace + 1
    file_struct%arr_dspace(file_struct%num_dspace) = space_id
    !+create Data-Access PL
    CALL h5pcreate_f(H5P_DATASET_CREATE_F, dcreation_pl_id, hdf_error)
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id

    ! set chunk 
    iN = 128
    pixels_per_block = 16  ! 8, 16 or 32
    ch_row = numfields
    ch_col = iN*pixels_per_block
    chunk_dims = (/ ch_row, ch_col /)
    CALL h5pset_chunk_f(dcreation_pl_id, RANK, chunk_dims, hdf_error)
    ! set SZIP
    option_mask=H5_SZIP_NN_OM_F;
    CALL H5Pset_szip_f (dcreation_pl_id, option_mask, pixels_per_block, hdf_error)    
    ! create dataset
    CALL h5dcreate_f(file_ID, "/time", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(2)
    CALL h5scopy_f(space_id, space_id_2, hdf_error)
    space_id = space_id_2
    CALL h5pcopy_f(dcreation_pl_id, dcreation_pl_id_2, hdf_error)
    dcreation_pl_id = dcreation_pl_id_2
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id
    ! create dataset
    CALL h5dcreate_f(file_ID, "/membrane_voltage", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(3)
    CALL h5scopy_f(space_id, space_id_2, hdf_error)
    space_id = space_id_2
    CALL h5pcopy_f(dcreation_pl_id, dcreation_pl_id_2, hdf_error)
    dcreation_pl_id = dcreation_pl_id_2
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id
    ! create dataset
    CALL h5dcreate_f(file_ID, "/calcium_myo", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(4)
    CALL h5scopy_f(space_id, space_id_2, hdf_error)
    space_id = space_id_2
    CALL h5pcopy_f(dcreation_pl_id, dcreation_pl_id_2, hdf_error)
    dcreation_pl_id = dcreation_pl_id_2
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id
    ! create dataset
    CALL h5dcreate_f(file_ID, "/calcium_nsr", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(5)
    numfields = NSFU
    current_dims = (/ numfields, allocfirst*BUFFER_SIZE /) !size at
    !creation time
    max_dims = (/ numfields, H5S_UNLIMITED_F /)
    !create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error, max_dims)
    file_struct%num_dspace = file_struct%num_dspace + 1
    file_struct%arr_dspace(file_struct%num_dspace) = space_id

    CALL h5pcreate_f(H5P_DATASET_CREATE_F, dcreation_pl_id, hdf_error)
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id

    ! set chunk 
    iN = 128
    pixels_per_block = 8  ! 8, 16 or 32
    ch_row = 500 ! read every 500 release site
    ch_col = iN*pixels_per_block
    chunk_dims = (/ ch_row, ch_col /)
    CALL h5pset_chunk_f(dcreation_pl_id, RANK, chunk_dims, hdf_error)
    ! set SZIP
    option_mask=H5_SZIP_NN_OM_F;
    CALL H5Pset_szip_f (dcreation_pl_id, option_mask, pixels_per_block, hdf_error)

    ! create dataset
    CALL h5dcreate_f(file_ID, "/isfu_state", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(6)
    CALL h5scopy_f(space_id, space_id_2, hdf_error)
    space_id = space_id_2
    CALL h5pcopy_f(dcreation_pl_id, dcreation_pl_id_2, hdf_error)
    dcreation_pl_id = dcreation_pl_id_2
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id
    ! create dataset
    CALL h5dcreate_f(file_ID, "/calcium_ds", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

    !(7)
    CALL h5scopy_f(space_id, space_id_2, hdf_error)
    space_id = space_id_2
    CALL h5pcopy_f(dcreation_pl_id, dcreation_pl_id_2, hdf_error)
    dcreation_pl_id = dcreation_pl_id_2
    file_struct%num_plist = file_struct%num_plist + 1
    file_struct%arr_plist(file_struct%num_plist) = dcreation_pl_id
    ! create dataset
    CALL h5dcreate_f(file_ID, "/calcium_jsr", H5T_NATIVE_REAL, space_id, &
         dset_id, hdf_error, dcreation_pl_id)
    file_struct%num_dset = file_struct%num_dset + 1
    file_struct%arr_dset(file_struct%num_dset) = dset_id
    ! optional
    CALL h5dextend_f(dset_id, current_dims, hdf_error)

  END SUBROUTINE hdf_spark_open_filefamily


  !! Close a file
  SUBROUTINE hdf_spark_close(file_struct)
    IMPLICIT NONE
    TYPE(HDFFILE) :: file_struct
    INTEGER :: hdf_error, ii

    DO ii = 1, file_struct%num_plist
       CALL h5pclose_f(file_struct%arr_plist(ii), hdf_error)
    END DO
    DO ii = 1, file_struct%num_dset
       CALL h5dclose_f(file_struct%arr_dset(ii), hdf_error)
    END DO
    DO ii = 1, file_struct%num_dspace
       CALL h5sclose_f(file_struct%arr_dspace(ii), hdf_error)
    END DO
    CALL h5fclose_f(file_struct%file_ID, hdf_error)
  END SUBROUTINE hdf_spark_close

  SUBROUTINE hdf_spark_write(file_struct, data_track, size, NSFU)
    TYPE(HDFFILE) :: file_struct
    REAL(KIND=real8), DIMENSION(:,:), ALLOCATABLE :: data_track
    INTEGER :: size, NSFU
    !local
    INTEGER, PARAMETER :: rank = 2
    INTEGER(HSIZE_T), DIMENSION(rank) :: dims
    INTEGER :: hdf_error

    dims = (/ 1, size /)
    ! time
    CALL h5dwrite_f(file_struct%arr_dset(1), &
         H5T_NATIVE_REAL, real(data_track(1,:)), &
         dims, hdf_error)
    ! Vm
    PRINT *, "Test 1"
    CALL h5dwrite_f(file_struct%arr_dset(2), &
         H5T_NATIVE_REAL, real(data_track(2,:)), &
         dims, hdf_error)
    ! Ca_myo
    PRINT *, "Test 2"
    CALL h5dwrite_f(file_struct%arr_dset(3), &
         H5T_NATIVE_REAL, real(data_track(3,:)), &
         dims, hdf_error)
    ! Ca_nsr
    PRINT *, "Test 3"
    CALL h5dwrite_f(file_struct%arr_dset(4), &
         H5T_NATIVE_REAL, real(data_track(4,:)), &
         dims, hdf_error)
    ! isfu
    dims = (/ NSFU, size /)
    PRINT *, "Test 4"
    CALL h5dwrite_f(file_struct%arr_dset(5), &
         H5T_NATIVE_REAL, real(data_track(5:NSFU+4,:)), &
         dims, hdf_error)
    ! Ca_ds
    PRINT *, "Test 5"
    CALL h5dwrite_f(file_struct%arr_dset(6), &
         H5T_NATIVE_REAL, real(data_track(NSFU+5:2*NSFU+4,:)), &
         dims, hdf_error)
    ! Ca_jsr
    PRINT *, "Test 6"
    CALL h5dwrite_f(file_struct%arr_dset(7), &
         H5T_NATIVE_REAL, real(data_track(2*NSFU+5:3*NSFU+4,:)), &
         dims, hdf_error)

  END SUBROUTINE hdf_spark_write

  SUBROUTINE hdf_write_simple(file_struct, data_buffer, numrecords, numfields)
    IMPLICIT NONE
    TYPE(HDFFILE) :: file_struct
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: data_buffer
    INTEGER :: numrecords, numfields
    !local
    INTEGER, PARAMETER :: rank = 2
    INTEGER(HSIZE_T), DIMENSION(rank) :: dims
    INTEGER :: hdf_error

    dims = (/ numfields, numrecords /)
    CALL h5dwrite_f(file_struct%arr_dset(1), &
         H5T_NATIVE_REAL, REAL(data_buffer(1:numfields,1:numrecords)), &
         dims, hdf_error)
  END SUBROUTINE hdf_write_simple


  
  !! Open file and setting up its structure
  !! No compression is supported
  SUBROUTINE hdf_spatial_open(file_struct, filename, num_vlayers, num_hlayers, outerDIMX, outerDIMY, outerDIMZ)
    !create file, first group, 2 data-spaces (one to deal with vertical layers, one to deal with horizontal layers)
    IMPLICIT NONE
    TYPE(HDFFILE2), INTENT(INOUT) :: file_struct
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: num_vlayers, num_hlayers
    INTEGER, INTENT(IN) :: outerDIMX, outerDIMY, outerDIMZ
    ! constant
    INTEGER, PARAMETER :: RANK = 3 ! 3D-dataspace
    ! local
    INTEGER (HID_T) :: file_ID, space_id, group_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER :: hdf_error, group_index
    CHARACTER(100) :: group_name
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0
    file_struct%num_group = 0
    file_struct%time_data = 0.0d0
    
    !create the HDF5 file
    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F, H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
    !create the first group
    group_index = file_struct%num_group
    CALL get_new_hdf_group_name(group_name, group_index)
    CALL h5gcreate_f(file_ID, TRIM(group_name), group_id, hdf_error)
    file_struct%num_group = group_index
    file_struct%arr_group(file_struct%num_group) = group_id

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !setting up dataspace + dataset
    ! We need totally 2 data spaces
    !+create first data space (fixed size) - for vertical layers
    current_dims = (/ outerDIMX, outerDIMY, num_vlayers /) !size at creation time
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  1
    file_struct%arr_dspace(1) = space_id
    !+create second dataspace (fixed size) - for horizontal layers
    current_dims = (/ outerDIMY, outerDIMZ, num_hlayers /)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  2
    file_struct%arr_dspace(2) = space_id
  END SUBROUTINE hdf_spatial_open

  SUBROUTINE get_new_hdf_group_name(group_name, group_index)
    !get new group name by automatically increase group_index
    IMPLICIT NONE
    CHARACTER*(*) :: group_name
    INTEGER, INTENT(INOUT) :: group_index
    
    group_index = group_index + 1
    WRITE(group_name, "('/group-', i3.3)") group_index
  END SUBROUTINE get_new_hdf_group_name
  
  !SUBROUTINE  hdf_spatial_write(file_struct, num_dataset, arrObj, num_vlayers, idx_vlayers, num_hlayers, idx_hlayers, &
  SUBROUTINE  hdf_spatial_write(file_struct, arrObj, num_vlayers, idx_vlayers, num_hlayers, idx_hlayers, &
       dimX_arr, dimY_arr, dimZ_arr, time, is_Full_File)
    IMPLICIT NONE
    TYPE(HDFFILE2), INTENT(INOUT) :: file_struct
    
    REAL(KIND=dp), DIMENSION(:,:,:), INTENT(IN) :: arrObj
    INTEGER, INTENT(IN):: num_vlayers, num_hlayers
    INTEGER, DIMENSION(:), INTENT(IN) :: idx_vlayers, idx_hlayers
    INTEGER, INTENT(IN):: dimX_arr, dimY_arr, dimZ_arr
    REAL(KIND=dp), INTENT(IN) :: time
    LOGICAL, INTENT(INOUT) :: is_full_file
    !local
    CHARACTER(100):: dataset_name, group_name
    INTEGER, PARAMETER :: rank = 3
    INTEGER(HSIZE_T), DIMENSION(rank) :: dims
    INTEGER(HID_T) :: dset_id, group_id
    INTEGER :: hdf_error, group_index
    INTEGER :: num_time_point, num_dataset

    group_index = file_struct%num_group
    IF (file_struct%num_dset(group_index) .EQ. MAX_NUM_HANDLE) THEN
       !create_new_group
       CALL get_new_hdf_group_name(group_name, group_index)
       CALL h5gcreate_f(file_struct%file_ID, TRIM(group_name), group_id, hdf_error)
       file_struct%num_group = group_index
       file_struct%arr_group(group_index) = group_id
       file_struct%num_dset(group_index) = 0
    END IF    
    num_dataset = file_struct%num_dset(group_index)
    IF (num_dataset .LT. MAX_NUM_HANDLE) THEN
       num_dataset = num_dataset + 1
       WRITE(dataset_name, "('vdataset-', i3.3)") num_dataset
       !put to the current group
       CALL h5dcreate_f(file_struct%arr_group(file_struct%num_group), TRIM(dataset_name), &
            H5T_NATIVE_REAL, file_struct%arr_dspace(1), &
            dset_id, hdf_error)
       file_struct%arr_dset(1, num_dataset, file_struct%num_group) = dset_id
       !write data
       !+ to horizontal layers
       dims = (/ dimX_arr, dimY_arr, num_vlayers /)
       CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, REAL(arrObj(:,:,idx_vlayers(1:num_vlayers))), dims, hdf_error)
       
       WRITE(dataset_name, "('hdataset-', i3.3)") num_dataset
       !put to the current group
       CALL h5dcreate_f(file_struct%arr_group(file_struct%num_group), TRIM(dataset_name), &
            H5T_NATIVE_REAL, file_struct%arr_dspace(2), &
            dset_id, hdf_error)
       file_struct%arr_dset(2, num_dataset, file_struct%num_group) = dset_id
       !+ to vertical layers
       dims = (/ num_hlayers, dimY_arr, dimZ_arr /)
       CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, REAL(arrObj(idx_hlayers(1:num_hlayers),:,:)), dims, hdf_error)
       
       file_struct%num_dset(file_struct%num_group) = num_dataset
       
       num_time_point = (file_struct%num_group-1) * MAX_NUM_HANDLE + num_dataset
       file_struct%time_data(num_time_point) = time
    END IF

    !check for full file
    IF (file_struct%num_group .EQ. MAX_NUM_GROUPS .AND. file_struct%num_dset(file_struct%num_group) .eq. MAX_NUM_HANDLE) THEN
       is_full_file = .TRUE.
    ENDIF

  END SUBROUTINE hdf_spatial_write

  !! Close a file
  SUBROUTINE hdf_spatial_close(file_struct)
    IMPLICIT NONE
    TYPE(HDFFILE2), INTENT(INOUT) :: file_struct
    INTEGER :: hdf_error, ii, jj
    INTEGER, PARAMETER :: RANK = 1
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    CHARACTER(100) :: dataset_name
    INTEGER :: num_time_point
    INTEGER(HID_T) :: space_id, dset_id
    
    !!create an attribute for root (/) group
    !CALL h5acreate_f("/",)
    !create a dataset for root (/) group that contains time points
    ! now we know how much data to record
    num_time_point = (file_struct%num_group-1) * MAX_NUM_HANDLE + file_struct%num_dset(file_struct%num_group)
    current_dims = (/ num_time_point /)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    dataset_name = "/time"
    CALL h5dcreate_f(file_struct%file_ID, TRIM(dataset_name), &
         H5T_NATIVE_REAL, space_id, dset_id, hdf_error)
    CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, REAL(file_struct%time_data(1:num_time_point)), current_dims, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)
    CALL h5dclose_f(dset_id, hdf_error)
    
    !close PL
    DO ii = 1, file_struct%num_plist
       CALL h5pclose_f(file_struct%arr_plist(ii), hdf_error)
    END DO
    !close dataset + group
    DO jj = 1, file_struct%num_group
       DO ii = 1, file_struct%num_dset(jj)
          CALL h5dclose_f(file_struct%arr_dset(1, ii, jj), hdf_error)
          CALL h5dclose_f(file_struct%arr_dset(2, ii, jj), hdf_error)
       END DO
       !close group
       CALL h5gclose_f(file_struct%arr_group(jj), hdf_error)
    END DO
    !close dataspace
    DO ii = 1, file_struct%num_dspace
       CALL h5sclose_f(file_struct%arr_dspace(ii), hdf_error)
    END DO
    
    CALL h5fclose_f(file_struct%file_ID, hdf_error)

    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0
    file_struct%num_group = 0
  END SUBROUTINE hdf_spatial_close
  
  !!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!
  !!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!-!!
  !! deal with spatial cell: write out whole cell
  SUBROUTINE hdf_spatial_open_2(file_struct, filename, dimX_arr, dimY_arr, dimZ_arr)
    !create file, 1 data-spaces
    IMPLICIT NONE
    TYPE(HDFFILE3), INTENT(INOUT) :: file_struct
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: dimX_arr, dimY_arr, dimZ_arr
    ! constant
    INTEGER, PARAMETER :: RANK = 3 ! 3D-dataspace
    ! local
    INTEGER (HID_T) :: file_ID, space_id, group_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER :: hdf_error, group_index
    CHARACTER(100) :: group_name
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0
    file_struct%time_data = 0.0d0
    
    !create the HDF5 file
    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F, H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !setting up dataspace
    ! We need only a single data spaces
    !+create the data space (fixed size) 
    current_dims = (/ dimX_arr, dimY_arr, dimZ_arr /) !size at creation time
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    file_struct%num_dspace =  1
    file_struct%arr_dspace(1) = space_id
  END SUBROUTINE hdf_spatial_open_2

  SUBROUTINE  hdf_spatial_write_2(file_struct, arrObj, &
       dimX_arr, dimY_arr, dimZ_arr, time, is_Full_File)
    !! data is organized as a series of datasets:
    !! 1 data set '/time'
    !! a series of dataset '/data-xxx' to contain either cmyo, cnsr, caf, depending
    !!             upon the data
    IMPLICIT NONE
    TYPE(HDFFILE3), INTENT(INOUT) :: file_struct    
    REAL(KIND=dp), DIMENSION(:,:,:), INTENT(IN) :: arrObj
    INTEGER, INTENT(IN):: dimX_arr, dimY_arr, dimZ_arr
    REAL(KIND=dp), INTENT(IN) :: time
    LOGICAL, INTENT(INOUT) :: is_full_file
    
    !local
    CHARACTER(100):: dataset_name, group_name
    INTEGER, PARAMETER :: rank = 3
    INTEGER(HSIZE_T), DIMENSION(rank) :: dims
    INTEGER(HID_T) :: dset_id, group_id
    INTEGER :: hdf_error, group_index
    INTEGER :: num_time_point, num_dataset

    num_dataset = file_struct%num_dset
    IF (num_dataset .LT. MAX_NUM_HANDLE) THEN
       num_dataset = num_dataset + 1
       WRITE(dataset_name, "('data-', i3.3)") num_dataset
       CALL h5dcreate_f(file_struct%file_ID, TRIM(dataset_name), &
            H5T_NATIVE_REAL, file_struct%arr_dspace(1), &
            dset_id, hdf_error)
       file_struct%arr_dset(num_dataset) = dset_id
       file_struct%num_dset = num_dataset
       !write data
       !+ to horizontal layers
       dims = (/ dimX_arr, dimY_arr, dimZ_arr /)
       CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, REAL(arrObj), dims, hdf_error)
       
       num_time_point = num_dataset
       file_struct%time_data(num_time_point) = time
    END IF

    !check for full file
    IF (file_struct%num_dset .eq. MAX_NUM_HANDLE) THEN
       is_full_file = .TRUE.
    ENDIF

  END SUBROUTINE hdf_spatial_write_2

  !! Close a file
  SUBROUTINE hdf_spatial_close_2(file_struct)
    IMPLICIT NONE
    TYPE(HDFFILE3), INTENT(INOUT) :: file_struct
    INTEGER :: hdf_error, ii, jj
    INTEGER, PARAMETER :: RANK = 1
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    CHARACTER(100) :: dataset_name
    INTEGER :: num_time_point
    INTEGER(HID_T) :: space_id, dset_id
    
    !create a dataset for root (/) group that contains time points
    ! now we know how much data to record
    num_time_point = file_struct%num_dset
    current_dims = (/ num_time_point /)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    dataset_name = "/time"
    CALL h5dcreate_f(file_struct%file_ID, TRIM(dataset_name), &
         H5T_NATIVE_REAL, space_id, dset_id, hdf_error)
    CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, REAL(file_struct%time_data(1:num_time_point)), current_dims, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)
    CALL h5dclose_f(dset_id, hdf_error)
    
    !close PL
    DO ii = 1, file_struct%num_plist
       CALL h5pclose_f(file_struct%arr_plist(ii), hdf_error)
    END DO
    !close dataset 
    DO ii = 1, file_struct%num_dset
       CALL h5dclose_f(file_struct%arr_dset(ii), hdf_error)
    END DO
   
    !close dataspace
    DO ii = 1, file_struct%num_dspace
       CALL h5sclose_f(file_struct%arr_dspace(ii), hdf_error)
    END DO
    
    CALL h5fclose_f(file_struct%file_ID, hdf_error)

    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0
  END SUBROUTINE hdf_spatial_close_2

  SUBROUTINE hdf_ecc_open(filename, NSFU, CHUNK_SIZE, CHUNKS_PER_FILE, file_struct, numfields)
    ! NSFU = number of SFU 
    ! CHUNK_SIZE = the number of record to keep before flushing out (A.K.A. a chunk)
    ! CHUNKS_PER_FILE = number of chunk in a file
    ! [OUT]file_struct = structure that hold all information about the HDF file
    ! [OUT]numfields = number of data per record 
    IMPLICIT NONE
    CHARACTER*(*) :: filename
    INTEGER, INTENT(IN) :: NSFU
    INTEGER, INTENT(IN) :: CHUNK_SIZE, CHUNKS_PER_FILE
    TYPE(HDFFILE), INTENT(INOUT) :: file_struct
    INTEGER, INTENT(IN) :: numfields
    ! local
    INTEGER (HID_T) :: file_ID
    INTEGER :: hdf_error
    !init
    file_struct%num_dspace = 0
    file_struct%num_dset = 0
    file_struct%num_plist = 0

    !create the HDF5 file

    CALL h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, file_ID, hdf_error, &
         H5P_DEFAULT_F)
    file_struct%file_ID = file_ID
  END SUBROUTINE hdf_ecc_open

  SUBROUTINE hdf_ecc_write(file_struct, data_buffer, numrecords, numfields, NSFU)
    IMPLICIT NONE
    INTEGER, PARAMETER :: rank = 2
    TYPE(HDFFILE) :: file_struct
    REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: data_buffer
    INTEGER, INTENT(IN) :: numrecords, numfields
    INTEGER, INTENT(IN) :: NSFU
    !local
    INTEGER (HID_T) ::  space_id, dset_id
    INTEGER(HSIZE_T), DIMENSION(RANK) :: current_dims
    INTEGER(HSIZE_T), DIMENSION(rank) :: dims
    INTEGER :: hdf_error, num_columns    


    !setting up dataspace + dataset
    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !! general data first
    num_columns = numfields - 3 * NSFU
    current_dims = (/ num_columns, numrecords /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    CALL h5dcreate_f(file_struct%file_ID, "/general", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    dims = current_dims
    CALL h5dwrite_f(dset_ID, &
         H5T_NATIVE_REAL, REAL(data_buffer(1:num_columns,1:numrecords)), &
         dims, hdf_error)
    CALL h5dclose_f(dset_ID, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)

    !! !!!
    !! SFU_state
    num_columns = NSFU
    current_dims = (/ num_columns, numrecords /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    CALL h5dcreate_f(file_struct%file_ID, "/sfu-state", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    dims = current_dims
    CALL h5dwrite_f(dset_ID, &
         H5T_NATIVE_REAL, REAL(data_buffer(numfields-3*NSFU+1:numfields-2*NSFU,1:numrecords)), &
         dims, hdf_error)
    CALL h5dclose_f(dset_ID, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)

    !! !!!!!!!!!!!
    !! Cds
    num_columns = NSFU
    current_dims = (/ num_columns, numrecords /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    CALL h5dcreate_f(file_struct%file_ID, "/cds", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    dims = current_dims
    CALL h5dwrite_f(dset_ID, &
         H5T_NATIVE_REAL, REAL(data_buffer(numfields-2*NSFU+1:numfields-NSFU,1:numrecords)), &
         dims, hdf_error)
    CALL h5dclose_f(dset_ID, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)

    !! !!!!!!!!!!!!!!!
    !! Cjsr
    num_columns = NSFU
    current_dims = (/ num_columns, numrecords /) !size at creation time
    !+create data space (extendable, unlimited)
    CALL h5screate_simple_f(RANK, current_dims, space_id, &
         hdf_error)
    CALL h5dcreate_f(file_struct%file_ID, "/cjsr", H5T_NATIVE_REAL, space_id, &
                 dset_id, hdf_error)
    dims = current_dims
    CALL h5dwrite_f(dset_ID, &
         H5T_NATIVE_REAL, REAL(data_buffer(numfields-NSFU+1:numfields,1:numrecords)), &
         dims, hdf_error)
    CALL h5dclose_f(dset_ID, hdf_error)
    CALL h5sclose_f(space_id, hdf_error)

  END SUBROUTINE hdf_ecc_write

  SUBROUTINE hdf_ecc_close(file_struct)
    IMPLICIT NONE
    TYPE(HDFFILE) :: file_struct
    INTEGER :: hdf_error, ii

    DO ii = 1, file_struct%num_plist
       CALL h5pclose_f(file_struct%arr_plist(ii), hdf_error)
    END DO
    DO ii = 1, file_struct%num_dset
       CALL h5dclose_f(file_struct%arr_dset(ii), hdf_error)
    END DO
    DO ii = 1, file_struct%num_dspace
       CALL h5sclose_f(file_struct%arr_dspace(ii), hdf_error)
    END DO
    CALL h5fclose_f(file_struct%file_ID, hdf_error)
  END SUBROUTINE hdf_ecc_close

END MODULE hdf5_utility
