/* 	$Id: myhdf5.h,v 1.4 2007/07/26 20:35:25 cvsictr1 Exp $	 */
#ifndef lint
static char myhdf5_vcid[] = "$Id: myhdf5.h,v 1.4 2007/07/26 20:35:25 cvsictr1 Exp $";
#endif /* lint */

#include <stdarg.h>

/* Global variables */
hid_t my_file,eu_file;
hid_t my_group,eu_group;
hid_t my_dataset,eu_dataset;
hid_t my_dataset_idx;
#define TRUE 1
#define FALSE 0

#define MYHDF5_VERSION 1

#define BubbleBaseType float

typedef struct {
  BubbleBaseType x,  y,  z;
} tri_field;

typedef struct {
  BubbleBaseType dxvx,  dxvy,  dxvz,  
                 dyvx,  dyvy,  dyvz,  
                 dzvx,  dzvy,  dzvz;
} grad_field;

uint *_idx=NULL;
void *_internals=NULL;

#define HDF5_NOT_SET 0x000000

/* Here globals */
int   DNS_SIZE=HDF5_NOT_SET;
float DNS_FORCING=0.0;
float PHYS_TAU_ETA=0.0;
float PHYS_RELAMBDA=0.0;
float PHYS_RE=0.0;
float PHYS_ETA=0.0;
float PHYS_EPS=0.0;
float PHYS_U_RMS_X=0.0;
float PHYS_U_RMS_Y=0.0;
float PHYS_U_RMS_Z=0.0;
float PHYS_L0=0.0;
float LAGR_T0=0.0;
float LAGR_DT=0.0;
float LAGR_TAU=0.0;
float LAGR_BETA=0.0;
float LAGR_STOKES=0.0;
int LAGR_PROPERTIES_PROP=0;
int LAGR_PROPERTIES_BASE=0;
int LAGR_IDX_TYPE=0;
int LAGR_TOTAL_TYPES=0;
int LAGR_THERMALIZED=0;

int LAGR_SLOWST=0;            /* number of slow particles per Stokes */
int SLOW_CUBE_SIZE =0;   /* size of the side of the cube in which slow particles are organized */
int SLOW_CUBE_NUM =0;    /* number of cubes (in one direction) in which slow particles are organized*/
int SLOW_LINEAR_SIZE =0; /* total number of such cubes*/
int EULER_CUBE_SIZE =0;   /* size of the side of the cube in which euler data are organized */
int EULER_CUBE_NUM =0;    /* number of cubes (in one direction) in which euler data are organized*/
int EULER_LINEAR_SIZE =0; /* total number of such cubes*/

/* ATTT: TODO
char EULER_ORIGINAL_FILE_TIME_STAMP[128];
*/



int myhdf5_iread(char *aname)
{
  hid_t  attribute;
  herr_t status;
  int tmp;
  /* Mandatory check of consistency, in order to satisfy regression tests */
  attribute = H5Aopen_name(my_group,aname);
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n",aname);
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &tmp);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n",aname);
    }
    status = H5Aclose(attribute);
  }
  return tmp;
}

float myhdf5_fread(char *aname)
{
  hid_t  attribute;
  herr_t status;
  float tmp;
  attribute = H5Aopen_name(my_group,aname);
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n",aname);
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &tmp);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n",aname);
    }
    status = H5Aclose(attribute);
  }
  return tmp;
}



int myhdf5_ireadEU(char *aname)
{
  hid_t  attribute;
  herr_t status;
  int tmp;
  /* Mandatory check of consistency, in order to satisfy regression tests */
  attribute = H5Aopen_name(eu_group,aname);
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n",aname);
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &tmp);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n",aname);
    }
    status = H5Aclose(attribute);
  }
  return tmp;
}

float myhdf5_freadEU(char *aname)
{
  hid_t  attribute;
  herr_t status;
  float tmp;
  attribute = H5Aopen_name(eu_group,aname);
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n",aname);
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &tmp);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n",aname);
    }
    status = H5Aclose(attribute);
  }
  return tmp;
}

#define fetch_int(A) A=myhdf5_iread(#A)
#define fetch_float(A) A=myhdf5_fread(#A)

#define fetchEU_int(A) A=myhdf5_ireadEU(#A)
#define fetchEU_float(A) A=myhdf5_freadEU(#A)

