#include <stdio.h>
#include <stdlib.h>
#include <gp.h>

#ifdef CRAY
void MACHINEID()
#else
#ifdef LINUX
void machineid_()
#else
void machineid()
#endif
#endif
{
/*   fprintf(stderr,"Running on: %s\n",hostname());*/
  system("hostname > hostname");
}

#ifdef LINUX
int srfintmax_()
#else
int srfintmax()
#endif
{
  return (int)(~0U>>1);
}

void datbfile_(double *u, double *v, double *c, double *p, int *nx, int *ny, int *nz, int *t, double *sigma, double *tau)
{
  int i, j, k, nvar = 4;
  float val;
  double h;
  char s[256];
  int n1,n2;
  FILE *fptr;
  double radius = GP_dfetch("radius",0.,"SEVERE"); 

  h=1./(*nx - 2);

  sprintf(s, "%s.%d", "plot.datb",*t);
  fptr = fopen(s, "wb");

  n1 = *nx-2;
  fwrite(&n1, sizeof(int), 1, fptr);
  n2 = *ny-2;
  fwrite(&n2, sizeof(int), 1, fptr);
  fwrite(&*nz, sizeof(int), 1, fptr);
  fwrite(&nvar, sizeof(int), 1, fptr);

  sprintf(s, "u");
  fwrite(s, sizeof(char), 256, fptr);
  sprintf(s, "v");
  fwrite(s, sizeof(char), 256, fptr);
  sprintf(s, "c");
  fwrite(s, sizeof(char), 256, fptr);
  sprintf(s, "p");
  fwrite(s, sizeof(char), 256, fptr);

  for (i = 1; i < *nx-1; i++) 
    for (j = 1; j < *ny-1; j++)
      for (k = 0; k < *nz; k++) {
        val = (u[i + (*nx)*j + (*nx)*(*ny)*k]+u[i + (*nx)*j + (*nx)*(*ny)*k + 1])/2;
        fwrite(&val, sizeof(float), 1, fptr);
      }

  for (i = 1; i < *nx-1; i++)  
    for (j = 1; j < *ny-1; j++)      
      for (k = 0; k < *nz; k++) {
        val = (v[i + (*nx)*j + (*nx)*(*ny)*k]+v[i + (*nx)*(j+1) + (*nx)*(*ny)*k])/2;
        fwrite(&val, sizeof(float), 1, fptr);
      }

  
    for (i = 1; i < *nx-1; i++)
      for (j = 1; j < *ny-1; j++) 
	for (k = 0; k < *nz; k++){
        val = c[i + (*nx)*j + (*nx)*(*ny)*k];
        fwrite(&val, sizeof(float), 1, fptr);
      }            

    /* u = u-nondim tau / h  -> p = p-nondim tau / h^2 */
  
    for (i = 1; i < *nx-1; i++)
      for (j = 1; j < *ny-1; j++) 
	for (k = 0; k < *nz; k++){
	  /*
for spurious current calculations
        val= (float) (h*h*radius * (p[i + (*nx)*j + (*nx)*(*ny)*k] - p[2 +  (*nx)*2])
	  /(*sigma * *tau * *tau) ) ;
	  */
        val=(float)(p[i + (*nx)*j + (*nx)*(*ny)*k] - p[2 + (*nx)*2]);
        fwrite(&val, sizeof(float), 1, fptr);
      }            

  fclose(fptr);
}