#include <stdio.h>
#include <stdlib.h>
#include <gp.h>
#include "isofile.h"
#include "mymalloc.h"


static void cutout( float *** c, int nx, int ny, int nz, int ncut)
{
  int i,j,k;
  if (ncut<=0) return;

  for(i=nx/2;i < nx;i++)
    for(j=0; j < ny;j++)
      for(k=0; k < nz;k++)
	{
	  if(i > nx - ncut - 1)
	    c[i][j][k]  = 0.;
	  if(j > ny - ncut - 1) 
	    c[i][j][k]  = 0.;
	  if(k> nz - ncut - 1) 
	    c[i][j][k]  = 0.;
	}
  return;
}

int main(int argc, char *argv[])
{
  int pprank=0;
  int parallel_output;
  int nx_proc,ny,nz,np,nporig;
  char filename[256];
  int ismerged;
  int mx,my,mz;
  int t;
  float *pflc, *flc, ***level;

  if (argc != 2) {
    fprintf(stderr, "usage: %s time < parfile.par \n", argv[0]);
    exit(1);
  }
  t = atoi(argv[1]); 
  nx_proc = echo_ifetch("nx",0,"SEVERE");
  ny = echo_ifetch("ny",0,"SEVERE");
  nz = echo_ifetch("nz",0,"SEVERE");
  nporig = echo_ifetch("np",0,"SEVERE");
  np = echo_ifetch("npcut",nporig,"DEFAULT");
  ismerged =  echo_ifetch("ismerged",0,"WARNING");
  parallel_output =  echo_ifetch("parallel_output",0,"WARNING");

  /*   if files where merged they where saved in a single datb file, so parallel output =0 . 
  if files where not merged they may originate from a parallel computation 
(so parallel output =1)   
 or from a scalar computation (parallel output =0).   
 In the last case, files must be treated as not 
 merged to remove the the outer band */

  if(ismerged == 1) parallel_output=0;

  mx = (nx_proc - 2)*np;
  my = ny - 2;
  mz = nz - 2;
  
  pflc = flc = (float *) malloc(mx * my * mz * sizeof(float));
  if(pflc == NULL) nrerror("zmerge: cannot allocate work array");

 	
  if(ismerged == 0 ) {
    for(pprank=0;pprank<np;pprank++)	{
      if(parallel_output == 1) 
	sprintf(filename, "%s.%d.%d", "plot.datb",t,pprank);
      else
	sprintf(filename, "%s.%d", "plot.datb",t);
      inigetflc(np,pflc,filename,ismerged);
      pflc += (nx_proc - 2)*my*mz;
    }
  }
  else    {
    sprintf(filename, "%s.%d", "plot.datb",t);
    inigetflc(np,pflc,filename,ismerged);
  }

#if 1
    /*  constructs a 3D array using the memory space allocated to flc */
      if (!(level = (float ***)mymalloc3D_with_p(flc,sizeof(float), mx, my, mz))) 	{
	  fprintf(stderr, "zmerge: cannot allocate array level\n");
	  exit(1);
	}
 #else
      if (!(level = (float ***)mymalloc3D(sizeof(float), mx, my, mz))) 	{
	  fprintf(stderr, "zmerge: cannot allocate array level\n");
	  exit(1);
	}
      { 	int i,j,k;
	for(i=0;i<mx;i++)
	  for(j=0;j<my;j++)
	    for(k=0;k<mz;k++)
	      level[i][j][k]  = (float) flc[i*my*mz + mz*j + k];
      }
 #endif
  cutout(level,mx,my,mz,echo_ifetch("cutout",0,"DEFAULT"));
  if(echo_ifetch("isolated",0,"DEFAULT")==1) OutputIsolated(level,mx,my,mz);
  callisofile(level,mx,my,mz); 
	
  return 1;
}