int myhdf5_open(char *fname, char *dname)
{
  herr_t (*old_func)(void*);
  void *old_client_data;

  char name[128];
  char *my_string;
  H5G_stat_t statbuf;

  hid_t  attribute;
  herr_t      status;
  H5Z_EDC_t   edc;
  hsize_t num_obj;
  int n, len, ival;

  H5Eget_auto(&old_func, &old_client_data);

  H5Eset_auto(NULL, NULL); 
  my_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
  my_group = H5Gopen(my_file,"DNS");
  sprintf(name,"DNS/%s",dname);
  my_dataset = H5Dopen(my_file, name);

  if (my_dataset<0) {
    fprintf(stderr,"Error opening dataset %s\n",name);
    return -1;
  }

  /* Mandatory check of consistency, in order to satisfy regression tests */
  attribute = H5Aopen_name(my_group,"MYHDF5_VERSION");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find MYHDF5_VERSION attribute in group.\n");
    return -2;
  }
  status = H5Aread(attribute, H5T_NATIVE_INT, &ival);  
#ifdef DEBUG_MYHDF5
  fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival);
#endif
  if (status<0) {
    fprintf(stderr,"Error cannot find attribute MYHDF5_VERSION\n");
    return -2;
  }
  status = H5Aclose(attribute);

  /* fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival); */

  if (ival>MYHDF5_VERSION) {
    fprintf(stderr,"ERROR: the version of myhdf5.h (%d) is too old for the dataset (requires >= %d)\n",MYHDF5_VERSION,ival);
    return -3;
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_TAU");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_TAU");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_TAU);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_TAU");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_BETA");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_BETA");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_BETA);
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_BETA");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_DT");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot fine %s attribute in group.\n","LAGR_DT");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_DT);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_DT");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_T0");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot fine %s attribute in group.\n","LAGR_T0");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_T0);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_T0");
    }
    status = H5Aclose(attribute);
  }

  fetch_float(DNS_FORCING);
  fetch_int(DNS_SIZE);
  fetch_float(PHYS_TAU_ETA);
  fetch_float(PHYS_RELAMBDA);
  fetch_float(PHYS_RE);
  fetch_float(PHYS_ETA);
  fetch_float(PHYS_EPS);
  fetch_float(PHYS_U_RMS_X);
  fetch_float(PHYS_U_RMS_Y);
  fetch_float(PHYS_U_RMS_Z);
  fetch_float(PHYS_L0);
  /*
    fetch_float(LAGR_T0);
    fetch_float(LAGR_DT);
  */

  H5Eset_auto(old_func, old_client_data);
  return 0;
}




int myhdf5_open_euler(char *fname, char *dname)
{
  herr_t (*old_func)(void*);
  void *old_client_data;
  char name[128],name_idx[128];
  char *my_string;
  H5G_stat_t statbuf;

  hid_t  attribute;
  herr_t      status;
  H5Z_EDC_t   edc;
  hsize_t num_obj;
  int n, len, ival;
  

  H5Eget_auto(&old_func, &old_client_data);

  H5Eset_auto(NULL, NULL); 
  eu_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
  eu_group = H5Gopen(eu_file,"DNS");
  sprintf(name,"DNS/%s",dname);
  eu_dataset = H5Dopen(eu_file, name);


  if (eu_dataset<0) {
    fprintf(stderr,"Error opening dataset %s\n",name);
    return -1;
  }

  /* Mandatory check of consistency, in order to satisfy regression tests */
  attribute = H5Aopen_name(eu_group,"MYHDF5_VERSION");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot fine MYHDF5_VERSION attribute in group.\n");
    return -2;
  }
  status = H5Aread(attribute, H5T_NATIVE_INT, &ival);  
#ifdef DEBUG_MYHDF5
  fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival);
