/*! This routine allocates memory for particle storage, both the
 *  collisionless and the SPH particles.
 */
void allocate_memory(void)
{
  size_t bytes;
  double bytes_tot = 0;

  if(All.MaxPart > 0)
    {
      if(!(P = malloc(bytes = All.MaxPart * sizeof(struct particle_data))))
	{
	  printf("failed to allocate memory for `P' (%g MB).\n", bytes / (1024.0 * 1024.0));
	  endrun(1);
	}
      bytes_tot += bytes;

      if(ThisTask == 0)
	printf("\nAllocated %g MByte for particle storage. %d\n\n", bytes_tot / (1024.0 * 1024.0), sizeof(struct particle_data));
    }

  if(All.MaxPartSph > 0)
    {
      bytes_tot = 0;

      if(!(SphP = malloc(bytes = All.MaxPartSph * sizeof(struct sph_particle_data))))
	{
	  printf("failed to allocate memory for `SphP' (%g MB) %d.\n", bytes / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
	  endrun(1);
	}
      bytes_tot += bytes;

      if(ThisTask == 0)
	printf("Allocated %g MByte for storage of SPH data. %d\n\n", bytes_tot / (1024.0 * 1024.0), sizeof(struct sph_particle_data));
    }
}

/* Load code and prototypes from file */
static void load(char *_filename) {
    FILE *fp;
    char *cp;
    word n, left;
    char filename[100]; /* should suffice on all systems */

    strcpy(filename, _filename);
    M = mallocate (memorysize + 1); /* allocate main memory array */
    if (M == NULL) abend("Memory size too large (use /m option)\n");

    addext(filename, ".ccd");
    if ((fp = fopen(filename, BINARYREAD)) == NULL) {
        fprintf(stderr, "Cannot open .ccd file\n");
        endrun(10);
    };

    ic = 0; /* read static data and code */
    left = memorysize + 1; /* from .ccd file */
    do {
        if (left == 0) abend("Memory size too small (use /m option)\n");
        n = min (IOBLOCK / sizeof(word), left);
        n = fread((char *) &M[ic], sizeof(word), (int) n, fp);
        ic += n;
        left -= n;
    } while (n != 0); /* now ic = number of words read */

    fclose(fp);
    /* Get various addresses passed by GENERATOR */
    ipradr = M[ic - 5]; /* primitive type desctriptions */
    temporary = M[ic - 4]; /* global temporary variables */
    strings = M[ic - 3]; /* string constants */
    lastprot = M[ic - 2]; /* last prototype number */
    freem = M[ic - 1]; /* first free word in memory */

    /* Read prototypes from .pcd file */
    addext(filename, ".pcd");
    if ((fp = fopen(filename, BINARYREAD)) == NULL) {
        fprintf(stderr, "Cannot open .pcd file\n");
        endrun(10);
    }
    for (n = MAINBLOCK; n <= lastprot; n++) {
        cp = ballocate (sizeof(protdescr));
        if (cp == NULL) abend("Memory size too large (use /m option)\n");
        prototype[n] = (protdescr *) cp;
        if (fread(cp, sizeof(protdescr), 1, fp) != 1)
            abend("Cannot read .pcd file\n");
    }
    fclose(fp);

    /* Open trace file */
    if (debug) {
        addext(filename, ".trd");
        if ((tracefile = fopen(filename, "w")) == NULL)
            abend("Cannot open .trd file\n");
    }
} /* end load */

/*!  This function opens various log-files that report on the status and
 *   performance of the simulstion. On restart from restart-files
 *   (start-option 1), the code will append to these files.
 */
void open_outputfiles(void)
{
  char mode[2], buf[200];

  if(ThisTask != 0)		/* only the root processor writes to the log files */
    return;

  if(RestartFlag == 0)
    strcpy(mode, "w");
  else
    strcpy(mode, "a");


  sprintf(buf, "%s%s", All.OutputDir, All.CpuFile);
  if(!(FdCPU = fopen(buf, mode)))
    {
      printf("error in opening file '%s'\n", buf);
      endrun(1);
    }

  sprintf(buf, "%s%s", All.OutputDir, All.InfoFile);
  if(!(FdInfo = fopen(buf, mode)))
    {
      printf("error in opening file '%s'\n", buf);
      endrun(1);
    }

  sprintf(buf, "%s%s", All.OutputDir, All.EnergyFile);
  if(!(FdEnergy = fopen(buf, mode)))
    {
      printf("error in opening file '%s'\n", buf);
      endrun(1);
    }

  sprintf(buf, "%s%s", All.OutputDir, All.TimingsFile);
  if(!(FdTimings = fopen(buf, mode)))
    {
      printf("error in opening file '%s'\n", buf);
      endrun(1);
    }

#ifdef FORCETEST
  if(RestartFlag == 0)
    {
      sprintf(buf, "%s%s", All.OutputDir, "forcetest.txt");
      if(!(FdForceTest = fopen(buf, "w")))
	{
	  printf("error in opening file '%s'\n", buf);
	  endrun(1);
	}
      fclose(FdForceTest);
    }
#endif
}

/*! This function allocates the memory neeed to compute the long-range PM
 *  force. Three fields are used, one to hold the density (and its FFT, and
 *  then the real-space potential), one to hold the force field obtained by
 *  finite differencing, and finally a workspace field, which is used both as
 *  workspace for the parallel FFT, and as buffer for the communication
 *  algorithm used in the force computation.
 */
void pm_init_periodic_allocate(int dimprod)
{
  static int first_alloc = 1;
  int dimprodmax;
  double bytes_tot = 0;
  size_t bytes;

  MPI_Allreduce(&dimprod, &dimprodmax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);

  /* allocate the memory to hold the FFT fields */
  #ifdef FFTW3
  if(!(rhogrid = fftw_alloc_real(bytes = fftsize_real)))
  #else
  if(!(rhogrid = (fftw_real *) malloc(bytes = fftsize * sizeof(fftw_real))))
  #endif
    {
      printf("failed to allocate memory for `FFT-rhogrid' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  #ifdef FFTW3
  if(!(forcegrid = fftw_alloc_real(bytes = imax(fftsize_real, dimprodmax) * 2)))
  #else
  if(!(forcegrid = (fftw_real *) malloc(bytes = imax(fftsize, dimprodmax) * sizeof(fftw_real))))
  #endif
    {
      printf("failed to allocate memory for `FFT-forcegrid' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  #ifdef FFTW3
  if(!(workspace = fftw_alloc_real(bytes = imax(maxfftsize, dimprodmax) * 2)))
  #else
  if(!(workspace = (fftw_real *) malloc(bytes = imax(maxfftsize, dimprodmax) * sizeof(fftw_real))))
  #endif
    {
      printf("failed to allocate memory for `FFT-workspace' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  if(first_alloc == 1)
    {
      first_alloc = 0;
      if(ThisTask == 0)
	printf("\nAllocated %g MByte for FFT data.\n\n", bytes_tot / (1024.0 * 1024.0));
    }

  fft_of_rhogrid = (fftw_complex *) & rhogrid[0];
}

void *mymalloc(size_t n)
{
  if((n % 8) > 0)
    n = (n / 8 + 1) * 8;

  if(n < 8)
    n = 8;

  if(Nblocks >= MAXBLOCKS)
    {
      printf("Task=%d: No blocks left in mymalloc().\n", ThisTask);
      endrun(813);
    }

#ifdef PEDANTIC_MEMORY_CEILING
  if(n > FreeBytes)
    {
      printf("Task=%d: Not enough memory in mymalloc(n=%g MB).  FreeBytes=%g MB\n",
	     ThisTask, n / (1024.0 * 1024.0), FreeBytes / (1024.0 * 1024.0));
      endrun(812);
    }
  Table[Nblocks] = Base + (TotBytes - FreeBytes);
  FreeBytes -= n;
#else
  Table[Nblocks] = malloc(n);
  if(!(Table[Nblocks]))
    {
      printf("failed to allocate %g MB of memory. (presently allocated=%g MB)\n",
	     n / (1024.0 * 1024.0), AllocatedBytes / (1024.0 * 1024.0));
      endrun(18);
    }
#endif

  AllocatedBytes += n;
  BlockSize[Nblocks] = n;

  Nblocks += 1;

  /*
     if(AllocatedBytes / (1024.0 * 1024.0) > Highmark)
     {
     Highmark = AllocatedBytes / (1024.0 * 1024.0);
     printf("Task=%d:   new highmark=%g MB\n",
     ThisTask, Highmark);
     fflush(stdout);
     }
   */

  return Table[Nblocks - 1];
}

void *myrealloc(void *p, size_t n)
{
  if((n % 8) > 0)
    n = (n / 8 + 1) * 8;

  if(n < 8)
    n = 8;

  if(Nblocks == 0)
    endrun(76879);

  if(p != Table[Nblocks - 1])
    {
      printf("Task=%d: Wrong call of myrealloc() - not the last allocated block!\n", ThisTask);
      fflush(stdout);
      endrun(815);
    }

  AllocatedBytes -= BlockSize[Nblocks - 1];
#ifdef PEDANTIC_MEMORY_CEILING
  FreeBytes += BlockSize[Nblocks - 1];
#endif


#ifdef PEDANTIC_MEMORY_CEILING
  if(n > FreeBytes)
    {
      printf("Task=%d: Not enough memory in myremalloc(n=%g MB). previous=%g FreeBytes=%g MB\n",
	     ThisTask, n / (1024.0 * 1024.0), BlockSize[Nblocks - 1] / (1024.0 * 1024.0),
	     FreeBytes / (1024.0 * 1024.0));
      endrun(812);
    }
  Table[Nblocks - 1] = Base + (TotBytes - FreeBytes);
  FreeBytes -= n;
#else
  Table[Nblocks - 1] = realloc(Table[Nblocks - 1], n);
  if(!(Table[Nblocks - 1]))
    {
      printf("failed to reallocate %g MB of memory. previous=%g FreeBytes=%g MB\n",
	     n / (1024.0 * 1024.0), BlockSize[Nblocks - 1] / (1024.0 * 1024.0),
	     FreeBytes / (1024.0 * 1024.0));
      endrun(18123);
    }
#endif

  AllocatedBytes += n;
  BlockSize[Nblocks - 1] = n;

  return Table[Nblocks - 1];
}

int main(int argc, char **argv)
{
  double s1;
  int i, irank, nrank;

  ARGC = argc;
  ARGV = argv;

  OUTPUT=0;
  if(argc==1)
    endrun("./QPM.mock qpm.bat_file > output");

  read_parameter_file(argv[1]);

  SIGMA_8Z0 = SIGMA_8;
  SIGMA_8 = SIGMA_8*growthfactor(REDSHIFT);
  fprintf(stdout,"SIGMA_8(Z=%.2f)= %.3f\n",REDSHIFT,SIGMA_8);
  RESET_COSMOLOGY++;

  if(argc>2)
    {
      if(atoi(argv[2])==999)
	test();
      else
	SUBFRAC = atof(argv[2]);
    }
  
  if(Task.create_halos)
    create_lognormal_halos();
  if(Task.populate_simulation)
    populate_simulation_hod();
  exit(0);

}

static inline void perp(const double x[], const double y[], double out[])
{
  // out is a unit vector perpendicular to both x and y
  double oo;
  out[0]= x[1]*y[2] - x[2]*y[1];
  out[1]= x[2]*y[0] - x[0]*y[2];
  out[2]= x[0]*y[1] - x[1]*y[0];
  
  oo= sqrt(out[0]*out[0] + out[1]*out[1] + out[2]*out[2]);

  if(oo == 0.0) {
    fprintf(stderr, "x= %e %e %e\n", x[0], x[1], x[2]);
    fprintf(stderr, "y= %e %e %e\n", y[0], y[1], y[2]);
    endrun(4003);
  }

  out[0] /= oo;
  out[1] /= oo;
  out[2] /= oo;

  /*
  if(out[0] != out[0] || out[1] != out[1] || out[2] != out[2])
    endrun(4004);
  */
}

void ReadIonizeParams(char *fname)
{
  int i;
  FILE *fdcool;

  if(!(fdcool = fopen(fname, "r")))
    {
      printf(" Cannot read ionization table in file `%s'\n", fname);
      endrun(456);
    }

  for(i = 0; i < TABLESIZE; i++)
    gH0[i] = 0;

  for(i = 0; i < TABLESIZE; i++)
    if(fscanf(fdcool, "%g %g %g %g %g %g %g",
	      &inlogz[i], &gH0[i], &gHe[i], &gHep[i], &eH0[i], &eHe[i], &eHep[i]) == EOF)
      break;

  fclose(fdcool);

  /*  nheattab is the number of entries in the table */

  for(i = 0, nheattab = 0; i < TABLESIZE; i++)
    if(gH0[i] != 0.0)
      nheattab++;
    else
      break;

  if(DEBUG)
    printf("\n\nread ionization table with %d entries in file `%s'.\n\n", nheattab, fname);
}

/*! Allocates memory for the neighbour list buffer.
 */
void ngb_treeallocate(int npart)
{
  double totbytes = 0;
  size_t bytes;

#ifdef PERIODIC
  boxSize = All.BoxSize;
  boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
  boxHalf_X = boxHalf * LONG_X;
  boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
  boxHalf_Y = boxHalf * LONG_Y;
  boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
  boxHalf_Z = boxHalf * LONG_Z;
  boxSize_Z = boxSize * LONG_Z;
#endif
#endif

  if(!(Ngblist = malloc(bytes = npart * (long) sizeof(int))))
    {
      printf("Failed to allocate %g MB for ngblist array\n", bytes / (1024.0 * 1024.0));
      endrun(78);
    }
  totbytes += bytes;

  if(ThisTask == 0)
    printf("allocated %g Mbyte for ngb search.\n", totbytes / (1024.0 * 1024.0));
}

void mymalloc_init(void)
{
#ifdef PEDANTIC_MEMORY_CEILING
  size_t n;
#endif

  BlockSize = (size_t *) malloc(MAXBLOCKS * sizeof(size_t));
  Table = (void **) malloc(MAXBLOCKS * sizeof(void *));

#ifdef PEDANTIC_MEMORY_CEILING
  n = PEDANTIC_MEMORY_CEILING * 1024.0 * 1024.0;

  if(!(Base = malloc(n)))
    {
      printf("Failed to allocate memory for `Base' (%d Mbytes).\n", (int) PEDANTIC_MEMORY_CEILING);
      endrun(122);
    }

  TotBytes = FreeBytes = n;
#endif

  AllocatedBytes = 0;
  Nblocks = 0;
  Highmark = 0;
}

/*! This function determines how many particles there are in a given block,
 *  based on the information in the header-structure.  It also flags particle
 *  types that are present in the block in the typelist array. If one wants to
 *  add a new output-block, this function should be augmented accordingly.
 */
int get_particles_in_block(enum iofields blocknr, int *typelist)
{
  int i, nall, ntot_withmasses, ngas, nstars;

  nall = 0;
  ntot_withmasses = 0;

  for(i = 0; i < 6; i++)
    {
      typelist[i] = 0;

      if(header.npart[i] > 0)
	{
	  nall += header.npart[i];
	  typelist[i] = 1;
	}

      if(All.MassTable[i] == 0)
	ntot_withmasses += header.npart[i];
    }

  ngas = header.npart[0];
  nstars = header.npart[4];


  switch (blocknr)
    {
    case IO_POS:
    case IO_VEL:
    case IO_ACCEL:
    case IO_TSTP:
    case IO_ID:
    case IO_POT:
      return nall;
      break;

    case IO_MASS:
      for(i = 0; i < 6; i++)
	{
	  typelist[i] = 0;
	  if(All.MassTable[i] == 0 && header.npart[i] > 0)
	    typelist[i] = 1;
	}
      return ntot_withmasses;
      break;

    case IO_U:
    case IO_RHO:
    case IO_HSML:
    case IO_DTENTR:
      for(i = 1; i < 6; i++)
	typelist[i] = 0;
      return ngas;
      break;
    }

  endrun(212);
  return 0;
}

/* This routine allocates memory for 
 * particle storage, both the collisionless and the SPH particles.
 * The memory for the ordered binary tree of the timeline
 * is also allocated.
 */
void allocate_memory_2d(void)
{
  int bytes,bytes_tot=0;

  printf("MaxPart %d\n",All.MaxPart);
  if(All.MaxPart>0)
    {
      if(!(Pn_data=malloc(bytes=All.MaxPart*sizeof(struct particle_data))))
	{
	  printf("failed to allocate memory for `Pn_data' (%d bytes).\n",bytes);
	  endrun(1);
	}
      bytes_tot+=bytes;


      if(!(PTimeTree=malloc(bytes=All.MaxPart*sizeof(struct timetree_data))))
	{
	  printf("failed to allocate memory for `PTimeTree' (%d bytes).\n",bytes);
	  endrun(1);
	}
      bytes_tot+=bytes;


      Pn = Pn_data-1;   /* start with offset 1 */
      PTimeTree--;

      printf("\nAllocated %g MByte for particle storage.\n\n",bytes_tot/(1024.0*1024.0));
    }

  if(All.MaxPartSph>0)
    {
      bytes_tot=0;

      if(!(SphPn_data=malloc(bytes=All.MaxPartSph*sizeof(struct  sph_particle_data))))
	{
	  printf("failed to allocate memory for `SphPn_data' (%d bytes).\n",bytes);
	  endrun(1);
	}
      bytes_tot+=bytes;

      SphPn= SphPn_data-1;   /* start with offset 1 */

      printf("Allocated %g MByte for storage of SPH data.\n\n",bytes_tot/(1024.0*1024.0));
    }
}

int MPI_Sizelimited_Sendrecv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
			     int dest, int sendtag, void *recvbuf, int recvcount,
			     MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm,
			     MPI_Status * status)
{
  int iter = 0, size_sendtype, size_recvtype, send_now, recv_now;
  int count_limit;


  if(dest != source)
    endrun(3);

  MPI_Type_size(sendtype, &size_sendtype);
  MPI_Type_size(recvtype, &size_recvtype);

  if(dest == ThisTask)
    {
      memcpy(recvbuf, sendbuf, recvcount * size_recvtype);
      return 0;
    }

  count_limit = (int) ((((long long) MPISENDRECV_SIZELIMIT) * 1024 * 1024) / size_sendtype);

  while(sendcount > 0 || recvcount > 0)
    {
      if(sendcount > count_limit)
	{
	  send_now = count_limit;
	  if(iter == 0)
	    {
	      printf("imposing size limit on MPI_Sendrecv() on task=%d (send of size=%d)\n",
		     ThisTask, sendcount * size_sendtype);
	      fflush(stdout);
	    }
	  iter++;
	}
      else
	send_now = sendcount;

      if(recvcount > count_limit)
	recv_now = count_limit;
      else
	recv_now = recvcount;

      MPI_Sendrecv(sendbuf, send_now, sendtype, dest, sendtag,
		   recvbuf, recv_now, recvtype, source, recvtag, comm, status);

      sendcount -= send_now;
      recvcount -= recv_now;

      sendbuf += send_now * size_sendtype;
      recvbuf += recv_now * size_recvtype;
    }

  return 0;
}

void myfree_msg(void *p, char *msg)
{
  if(Nblocks == 0)
    endrun(76878);

  if(p != Table[Nblocks - 1])
    {
      printf("Task=%d: Wrong call of myfree() - '%s' not the last allocated block!\n", ThisTask, msg);
      fflush(stdout);
      endrun(8141);
    }

  Nblocks -= 1;
  AllocatedBytes -= BlockSize[Nblocks];
#ifdef PEDANTIC_MEMORY_CEILING
  FreeBytes += BlockSize[Nblocks];
#else
  free(p);
#endif
}

/*! This catches I/O errors occuring for my_fwrite(). In this case we
 *  better stop.
 */
size_t my_fwrite(void *ptr, size_t size, size_t nmemb, FILE * stream)
{
  size_t nwritten;

  if((nwritten = fwrite(ptr, size, nmemb, stream)) != nmemb)
    {
      printf("I/O error (fwrite) on task=%d has occured: %s\n", ThisTask, strerror(errno));
      fflush(stdout);
      endrun(777);
    }
  return nwritten;
}

/*! This catches I/O errors occuring for fread(). In this case we
 *  better stop.
 */
size_t my_fread(void *ptr, size_t size, size_t nmemb, FILE * stream)
{
  size_t nread;

  if((nread = fread(ptr, size, nmemb, stream)) != nmemb)
    {
      printf("I/O error (fread) on task=%d has occured: %s\n", ThisTask, strerror(errno));
      fflush(stdout);
      endrun(778);
    }
  return nread;
}

double esys_half(double r)
{
  static int flag=1,n;
  static double *x,*y,*y2;
  FILE *fp;
  int i;
  double e,m,b;
  float x1,x2,x3,x4,x5;

  if(Work.SysErrFlag==0)return(0);

  return 0.03;
  return 0.04;

  if(flag)
    {
      flag=0;
      if(!(fp=fopen(Work.esysfile,"r")))
	{
	  fprintf(stdout,"ERROR opening [%s]\n",Work.esysfile);
	  endrun("error opening esys file");
	}
      n=filesize(fp);
      x=dvector(1,n);
      y=dvector(1,n);
      y2=dvector(1,n);

      for(i=1;i<=n;++i)
	{
	  fscanf(fp,"%f %f %f %f %f",&x1,&x2,&x3,&x4,&x5);
	  x[i]=x4;
	  y[i]=x5;
	  /*
	  y[i]*=5./3.;
	  if(y[i]<0.07)y[i]=0.07;
	  */
	}
      spline(x,y,n,1.0E+30,1.0E+30,y2);
      fclose(fp);
    }
  /*
  for(i=1;i<=n;++i)
    if(r>x[i])break;
  if(i==1)return(y[1]);
  if(i>=n)return(y[n]);
  m=(y[i-1]-y[i])/(x[i-1]-x[i]);
  b=y[i]-m*x[i];
  return(m*r+b);
  */
  splint(x,y,y2,n,r,&e);  
  return(e);
}

/*! This function allocates the workspace needed for the non-periodic FFT
 *  algorithm. Three fields are used, one for the density/potential fields,
 *  one to hold the force field obtained by finite differencing, and finally
 *  an additional workspace which is used both in the parallel FFT itself, and
 *  as a buffer for the communication algorithm.
 */
void pm_init_nonperiodic_allocate(int dimprod)
{
  static int first_alloc = 1;
  int dimprodmax;
  double bytes_tot = 0;
  size_t bytes;

  MPI_Allreduce(&dimprod, &dimprodmax, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);

  if(!(rhogrid = (fftw_real *) malloc(bytes = fftsize * sizeof(fftw_real))))
    {
      printf("failed to allocate memory for `FFT-rhogrid' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  fft_of_rhogrid = (fftw_complex *) rhogrid;

  if(!(forcegrid = (fftw_real *) malloc(bytes = imax(fftsize, dimprodmax) * sizeof(fftw_real))))
    {
      printf("failed to allocate memory for `FFT-forcegrid' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  if(!(workspace = (fftw_real *) malloc(bytes = imax(maxfftsize, dimprodmax) * sizeof(fftw_real))))
    {
      printf("failed to allocate memory for `FFT-workspace' (%g MB).\n", bytes / (1024.0 * 1024.0));
      endrun(1);
    }
  bytes_tot += bytes;

  if(first_alloc == 1)
    {
      first_alloc = 0;
      if(ThisTask == 0)
	printf("\nUsing %g MByte for non-periodic FFT computation.\n\n", bytes_tot / (1024.0 * 1024.0));
    }
}

/*! Allocates memory for the neighbour list buffer.
 */
void ngb_treeallocate(int npart)
{
  double totbytes = 0;
  size_t bytes;

#ifdef PERIODIC
  boxSize = All.BoxSize;
  boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
  boxHalf_X = boxHalf * LONG_X;
  boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
  boxHalf_Y = boxHalf * LONG_Y;
  boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
  boxHalf_Z = boxHalf * LONG_Z;
  boxSize_Z = boxSize * LONG_Z;
#endif
#endif

#ifdef _OPENMP
  if(MAXTHREADS < omp_get_max_threads()){
    printf("Can't allocate memory, there will be an overrun! Change MAXTHREADS\n");
    endrun(788888);
  }
#endif
  if(!(Ngblist = malloc(bytes = MAXTHREADS * npart * (long) sizeof(int)))) //npart = MAXN
    {
      printf("Failed to allocate %g MB for ngblist array\n", bytes / (1024.0 * 1024.0));
      endrun(78);
    }
  totbytes += bytes;

  if(ThisTask == 0)
    printf("allocated %g Mbyte for ngb search. %d \n", totbytes / (1024.0 * 1024.0),MAXTHREADS);
}