#endif
  if (status<0) {
    fprintf(stderr,"Error cannot find attribute MYHDF5_VERSION\n");
    return -2;
  }
  status = H5Aclose(attribute);

  /* fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival); */

  if (ival>MYHDF5_VERSION) {
    fprintf(stderr,"ERROR: the version of myhdf5.h (%d) is too old for the dataset (requires >= %d)\n",MYHDF5_VERSION,ival);
    return -3;
  }

 
  attribute = H5Aopen_name(eu_dataset,"EULER_CUBE_SIZE");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","EULER_CUBE_SIZE");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &EULER_CUBE_SIZE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","EULER_CUBE_SIZE");
    }
    status = H5Aclose(attribute);
  }
 
  attribute = H5Aopen_name(eu_dataset,"EULER_CUBE_NUM");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","EULER_CUBE_NUM");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &EULER_CUBE_NUM);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","EULER_CUBE_NUM");
    }
    status = H5Aclose(attribute);
  }
 
  attribute = H5Aopen_name(eu_dataset,"EULER_LINEAR_SIZE");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","EULER_LINEAR_SIZE");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &EULER_LINEAR_SIZE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","EULER_LINEAR_SIZE");
    }
    status = H5Aclose(attribute);
  }

  /* ATTT::::  TODO  ::::::::
  attribute = H5Aopen_name(eu_dataset,"EULER_ORIGINAL_FILE_TIME_STAMP");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","EULER_ORIGINAL_FILE_TIME_STAMP");
  } else {
    status = H5Aread(attribute, ?????, &EULER_ORIGINAL_FILE_TIME_STAMP);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","EULER_ORIGINAL_FILE_TIME_STAMP");
    }
    status = H5Aclose(attribute);
  }
  */
  

  fetchEU_float(DNS_FORCING);
  fetchEU_int(DNS_SIZE);
  fetchEU_float(PHYS_TAU_ETA);
  fetchEU_float(PHYS_RELAMBDA);
  fetchEU_float(PHYS_RE);
  fetchEU_float(PHYS_ETA);
  fetchEU_float(PHYS_EPS);
  fetchEU_float(PHYS_U_RMS_X);
  fetchEU_float(PHYS_U_RMS_Y);
  fetchEU_float(PHYS_U_RMS_Z);
  fetchEU_float(PHYS_L0);

  H5Eset_auto(old_func, old_client_data);
  return 0;
}





/* this subroutine open the dataset slow */
int myhdf5_open_slow(char *fname, char *dname)
{
  herr_t (*old_func)(void*);
  void *old_client_data;
  char name[128],name_idx[128];
  char *my_string;
  H5G_stat_t statbuf;

  hid_t  attribute;
  herr_t      status;
  H5Z_EDC_t   edc;
  hsize_t num_obj;
  int n, len, ival;
  

  H5Eget_auto(&old_func, &old_client_data);

  H5Eset_auto(NULL, NULL); 
  my_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
  my_group = H5Gopen(my_file,"DNS");
  sprintf(name,"DNS/%s",dname);
  my_dataset = H5Dopen(my_file, name);
  sprintf(name_idx,"DNS/%s_IDX",dname);
  my_dataset_idx = H5Dopen(my_file, name_idx);

  if (my_dataset<0) {
    fprintf(stderr,"Error opening dataset %s\n",name);
    return -1;
  }
  if (my_dataset_idx<0) {
    fprintf(stderr,"Error opening dataset %s\n",name_idx);
    return -1;
  }

  /* Mandatory check of consistency, in order to satisfy regression tests */
  attribute = H5Aopen_name(my_group,"MYHDF5_VERSION");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot fine MYHDF5_VERSION attribute in group.\n");
    return -2;
  }
  status = H5Aread(attribute, H5T_NATIVE_INT, &ival);  
#ifdef DEBUG_MYHDF5
  fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival);
#endif
  if (status<0) {
    fprintf(stderr,"Error cannot find attribute MYHDF5_VERSION\n");
    return -2;
  }
  status = H5Aclose(attribute);

  /* fprintf(stderr,"MYHDF5_VERSION requested >= %d\n",ival); */

  if (ival>MYHDF5_VERSION) {
    fprintf(stderr,"ERROR: the version of myhdf5.h (%d) is too old for the dataset (requires >= %d)\n",MYHDF5_VERSION,ival);
    return -3;
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_TAU");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_TAU");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_TAU);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_TAU");
    }
    status = H5Aclose(attribute);
  }


  attribute = H5Aopen_name(my_dataset,"LAGR_PROPERTIES_PROP");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_PROPERTIES_PROP");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_PROPERTIES_PROP);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_PROPERTIES_PROP");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_TAU");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_TAU");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_TAU);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_TAU");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_BETA");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_BETA");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_BETA);
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_BETA");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_STOKES");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_STOKES");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT, &LAGR_STOKES);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_STOKES");
    }
    status = H5Aclose(attribute);
  }


  attribute = H5Aopen_name(my_dataset,"LAGR_PROPERTIES_BASE");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_PROPERTIES_BASE");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_PROPERTIES_BASE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_PROPERTIES_BASE");
    }
    status = H5Aclose(attribute);
  }


  

  attribute = H5Aopen_name(my_dataset,"LAGR_IDX_TYPE");
  if (attribute<0) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_IDX_TYPE");
#endif
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_IDX_TYPE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_IDX_TYPE");
    }
    status = H5Aclose(attribute);
  }
 

  attribute = H5Aopen_name(my_dataset,"LAGR_TOTAL_TYPES");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_TOTAL_TYPES");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_TOTAL_TYPES);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_TOTAL_TYPES");
    }
    status = H5Aclose(attribute);
  }
 

 
  attribute = H5Aopen_name(my_dataset,"LAGR_T0");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_T0");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT,  &LAGR_T0);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n"," LAGR_T0");
    }
    status = H5Aclose(attribute);
  }
 
 
  attribute = H5Aopen_name(my_dataset,"LAGR_DT");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_DT");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_FLOAT,  &LAGR_DT);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_DT");
    }
    status = H5Aclose(attribute);
  }
 

  attribute = H5Aopen_name(my_dataset,"LAGR_THERMALIZED");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_THERMALIZED");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_THERMALIZED);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_THERMALIZED");
    }
    status = H5Aclose(attribute);
  }

  attribute = H5Aopen_name(my_dataset,"LAGR_SLOWST");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","LAGR_SLOWST");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &LAGR_SLOWST);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","LAGR_SLOWST");
    }
    status = H5Aclose(attribute);
  }
 
  attribute = H5Aopen_name(my_dataset,"SLOW_CUBE_SIZE");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","SLOW_CUBE_SIZE");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &SLOW_CUBE_SIZE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","SLOW_CUBE_SIZE");
    }
    status = H5Aclose(attribute);
  }
 
  attribute = H5Aopen_name(my_dataset,"SLOW_CUBE_NUM");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","SLOW_CUBE_NUM");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &SLOW_CUBE_NUM);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","SLOW_CUBE_NUM");
    }
    status = H5Aclose(attribute);
  }
 
  attribute = H5Aopen_name(my_dataset,"SLOW_LINEAR_SIZE");
  if (attribute<0) {
    fprintf(stderr,"ERROR cannot find %s attribute in group.\n","SLOW_LINEAR_SIZE");
  } else {
    status = H5Aread(attribute, H5T_NATIVE_INT, &SLOW_LINEAR_SIZE);  
    if (status<0) {
      fprintf(stderr,"Error cannot read attribute %s\n","SLOW_LINEAR_SIZE");
    }
    status = H5Aclose(attribute);
  }
 
  _idx = (uint *) malloc(sizeof(uint)*SLOW_LINEAR_SIZE);
  status = H5Dread( my_dataset_idx, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, _idx);
  
  fetch_float(DNS_FORCING);
  fetch_int(DNS_SIZE);
  fetch_float(PHYS_TAU_ETA);
  fetch_float(PHYS_RELAMBDA);
  fetch_float(PHYS_RE);
  fetch_float(PHYS_ETA);
  fetch_float(PHYS_EPS);
  fetch_float(PHYS_U_RMS_X);
  fetch_float(PHYS_U_RMS_Y);
  fetch_float(PHYS_U_RMS_Z);
  fetch_float(PHYS_L0);
  /*
    fetch_float(LAGR_T0);
    fetch_float(LAGR_DT);
  */

  H5Eset_auto(old_func, old_client_data);
  return 0;
}



const static void * _MY_VA_CLOSE_TAG_= (void *) 0xfede;
#define myhdf5_read_cube(a,b,c,...) va_myhdf5_read_cube(a,b,c,__VA_ARGS__,_MY_VA_CLOSE_TAG_)

int va_myhdf5_read_cube(int ix,int iy,int iz, ...)
{
  /* Must deal with vargs */
  va_list args;
  tri_field **p   = NULL;
  tri_field **u   = NULL;
  tri_field **v   = NULL;
  grad_field **dv = NULL;
  tri_field **Dv = NULL;
  hid_t filespace;
  hid_t memspace;
  hid_t hdf5_pos;
  hid_t hdf5_u;
  hid_t hdf5_v;
  hid_t hdf5_dv;
  hid_t hdf5_Dv;
  hsize_t offset[1], count[1];
  hsize_t moffset[1], mcount[1];
  herr_t      status, status_n, herr;

  char name[128];
  int icube,Cstart;
  int npart;

  int x=0;
  void *val;
  
  /* variadic walk-through */
  va_start( args, iz );
  while ((val= (void *) va_arg( args, double ** )) != _MY_VA_CLOSE_TAG_ ) {
    if (x==0)  p = ( tri_field **) val;
    if (x==1)  u = ( tri_field **) val;
    if (x==2)  v = ( tri_field **) val;
    if (x==3) dv = (grad_field **) val;
    if (x==4) Dv = (tri_field **) val;
    x++;
  }
  va_end(args);


  //  fprintf(stderr,"%d\n",x);

#ifdef DEBUG_MYHDF5
  if (p  != NULL) fprintf(stderr,"INFO: requested to fetch  p  tri_field\n");
  if (u  != NULL) fprintf(stderr,"INFO: requested to fetch  u  tri_field\n");
  if (v  != NULL) fprintf(stderr,"INFO: requested to fetch  v  tri_field\n");
  if (dv != NULL) fprintf(stderr,"INFO: requested to fetch dv grad_field\n");
  if (Dv != NULL) fprintf(stderr,"INFO: requested to fetch Dv tri_field\n");
#endif

  if ((ix<SLOW_CUBE_NUM) && (iy<SLOW_CUBE_NUM) && (iz<SLOW_CUBE_NUM)) {
    icube = iz + SLOW_CUBE_NUM*(iy+SLOW_CUBE_NUM*ix);  
    if (icube<SLOW_LINEAR_SIZE-1) {
      Cstart = *(_idx+icube);
      npart  = *(_idx+icube+1) - Cstart;
    } else {
      Cstart = *(_idx+icube);
      npart =  (LAGR_SLOWST-1)-Cstart;
    }
#ifdef DEBUG_MYHDF5
        fprintf(stderr,"Cube[%d][%d][%d]=%d start at %d and has %d PART\n",ix,iy,iz,icube,Cstart,npart); 
#endif
  } else {
    fprintf(stderr,"ERR: Cube you asked for DOES NOT EXIST!!\n");
  }
#ifdef DEBUG_MYHDF5
  fprintf(stderr,"npart %d cstart %d\n", npart, Cstart);
#endif   
  offset[0]=Cstart; count[0]=npart;
  moffset[0]=0;    mcount[0]=npart;

  filespace = H5Dget_space(my_dataset);
  memspace  = H5Dget_space(my_dataset);
  H5Sselect_hyperslab(memspace, H5S_SELECT_SET, moffset, NULL, mcount, NULL);
  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);
  
  if (p==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"p not requested\n");
#endif
  } else {
    *p = (tri_field *) malloc(sizeof(tri_field)*npart); 
    hdf5_pos = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_pos, "x", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_pos, "y", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_pos, "z", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
    status = H5Dread(my_dataset, hdf5_pos , memspace, filespace, H5P_DEFAULT, *p);
  }

  if (u==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"u not requested\n");
#endif
  } else {
    *u = (tri_field *) malloc(sizeof(tri_field)*npart); 
    hdf5_u = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_u, "ux", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_u, "uy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_u, "uz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
    status = H5Dread(my_dataset, hdf5_u , memspace, filespace, H5P_DEFAULT, *u);
  }

  if (v==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"v not requested\n");
#endif
  } else {
    *v = (tri_field *) malloc(sizeof(tri_field)*npart); 
    hdf5_v = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_v, "vx", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_v, "vy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_v, "vz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
    status = H5Dread(my_dataset, hdf5_v , memspace, filespace, H5P_DEFAULT, *v);
  }

  if (dv==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"dv not requested\n");
#endif
  } else {
    *dv = (grad_field *) malloc(sizeof(grad_field)*npart); 
    hdf5_dv = H5Tcreate (H5T_COMPOUND, sizeof(grad_field));
    H5Tinsert(hdf5_dv, "dxvx", HOFFSET(grad_field, dxvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dxvy", HOFFSET(grad_field, dxvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dxvz", HOFFSET(grad_field, dxvz), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dyvx", HOFFSET(grad_field, dyvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dyvy", HOFFSET(grad_field, dyvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dyvz", HOFFSET(grad_field, dyvz), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dzvx", HOFFSET(grad_field, dzvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dzvy", HOFFSET(grad_field, dzvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dzvz", HOFFSET(grad_field, dzvz), H5T_NATIVE_FLOAT);
    status = H5Dread(my_dataset, hdf5_dv , memspace, filespace, H5P_DEFAULT, *dv);
  }
 
  if (Dv==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"Dv not requested\n");
#endif
  } else {
    *Dv = (tri_field *) malloc(sizeof(tri_field)*npart); 
    hdf5_Dv = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_Dv, "Dvx", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_Dv, "Dvy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_Dv, "Dvz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
    status = H5Dread(my_dataset, hdf5_Dv , memspace, filespace, H5P_DEFAULT, *Dv);
  }


  if(p!=NULL)H5Tclose(hdf5_pos);
  if(u!=NULL)H5Tclose(hdf5_u);
  if(v!=NULL)H5Tclose(hdf5_v);
  if(dv!=NULL)H5Tclose(hdf5_dv);
  if(Dv!=NULL)H5Tclose(hdf5_Dv);

  //  H5Sclose(memspace);
  //  H5Sclose(filespace);
  return npart;
}


int myhdf5_euler_read_cube(int ix,int iy,int iz, tri_field **v)
{
  /* Must deal with vargs */
  hid_t filespace;
  hid_t memspace;
  hid_t hdf5_pos;
  hid_t hdf5_v;
  hsize_t offset[4], count[4];
  hsize_t moffset[4], mcount[4];
  herr_t  status, status_n, herr;

  char name[128];
  int icube,Cstart;
  int nfields;


  if ((ix<EULER_CUBE_NUM) && (iy<EULER_CUBE_NUM) && (iz<EULER_CUBE_NUM)) {
    icube = iz + EULER_CUBE_NUM*(iy+EULER_CUBE_NUM*ix);  
    nfields  = EULER_CUBE_SIZE*EULER_CUBE_SIZE*EULER_CUBE_SIZE;
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"Cube[%d][%d][%d]=%d  totfields=%d\n",ix,iy,iz,icube,nfields); 
#endif
  } else {
    fprintf(stderr,"ERR: Cube you asked for DOES NOT EXIST!!\n");
  }


  offset[0]=icube; 
  offset[1]=0;
  offset[2]=0;
  offset[3]=0;

  count[0]=1;
  count[1] = EULER_CUBE_SIZE;
  count[2] = EULER_CUBE_SIZE;
  count[3] = EULER_CUBE_SIZE;

  moffset[0]=0; 
  moffset[1]=0;
  moffset[2]=0;
  moffset[3]=0;

  mcount[0]=1;
  mcount[1] = EULER_CUBE_SIZE;
  mcount[2] = EULER_CUBE_SIZE;
  mcount[3] = EULER_CUBE_SIZE;


  filespace = H5Dget_space(eu_dataset);
  memspace  = H5Dget_space(eu_dataset);

  H5Sselect_hyperslab(memspace, H5S_SELECT_SET, moffset, NULL, mcount, NULL);
  H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);


  *v = (tri_field *) malloc(sizeof(tri_field)*nfields); 

  hdf5_v = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
  H5Tinsert(hdf5_v, "vx", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
  H5Tinsert(hdf5_v, "vy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
  H5Tinsert(hdf5_v, "vz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
  status = H5Dread(eu_dataset, hdf5_v , memspace, filespace, H5P_DEFAULT, *v);

  H5Tclose(hdf5_v);
  H5Sclose(memspace);
  H5Sclose(filespace);
  return nfields;
  }


/* This subroutine just print some information */
void myhdf5_info(void)
{
  char name[128];
  char *my_string;
  H5G_stat_t statbuf;

  herr_t      status, status_n, herr;
  H5Z_EDC_t   edc;
  hsize_t num_obj;
  int n, len;
  hid_t sid, tid;
  H5T_class_t t_class;

  sprintf(name,"DNS/St0");
  H5Gget_objinfo(my_file, name, FALSE, &statbuf);
  switch (statbuf.type) {
  case H5G_GROUP:
    printf("Object with name %s is a group \n", name); break;
  case H5G_DATASET:
    printf("Object with name %s is a dataset \n", name); break;
  case H5G_TYPE:
    printf("Object with name %s is a named datatype \n", name); break;
  case H5G_LINK:
    my_string = (char *) malloc(statbuf.linklen);

    H5Gget_linkval(my_file, name, statbuf.linklen, my_string);
    printf(" Object with name %s is a link to %s \n", name, my_string);
    H5Gget_objinfo(my_file, name, TRUE, &statbuf);
    switch (statbuf.type) {
    case H5G_GROUP:
      printf(" Target of link name %s is a group \n", name); break;
    case H5G_DATASET:
      printf(" Target of link name %s is a dataset \n", name); break;
    case H5G_TYPE:
      printf(" Target of link name %s is a named datatype \n", name); break;
    case H5G_LINK:
      printf(" Target of link name %s is a soft link \n", name); break;
    default:
      printf(" Unable to identify target ");
    }
    free(my_string); break;
  default:
    printf(" Unable to identify an object ");
  }

  herr  = H5Gget_num_objs(my_group,  &num_obj);
  fprintf(stderr,"Inside the group there is %d %d object\n",(int) num_obj,(int) herr);

  for (n=0;n<(int) num_obj;n++) {
    fprintf(stderr,"%d\n",n);
    len = H5Gget_objname_by_idx(my_group,(hsize_t) n, name,(size_t) 1024);
    fprintf(stderr,"%d %d %s\n",n,len,name);
  }


  sid = H5Dget_space(my_dataset); /* the dimensions of the dataset (not shown) */
  tid = H5Dget_type(my_dataset);

  t_class = H5Tget_class(tid);
  if (t_class < 0){
    puts(" Invalid datatype.\n");
  } else {
    if(t_class == H5T_COMPOUND) {
      puts("Datatype is 'H5T_COMPOUND'.\n");
    }
  }
  
  num_obj  = H5Tget_nmembers(tid);
  fprintf(stderr,"Inside the COMPOUND there are %d object\n",(int) num_obj);
  
  for (n=0; n<num_obj; n++) {
    my_string = H5Tget_member_name(tid, n);
    fprintf(stderr,"%d %s\n",n,my_string);
  }  
  H5Tclose(tid);
}

void myhdf5_close()
{
  H5Dclose(my_dataset);
  H5Gclose(my_group);
  H5Fclose(my_file);
}

void myhdf5_close_slow()
{
  if (_idx!=NULL) free(_idx);
  H5Dclose(my_dataset_idx);
  H5Dclose(my_dataset);
  H5Gclose(my_group);
  H5Fclose(my_file);
}



void myhdf5_close_euler()
{
  H5Dclose(eu_dataset);
  H5Gclose(eu_group);
  H5Fclose(eu_file);
}

/* #define _MY_VA_CLOSE_TAG_ 0xfede */
#define myhdf5_read(a,b,c,...) va_myhdf5_read(a,b,c,__VA_ARGS__,_MY_VA_CLOSE_TAG_)
void va_myhdf5_read(int *beam_traj, int *time_step, tri_field **p, ...)
{
  /* Must deal with vargs */
  va_list args;
  
  tri_field **u   = NULL;
  tri_field **v   = NULL;
  grad_field **dv = NULL;
  tri_field **Dv = NULL;

  char name[128];
  char *my_string;
  H5G_stat_t statbuf;

  double sum = 0;

  herr_t status, status_n;                             
  H5Z_EDC_t   edc;
  hsize_t num_obj;
  int n, len;
  hid_t filespace;

  int rank;
  hsize_t     dims[2];                     /* dataset and chunk dimensions*/ 
  hid_t hdf5_pos;
  hid_t hdf5_u;
  hid_t hdf5_v;
  hid_t hdf5_dv;
  hid_t hdf5_Dv;

  int x=0;
  void *val;
  
  /* variadic walk-through */
  va_start( args, p );
  while ((val= (void *) va_arg( args, double ** )) != _MY_VA_CLOSE_TAG_ ) {
    if (x==0)  u = ( tri_field **) val;
    if (x==1)  v = ( tri_field **) val;
    if (x==2) dv = (grad_field **) val;
    if (x==3) Dv = ( tri_field **) val;

    x++;
  }
  va_end(args);

  /* fprintf(stderr,"x=%d\n",x); */

#ifdef DEBUG_MYHDF5
  if (u  != NULL) fprintf(stderr,"INFO: requested to fetch  u  tri_field\n");
  if (v  != NULL) fprintf(stderr,"INFO: requested to fetch  v  tri_field\n");
  if (dv != NULL) fprintf(stderr,"INFO: requested to fetch dv grad_field\n");
  if (Dv != NULL) fprintf(stderr,"INFO: requested to fetch Dv  tri_field\n");
#endif

  filespace = H5Dget_space(my_dataset);    
  rank      = H5Sget_simple_extent_ndims(filespace);
  status_n  = H5Sget_simple_extent_dims(filespace, dims, NULL);

  *p = (tri_field *) malloc(sizeof(tri_field)*dims[0]*dims[1]); 

  *beam_traj = dims[1];
  *time_step = dims[0];
  
  hdf5_pos = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
  H5Tinsert(hdf5_pos, "x", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
  H5Tinsert(hdf5_pos, "y", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
  H5Tinsert(hdf5_pos, "z", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);
  status = H5Dread(my_dataset, hdf5_pos , H5S_ALL, H5S_ALL, H5P_DEFAULT, *p);
  
  if (u==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"u not requested\n");
#endif
  } else {
    
    *u = (tri_field *) malloc(sizeof(tri_field)*dims[0]*dims[1]); 
    hdf5_u = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_u, "ux", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_u, "uy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_u, "uz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);

    status = H5Dread(my_dataset, hdf5_u   , H5S_ALL, H5S_ALL, H5P_DEFAULT, *u);
  }

  if (v==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"v not requested\n");
#endif
  } else {

    *v = (tri_field *) malloc(sizeof(tri_field)*dims[0]*dims[1]); 
    if (*v == NULL) fprintf(stderr,"Problema\n");

    hdf5_v = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_v, "vx", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_v, "vy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_v, "vz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);

    status = H5Dread(my_dataset, hdf5_v   , H5S_ALL, H5S_ALL, H5P_DEFAULT, *v);
  }

  if (dv==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"dv not requested\n");
#endif
  } else {

    *dv = (grad_field *) malloc(sizeof(grad_field)*dims[0]*dims[1]);

    hdf5_dv = H5Tcreate (H5T_COMPOUND, sizeof(grad_field));
    H5Tinsert(hdf5_dv, "dxvx", HOFFSET(grad_field, dxvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dxvy", HOFFSET(grad_field, dxvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dxvz", HOFFSET(grad_field, dxvz), H5T_NATIVE_FLOAT);

    H5Tinsert(hdf5_dv, "dyvx", HOFFSET(grad_field, dyvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dyvy", HOFFSET(grad_field, dyvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dyvz", HOFFSET(grad_field, dyvz), H5T_NATIVE_FLOAT);

    H5Tinsert(hdf5_dv, "dzvx", HOFFSET(grad_field, dzvx), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dzvy", HOFFSET(grad_field, dzvy), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_dv, "dzvz", HOFFSET(grad_field, dzvz), H5T_NATIVE_FLOAT);

    status = H5Dread(my_dataset, hdf5_dv   , H5S_ALL, H5S_ALL, H5P_DEFAULT, *dv);
  }

  if (Dv==NULL) {
#ifdef DEBUG_MYHDF5
    fprintf(stderr,"Dv not requested\n");
#endif
  } else {

    *Dv = (tri_field *) malloc(sizeof(tri_field)*dims[0]*dims[1]);
    if (*Dv == NULL) fprintf(stderr,"Problema\n");

    hdf5_Dv = H5Tcreate (H5T_COMPOUND, sizeof(tri_field));
    H5Tinsert(hdf5_Dv, "Dvx", HOFFSET(tri_field, x), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_Dv, "Dvy", HOFFSET(tri_field, y), H5T_NATIVE_FLOAT);
    H5Tinsert(hdf5_Dv, "Dvz", HOFFSET(tri_field, z), H5T_NATIVE_FLOAT);

    status = H5Dread(my_dataset, hdf5_Dv   , H5S_ALL, H5S_ALL, H5P_DEFAULT, *Dv);
  }


 

  if ( p != NULL ) H5Tclose(hdf5_pos);
  if ( u != NULL ) H5Tclose(hdf5_u);
  if ( v != NULL ) H5Tclose(hdf5_v);
  if (dv != NULL ) H5Tclose(hdf5_dv);
  if (Dv != NULL ) H5Tclose(hdf5_Dv);
  //  H5Sclose(rank);
  //  H5Dclose(filespace);
}

void myhdf5_stat(int *beam_traj, int *time_step)
{
  /* Must deal with vargs */
  char name[128];
  char *my_string;
  H5G_stat_t statbuf;

  double sum = 0;

  herr_t status;
  hid_t filespace;

  int rank;
  hsize_t     dims[2];                     /* dataset and chunk dimensions*/ 

  int x=0;

  filespace = H5Dget_space(my_dataset);    
  rank      = H5Sget_simple_extent_ndims(filespace);
  status    = H5Sget_simple_extent_dims(filespace, dims, NULL);

  *beam_traj = dims[1];
  *time_step = dims[0];

  //  H5Sclose(rank);
  //  H5Dclose(filespace);
}