/* io.c
 *          *
 *        Emilio            4/18/2001
 *          In ReadT() and ReadTXYbl(), fixed the use of fseek(), so that the
 *          2nd argument is a position index multiplied by the size of the data
 *          type; previously, the data type size was not accounted for, and that
 *          may have been yielding incorrect data access results.
 *          *
 ******************************************************************************/

#include "emu.h"


/* functions defined in function.c */
extern BOOL GetFileName(VAR * pVar, char * pszDest);
extern void ReplaceDynAlias(char * pszDest, char * pszAlias, float fData,
                            emu_type EmuType);

/* LOCAL FUNCTIONS: Important "system" functions */
FILE *OpenFile(VAR var, int nYear, int nMonth, char * pszModeString, char * pszFile);


/* ================================================================= */
/* READ XY VARIABLE (maps/images), WRITE ARRAY INTO FLOAT POINTER pfData */

void ReadXY(float *pfData, VAR var, int year, int month)
{
  FILE *fp;
  BOOL bSwap = FALSE;
  int i, irev, col, row, nElements;      int bug=0;
  char szCurrentFile[FILENAMELEN];
  short *psData;
  float *x4bytetmp;
  unsigned short *pusData;
  unsigned char *pbyData;
  static size_t count[] = {0,0,0}, start[] = {0,0,0};
  char s[FILENAMELEN];


  nElements = var.nRows * var.nCols;

  /* ================== NETCDF FILE TYPE ====================== */
  if (istype(var.filetype, "nc"))
  {
    if (year == NONE && month == NONE)
    {
      count[0] = var.nRows;
      count[1] = var.nCols;
    }
    else
    {
      count[0] = 1;
      count[1] = var.nRows;
      count[2] = var.nCols;
    }

    /* set the time index */
    if (var.tTimeType == MONTHLY)
      start[0] = (12*year + month) - ((12*var.YearOrigin + 11) + (long)var.t1);
    else if (var.tTimeType == MONMEAN)
    {
      /*if (month < 0 || month > 11)
      {
        monthMult = (int) fmod((double)(month+1), (double)12);
        monthRem  = (month+1) % 12;
      } */
      start[0] = month;
    }
    else if (var.tTimeType == ANNUAL)
      start[0] = year - (var.YearOrigin + 1);
    else if (var.tTimeType == NOTIME)
      start[0] = 0;

    /* read array from the nc file */
    x4bytetmp = alloc1d_f(0,nElements-1);
    errnc = nc_get_vara_float(var.ncid, var.var_id, start, count, x4bytetmp);

    for (irev = 0; irev < nElements; irev++)
    {
      i = getRowRevCol(irev, &col, &row);
      pfData[i] = x4bytetmp[irev] * var.scale_factor + var.add_offset;
    }

    free1d_f(x4bytetmp,0,nElements-1);
  }
  /* ================== BINARY FILE TYPE ================== */
  else if (istype(var.filetype, "bin"))
  {
    bSwap = var.bByteSwap;

    /* construct file paths and names, open file */
    fp = OpenFile(var, year, month, "rb", szCurrentFile);
    if (bDebug)  printf("Opened binary file: %s\n", szCurrentFile);

    /* open the files, convert to float if necessary */
    if (var.tEmuType == EMU_FLOAT)
    {
      if (fread(pfData, sizeof(float), nElements, fp) != (size_t)(nElements))
        Warning("Could not do fread in %s \n", szCurrentFile);

      if (bSwap)  SwapBytes((BYTE *) pfData, 4, sizeof(float) * nElements);
    }
    else if (var.tEmuType == EMU_SHORT)
    {
      psData = alloc1d_s(0,nElements-1);
      if (fread(psData, sizeof(short), nElements, fp) != (size_t)(nElements))
        Warning("Could not do fread in %s \n", szCurrentFile);

      if (bSwap)  SwapBytes((BYTE *) psData, 2, sizeof(short) * nElements);
      for (i=0; i < nElements; i++)
        pfData[i] = (float) psData[i];

      free1d_s(psData,0,nElements-1);
    }
    else if (var.tEmuType == EMU_USHORT)
    {
      pusData = alloc1d_us(0,nElements-1);
      if (fread(pusData, sizeof(unsigned short), nElements, fp) != (size_t)(nElements))
        Warning("Could not do fread in %s \n", szCurrentFile);

      if (bSwap)  SwapBytes((BYTE *) pusData, 2, sizeof(short) * nElements);
      for (i=0; i < nElements; i++)
        pfData[i] = (float) pusData[i];

      free1d_us(pusData,0,nElements-1);
    }
    else if (var.tEmuType == EMU_BYTE)
    {
      pbyData = alloc1d_uc(0,nElements-1);
      if (fread(pbyData, sizeof(BYTE), nElements, fp) != (size_t)(nElements))
        Warning("Could not do fread in %s \n", szCurrentFile);

      for (i=0; i < nElements; i++)
        pfData[i] = (float) pbyData[i];

      free1d_uc(pbyData,0,nElements-1);
    }
    else
    {
      sprintf(s, "%d", var.tEmuType);
      Warning("Binary data type %s currently unsupported\n", s);
    }

    fclose(fp);
  }
  else if (istype(var.filetype, "ascii"))
    Warning("ASCII data format %s currently unsupported\n", var.filetype);
  else
    Warning("Data format %s currently unsupported\n", var.filetype);
}

/* ================================================================= */
/* READ T VARIABLE (time series), WRITE ARRAY INTO FLOAT POINTER pfData */
/* Note: This function is extremely inefficient, especially when reading sets
 * of individual XY binary files. In general, it's better to use ReadTXYbl() */

void ReadT(float *pfData, VAR var, long int i)
{
  FILE *binfin;
  BOOL bSwap = FALSE;
  int t, col, row, year, month, nElements;      int bug=0;
  char szCurrentFile[FILENAMELEN];
  short *ts2byte;
  float *ts4bytetmp;
  static size_t count[] = {0,0,0}, start[] = {0,0,0};


  nElements = nTimeSteps;

  /* ================== NETCDF FILE TYPE ====================== */
  if (istype(var.filetype, "nc"))
  {
    /* set the (reversed) row and (normal) col index */
    (void) getRowRevCol(i, &col, &row);
    row = (nRows - 1) - row;

    start[1] = row;
    start[2] = col;
    count[1] = 1;
    count[2] = 1;

    /* set the time index and count */
    if (var.tTimeType == MONTHLY)
    {
      start[0] = (12*nFirstYear + (nFirstMonth-1)) - ((12*var.YearOrigin + 11) + (long)var.t1);
      count[0] = nElements;
    }
    else if (var.tTimeType == MONMEAN)
    {
/* FOR THIS TYPE, THE start[0] and count[0] VALUES ARE COMPLICATED ... */
      start[0] = nFirstMonth-1;
      count[0] = YEAR2MONTH;
    }
    else if (var.tTimeType == ANNUAL)
    {
      start[0] = nFirstYear - (var.YearOrigin + 1);
      count[0] = nTotalYears;
    }

    /* read array from the nc file */
    ts4bytetmp = alloc1d_f(0,nElements-1);
    errnc = nc_get_vara_float(var.ncid, var.var_id, start, count, ts4bytetmp);

    for (t = 0; t < nElements; t++)
      pfData[t] = ts4bytetmp[t] * var.scale_factor + var.add_offset;

    free1d_f(ts4bytetmp,0,nElements-1);
  }
  /* ================== BINARY FILE TYPE ================== */
  else if (istype(var.filetype, "bin"))
  {
    bSwap = var.bByteSwap;

    /* open the files, convert to float if necessary */
    if (var.tEmuType == EMU_FLOAT)
    {
      for (t = 0; t < nElements; t++)
      {
        getyrmonth(t, &year, &month);
        binfin = OpenFile(var, year, month, "rb", szCurrentFile);
        if (fseek(binfin, i * sizeof(float), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", szCurrentFile);
        if (fread(&pfData[t], sizeof(float), 1, binfin) != 1)
          nrerror("Could not do fread in %s \n", szCurrentFile);
        fclose(binfin);

        if (bSwap)  SwapBytes((BYTE *) pfData, 4, sizeof(float) * nElements);
      }
    }
    else if (var.tEmuType == EMU_SHORT)
    {
      ts2byte = alloc1d_s(0,nElements-1);
      for (t = 0; t < nElements; t++)
      {
        getyrmonth(t, &year, &month);
        binfin = OpenFile(var, year, month, "rb", szCurrentFile);
        if (fseek(binfin, i * sizeof(short), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", szCurrentFile);
        if (fread(&ts2byte[t], sizeof(short), 1, binfin) != 1)
          nrerror("Could not do fread in %s \n", szCurrentFile);
        fclose(binfin);
      }

      if (bSwap)  SwapBytes((BYTE *) ts2byte, 2, sizeof(short) * nElements);
      for (t = 0; t < nElements; t++)
        pfData[t] = (float) ts2byte[t];

      free1d_s(ts2byte,0,nElements-1);
    }
  }
  else if (istype(var.filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", var.filetype);
  else
    Warning("Data type %s currently unsupported\n", var.filetype);
}

/* ================================================================= */
void ReadTXYbl(float **TsXYblock4byte, VAR var, int nyblock)
{
  FILE *binfin;
  int t, tvar, yr, mo, i, i0, irev, blockNCindx, col, row;  int bug=0;
  char currentfile[FILENAMELEN];
  short *TsXYblock2byte;
  float *TsXYblock4bytetmp;
  unsigned short *TsXYblocku2byte;
  unsigned char *TsXYblocku1byte;
  static size_t count[3], start[] = {0,0,0};


  /* ================== NETCDF FILE TYPE ====================== */
  if (istype(var.filetype, "nc"))
  {
    count[1] = nBlockRows;
    count[2] = nCols;
    start[1] = nRows - (nyblock+1)* nBlockRows;

    /* set the time index and count */
    if (var.tTimeType == MONTHLY)
    {
      start[0] = (12*nFirstYear + (nFirstMonth-1)) - ((12*var.YearOrigin + 11) + (long)var.t1);
      count[0] = nTimeSteps;
    }
    else if (var.tTimeType == MONMEAN)
    {
/* FOR THIS TYPE, THE start[0] and count[0] VALUES ARE COMPLICATED ... */
      start[0] = nFirstMonth-1;
      count[0] = YEAR2MONTH;
    }
    else if (var.tTimeType == ANNUAL)
    {
      start[0] = nFirstYear - (var.YearOrigin + 1);
      count[0] = nTotalYears;
    }


    if (bDebug)  printf("HERE: in readTseriesXYblock(), %d: start[0]=%d,start[1]=%d  \n", bug++,start[0],start[1]);

    /* read array from the nc file */
    TsXYblock4bytetmp = alloc1d_f(0,(count[0]*nBlockRows*nCols)-1);
    errnc = nc_get_vara_float(var.ncid, var.var_id, start, count, TsXYblock4bytetmp);

    i0 = nyblock*nBlockRows*nCols;
    for (t = 0; t < nTimeSteps; t++)
      for (i = i0; i < i0+nBlockRows*nCols; i++)
      {
        irev = getRowRevCol(i, &col, &row);
        getyrmonth(t, &yr, &mo);
        if (var.tTimeType == MONTHLY)
          tvar = t;
        else if (var.tTimeType == MONMEAN)
          tvar = mo;
        else if (var.tTimeType == ANNUAL)
          tvar = yr - (var.YearOrigin + 1);

        blockNCindx = (irev-start[1]*nCols) + (nBlockRows*nCols) * tvar;
        TsXYblock4byte[t][i-i0] = TsXYblock4bytetmp[blockNCindx] * var.scale_factor + var.add_offset;
      }   /* END i FOR-LOOP (and t for-loop afterwards) */

    free1d_f(TsXYblock4bytetmp,0,(count[0]*nBlockRows*nCols)-1);
  }
  /* ================== BINARY FILE TYPE ================== */
  else if (istype(var.filetype, "bin"))
  {
    if (var.tEmuType == EMU_FLOAT)
    {
      i0 = nyblock*nBlockRows*nCols;
      for (tvar = 0; tvar < nTimeSteps;tvar++)
      {
        getyrmonth(tvar, &yr, &mo);
        binfin = OpenFile(var, yr, mo, "r", currentfile);
        if (fseek(binfin, i0 * sizeof(float), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", currentfile);
        if (fread(TsXYblock4byte[tvar], sizeof(float), nBlockRows*nCols, binfin) != (size_t) nBlockRows*nCols)
          nrerror("Could not do fread in %s \n", currentfile);
        fclose(binfin);
      }
    }
    else if (var.tEmuType == EMU_SHORT)
    {
      TsXYblock2byte = alloc1d_s(0,nBlockRows*nCols-1);
      i0 = nyblock*nBlockRows*nCols;
      for (tvar = 0; tvar < nTimeSteps;tvar++)
      {
        getyrmonth(tvar, &yr, &mo);
        binfin = OpenFile(var, yr, mo, "r", currentfile);
        if (fseek(binfin, i0 * sizeof(short), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", currentfile);
        if (fread(TsXYblock2byte, sizeof(short), nBlockRows*nCols, binfin) != (size_t) nBlockRows*nCols)
          nrerror("Could not do fread in %s \n", currentfile);
        fclose(binfin);
        for (i = 0; i < nBlockRows*nCols; i++)
          TsXYblock4byte[tvar][i] = (float) TsXYblock2byte[i];
      }
      free1d_s(TsXYblock2byte,0,nBlockRows*nCols-1);
    }
    else if (var.tEmuType == EMU_BYTE)
    {
      TsXYblocku1byte = alloc1d_uc(0,nBlockRows*nCols-1);
      i0 = nyblock*nBlockRows*nCols;
      for (tvar = 0; tvar < nTimeSteps;tvar++)
      {
        getyrmonth(tvar, &yr, &mo);
        binfin = OpenFile(var, yr, mo, "r", currentfile);
        if (fseek(binfin, i0 * sizeof(BYTE), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", currentfile);
        if (fread(TsXYblocku1byte, sizeof(char), nBlockRows*nCols, binfin) != (size_t) nBlockRows*nCols)
          nrerror("Could not do fread in %s \n", currentfile);
        fclose(binfin);
        for (i = 0; i < nBlockRows*nCols; i++)
          TsXYblock4byte[tvar][i] = (float) TsXYblocku1byte[i];
      }
      free1d_uc(TsXYblocku1byte,0,nBlockRows*nCols-1);
    }
    else if (var.tEmuType == EMU_USHORT)
    {
      TsXYblocku2byte = alloc1d_us(0,nBlockRows*nCols-1);
      i0 = nyblock*nBlockRows*nCols;
      for (tvar = 0; tvar < nTimeSteps;tvar++)
      {
        getyrmonth(tvar, &yr, &mo);
        binfin = OpenFile(var, yr, mo, "r", currentfile);
        if (fseek(binfin, i0 * sizeof(short), SEEK_SET) != 0)
          nrerror("Error w/ fseek() while reading %s \n", currentfile);
        if (fread(TsXYblocku2byte, sizeof(short), nBlockRows*nCols, binfin) != (size_t) nBlockRows*nCols)
          nrerror("Could not do fread in %s \n", currentfile);
        fclose(binfin);
        for (i = 0; i < nBlockRows*nCols; i++)
          TsXYblock4byte[tvar][i] = (float) TsXYblocku2byte[i];
      }
      free1d_us(TsXYblocku2byte,0,nBlockRows*nCols-1);
    }
  }
  else if (istype(var.filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", var.filetype);
  else
    Warning("Data type %s currently unsupported\n", var.filetype);
}

/* ================================================================= */
/* XB Read and Write are designed for line-by-line processing of band data,
so it reads (or writes) the line of data specified in the argument list, and
places it (or reads from) into the array in the argument list.  The array
is assumed to be in [x][b] format. */

void ReadXBData(float * * ppfData, VAR * pVar, int y)
{
  FILE *fp;
/*  short int * psData;  */
  float * pfData;
  double * pdData;
  BYTE * pbyData;
  BOOL bSwap = FALSE;
/*  int i;  */
  char s[FILENAMELEN];
  char szFile[FILENAMELEN];
  int x, b;
  int nBands;


  nBands = pVar->nBands;

  if (istype(pVar->filetype, "nc"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else if (istype(pVar->filetype, "bin"))
  {                             /* READ BINARY FILES */
    bSwap = pVar->bByteSwap;

    if (pVar->tDataOrg == FSQ)
      Warning("Data org FSQ currently unsupported for input\n", "");
    else
    {
      if (GetFileName(pVar, szFile))
        Warning("Undefined Dynamic Alias in %s \n", szFile);

      if ((fp = fopen(szFile, "rb")) == NULL)
        Warning("Couldn't open file %s for reading\n", szFile);

      if (pVar->tDataOrg == BSQ)
      {
        if (pVar->tEmuType == EMU_FLOAT)
        {
          pfData = alloc1d_f(0,(nBands*nCols)-1);
          for (b=0; b < nBands; b++)
          {
            fseek(fp, (b * nRows * nCols + y * nCols) * sizeof(float), SEEK_SET);

            if (fread(&pfData[b * nCols], sizeof(float), nCols, fp) != (size_t)nCols)
              Warning("Could not do fread in %s \n", szFile);
          }

          if (bSwap)
            SwapBytes((BYTE *) pfData, 4, sizeof(float) * nBands * nCols);

          for (x=0; x < nCols; x++)
            for (b=0; b < nBands; b++)
              ppfData[x][b] = (float) pfData[x+b*nCols];

          free1d_f(pfData,0,nBands*nCols-1);
        }
        else if (pVar->tEmuType == EMU_BYTE)
        {
          pbyData = alloc1d_uc(0,(nBands*nCols)-1);
          for (b=0; b < nBands; b++)
          {
            fseek(fp, (b * nRows * nCols + y * nCols) * sizeof(BYTE), SEEK_SET);

            if (fread(&pbyData[b * nCols], sizeof(BYTE), nCols, fp) != (size_t)nCols)
              Warning("Could not do fread in %s \n", szFile);
          }

          for (x=0; x < nCols; x++)
            for (b=0; b < nBands; b++)
              ppfData[x][b] = (float) pbyData[x+b*nCols];

          free1d_uc(pbyData,0,nBands*nCols-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }

      }
      else if (pVar->tDataOrg == BIL)
        Warning("Data org BIL currently unsupported for input\n", "");
      else if (pVar->tDataOrg == BIP)
      {
        if (pVar->tEmuType == EMU_DOUBLE)
        {
          fseek(fp, y * nCols * nBands * sizeof(double), SEEK_SET);

          pdData = alloc1d_d(0,(nBands*nCols)-1);
          if (fread(pdData, sizeof(double), nCols*nBands, fp) != (size_t)(nCols*nBands))
            Warning("Could not do fread in %s \n", szFile);

          if (bSwap)  SwapBytes((BYTE *) pdData, 8, sizeof(double) * nBands * nCols);

          for (x=0; x < nCols; x++)
            for (b=0; b < nBands; b++)
              ppfData[x][b] = (float) pdData[x*nBands+b];

          free1d_d(pdData,0,nBands*nCols-1);
        }
        else if (pVar->tEmuType == EMU_FLOAT)
        {
          fseek(fp, y * nCols * nBands * sizeof(float), SEEK_SET);

          pfData = alloc1d_f(0,(nBands*nCols)-1);
          if (fread(pfData, sizeof(float), nCols*nBands, fp) != (size_t)(nCols*nBands))
            Warning("Could not do fread in %s \n", szFile);

          if (bSwap)
            SwapBytes((BYTE *) pfData, 4, sizeof(float) * nBands * nCols);

          for (x=0; x < nCols; x++)
            for (b=0; b < nBands; b++)
              ppfData[x][b] = (float) pfData[x*nBands+b];

          free1d_f(pfData,0,nBands*nCols-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }
      }
      fclose(fp);
    }
  }
  else if (istype(pVar->filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else
    Warning("Data type %s currently unsupported\n", pVar->filetype);
}

/* ================================================================= */
/* YXB Read and Write are designed for block processing of band data,
so it reads (or writes) the block of data specified in the argument list, and
places it (or reads from) into the array in the argument list.  The array
is assumed to be in [y][x][b] format. */

void ReadYXBData(float * * * pppfData, VAR * pVar, int yStart, int yExtent, int xStart, int xExtent)
{
  FILE *fp;
  float * pfData;
/*  double * pdData;
  BYTE * pbyData; */
  BOOL bSwap = FALSE;
  char s[FILENAMELEN];
  char szFile[FILENAMELEN];
  int y, x, b;
  int nBands;

  nBands = pVar->nBands;

  if (istype(pVar->filetype, "nc"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else if (istype(pVar->filetype, "bin"))
  {                             /* READ BINARY FILES */
    bSwap = pVar->bByteSwap;

    if (pVar->tDataOrg == FSQ)
      Warning("Data org FSQ currently unsupported for block input\n", "");
    else
    {
      if (GetFileName(pVar, szFile))
        Warning("Undefined Dynamic Alias in %s \n", szFile);

      if ((fp = fopen(szFile, "rb")) == NULL)
        Warning("Couldn't open file %s for reading\n", szFile);

      if (pVar->tDataOrg == BSQ)
        Warning("Data org BSQ currently unsupported for block input\n", "");
      else if (pVar->tDataOrg == BIL)
        Warning("Data org BIL currently unsupported for block input\n", "");
      else if (pVar->tDataOrg == BIP)
      {
        if (pVar->tEmuType == EMU_FLOAT)
        {

          pfData = alloc1d_f(0,(nBands*xExtent)-1);
          for (y=yStart; y<yStart+yExtent; y++)
          {
            fseek(fp, (y * nCols * nBands + xStart * nBands) * sizeof(float), SEEK_SET);

            if (fread(pfData, sizeof(float), xExtent*nBands, fp) != (size_t)(xExtent*nBands))
              Warning("Wrong Data size read in %s \n", szFile);

            if (bSwap)
              SwapBytes((BYTE *) pfData, 4, sizeof(float) * nBands * xExtent);

            for (x=0; x < xExtent; x++)
            {
              for (b=0; b < nBands; b++)
              {
                pppfData[y][x][b] = pfData[x*nBands+b];
              }
            }

          }

          free1d_f(pfData,0,(nBands*xExtent)-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }
      }
      fclose(fp);
    }
  }
  else if (istype(pVar->filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else
    Warning("Data type %s currently unsupported\n", pVar->filetype);
}

/* ================================================================= */
void WriteXY(float * pfData, VAR var, int nYear, int nMonth)
{
  FILE * fp;
  char szCurrentFile[FILENAMELEN];


  if (istype(var.filetype, "nc"))
  {   /* Netcdf */
    Warning("Output type %s currently unsupported\n", var.filetype);
  }
  else if (istype(var.filetype, "bin"))
  {       /* WRITE BINARY FILES */
    fp = OpenFile(var, nYear, nMonth, "wb", szCurrentFile);
    WriteBinaryOutput(fp, pfData, var.tEmuType, nRows*nCols);

    fclose(fp);
  }
  else if (istype(var.filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", var.filetype);
  else
    Warning("Data type %s currently unsupported\n", var.filetype);
}

/* ================================================================= */
void WriteXBData(float * * ppfData, VAR * pVar, int y)
{
  FILE *fp;
  short int * psData;
/*  double * pdData; */
  float * pfData;
  BYTE * pbyData;
  BOOL bSwap = FALSE;
/*  int i;  */
  char s[FILENAMELEN];
  char szFile[FILENAMELEN];
  int x, b;
  int nBands;
  int l_nCols, l_nRows;

  nBands = pVar->nBands;
  l_nCols = pVar->nCols;
  l_nRows = pVar->nRows;

  if (istype(pVar->filetype, "nc"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else if (istype(pVar->filetype, "bin"))
  {                             /* READ BINARY FILES */
    bSwap = pVar->bByteSwap;

    if (pVar->tDataOrg == FSQ)
    {
      for (b=0; b<nBands; b++)
      {
        if (GetFileName(pVar, szFile))
        {
          ReplaceDynAlias(szFile, "BAND", (float)b, EMU_INT);
          strcpy(s, szFile);
        }
        else
          sprintf(s, "%sb%d", szFile, b);  /* THIS IS A NAMING CONVENTION ! */

        if ((fp = fopen(s, "rb+")) == NULL)
        {
          if ((fp = fopen(s, "wb")) == NULL)
            Warning("Couldn't open file %s for writing\n", s);

          pbyData = alloc1d_uc(0,l_nCols-1);
          for (x=0; x < l_nRows; x++)
            fwrite(pbyData, EmuTypeSize(pVar->tEmuType), l_nCols, fp);
          free1d_uc(pbyData,0,l_nCols-1);
          fclose(fp);

          if ((fp = fopen(s, "rb+")) == NULL)
            Warning("Couldn't open file %s for writing\n", s);
        }
        if (pVar->tEmuType == EMU_BYTE)
        {
          pbyData = alloc1d_uc(0,l_nCols-1);

          for (x=0; x < l_nCols; x++)
            pbyData[x] = (BYTE) ppfData[x][b];

          fseek(fp, y * l_nCols, SEEK_SET);
          fwrite(pbyData, 1, l_nCols, fp);

          free1d_uc(pbyData,0,l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_SHORT)
        {
          psData = alloc1d_s(0,l_nCols-1);

          for (x=0; x < l_nCols; x++)
            psData[x] = (short) ppfData[x][b];

          fseek(fp, y * l_nCols * 2, SEEK_SET);
          fwrite(psData, 2, l_nCols, fp);

          free1d_s(psData,0,l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_FLOAT)
        {
          pfData = alloc1d_f(0,l_nCols-1);

          for (x=0; x < l_nCols; x++)
            pfData[x] = ppfData[x][b];

          fseek(fp, y * l_nCols * 4, SEEK_SET);
          fwrite(pfData, 4, l_nCols, fp);

          free1d_f(pfData,0,l_nCols-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }

        fclose(fp);
      }
    }
    else
    {
      if (GetFileName(pVar, szFile))
        Warning("Undefined Dynamic Alias in %s \n", szFile);

      if ((fp = fopen(szFile, "rb+")) == NULL)
      {
        if ((fp = fopen(szFile, "wb")) == NULL)
          Warning("Couldn't open file %s for writing\n", szFile);

        pbyData = alloc1d_uc(0,(nBands*l_nCols)-1);
        for (x=0; x < l_nRows; x++)
          fwrite(pbyData, EmuTypeSize(pVar->tEmuType), nBands*l_nCols, fp);
        free1d_uc(pbyData,0,nBands*l_nCols-1);
        fclose(fp);

        if ((fp = fopen(szFile, "rb+")) == NULL)
          Warning("Couldn't open file %s for writing\n", szFile);
      }
      if (pVar->tDataOrg == BSQ)
      {
        if (pVar->tEmuType == EMU_BYTE)
        {
          pbyData = alloc1d_uc(0,l_nCols-1);

          for (b=0; b < nBands; b++)
          {
            fseek(fp, b * l_nRows * l_nCols + y * l_nCols, SEEK_SET);
            for (x=0; x < l_nCols; x++)
              pbyData[x] = (BYTE) ppfData[x][b];
            fwrite(pbyData, 1, l_nCols, fp);
          }

          free1d_uc(pbyData,0,l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_SHORT)
        {
          psData = alloc1d_s(0,l_nCols-1);

          for (b=0; b < nBands; b++)
          {
            fseek(fp, (b * l_nRows * l_nCols + y * l_nCols)*sizeof(short), SEEK_SET);
            for (x=0; x < l_nCols; x++)
              psData[x] = (short) ppfData[x][b];
            fwrite(psData, sizeof(short), l_nCols, fp);
          }

          free1d_s(psData,0,l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_FLOAT)
        {
          pfData = alloc1d_f(0,l_nCols-1);

          for (b=0; b < nBands; b++)
          {
            fseek(fp, (b * l_nRows * l_nCols + y * l_nCols)*sizeof(float), SEEK_SET);
            for (x=0; x < l_nCols; x++)
               pfData[x] = ppfData[x][b];
            fwrite(pfData, sizeof(float), l_nCols, fp);
          }
          free1d_f(pfData,0,l_nCols-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }
      }
      else if (pVar->tDataOrg == BIL)
      {
        if (pVar->tEmuType == EMU_BYTE)
        {

         pbyData = alloc1d_uc(0,(nBands*l_nCols)-1);

          for (b=0; b < nBands; b++)
            for (x=0; x < l_nCols; x++)
               pbyData[b*l_nCols+x] = (BYTE) ppfData[x][b];

          fseek(fp, y * l_nCols * nBands, SEEK_SET);
          fwrite(pbyData, 1, nBands*l_nCols, fp);

          free1d_uc(pbyData,0,nBands*l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_SHORT)
        {
          psData = alloc1d_s(0,(nBands*l_nCols)-1);

          for (b=0; b < nBands; b++)
            for (x=0; x < l_nCols; x++)
               psData[b*l_nCols+x] = (short) ppfData[x][b];

          fseek(fp, y * l_nCols * nBands * 2, SEEK_SET);
          fwrite(psData, 2, nBands*l_nCols, fp);

          free1d_s(psData,0,nBands*l_nCols-1);
        }
        else if (pVar->tEmuType == EMU_FLOAT)
        {
          pfData = alloc1d_f(0,(nBands*l_nCols)-1);

          for (b=0; b < nBands; b++)
            for (x=0; x < l_nCols; x++)
               pfData[b*l_nCols+x] = ppfData[x][b];

          fseek(fp, y * l_nCols * nBands * 4, SEEK_SET);
          fwrite(pfData, 4, nBands*l_nCols, fp);

          free1d_f(pfData,0,nBands*l_nCols-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }
      }
      else if (pVar->tDataOrg == BIP)
        Warning("Data org BIP currently unsupported for output\n", "");

      fclose(fp);
    }
  }
  else if (istype(pVar->filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else
    Warning("Data type %s currently unsupported\n", pVar->filetype);
}

/* ================================================================= */
void WriteYXBData(float * * * pppfData, VAR * pVar, int yStart, int yExtent, int xStart, int xExtent)
{
  FILE *fp;
  float * pfData;
/*  double * pdData;
  BYTE * pbyData; */
  BOOL bSwap = FALSE;
  char s[FILENAMELEN];
  char szFile[FILENAMELEN];
  int y, x, b;
  int nBands;

  nBands = pVar->nBands;

  if (istype(pVar->filetype, "nc"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else if (istype(pVar->filetype, "bin"))
  {                             /* READ BINARY FILES */
    bSwap = pVar->bByteSwap;

    if (pVar->tDataOrg == FSQ)
      Warning("Data org FSQ currently unsupported for block output\n", "");
    else
    {
      if (GetFileName(pVar, szFile))
        Warning("Undefined Dynamic Alias in %s \n", szFile);

      if ((fp = fopen(szFile, "rb+")) == NULL)
        Warning("Couldn't open file %s for writing\n", szFile);

      if (pVar->tDataOrg == BSQ)
        Warning("Data org BSQ currently unsupported for block output\n", "");
      else if (pVar->tDataOrg == BIL)
        Warning("Data org BIL currently unsupported for block output\n", "");
      else if (pVar->tDataOrg == BIP)
      {
        if (pVar->tEmuType == EMU_FLOAT)
        {
          pfData = alloc1d_f(0,(nBands*xExtent)-1);
          for (y=yStart; y<yStart+yExtent; y++)
          {
            for (x=0; x < xExtent; x++)
            {
              for (b=0; b < nBands; b++)
              {
                pfData[x*nBands+b] = pppfData[y][x][b];
              }
            }
            if (bSwap)
              SwapBytes((BYTE *) pfData, 4, sizeof(float) * nBands * xExtent);

            fseek(fp, (y * nCols * nBands + xStart * nBands) * sizeof(float), SEEK_SET);
            if (fwrite(pfData, sizeof(float), xExtent*nBands, fp) != (size_t)(xExtent*nBands))
              Warning("Wrong Data size written in %s \n", szFile);
          }

          free1d_f(pfData,0,(nBands*xExtent)-1);
        }
        else
        {
          sprintf(s, "%d", pVar->tEmuType);
          Warning("Data type %s currently unsupported\n", s);
        }
      }
      fclose(fp);
    }
  }
  else if (istype(pVar->filetype, "ascii"))
    Warning("Data type %s currently unsupported\n", pVar->filetype);
  else
    Warning("Data type %s currently unsupported\n", pVar->filetype);
}

/* ================================================================= */
void WriteBMP(VAR * pVar, BYTE * * ppbyData)
{
  FILE * fp;
  int y;
  BITMAPFILEHEADER bmfh;
  BITMAPINFOHEADER bmih;
  int nBytesPerLine;
  char szFile[FILENAMELEN];

  if (GetFileName(pVar, szFile))
    Warning("Undefined Dynamic Alias in %s \n", szFile);

  if ((fp = fopen(szFile, "wb")) == NULL)
    Warning("Couldn't open file %s for writing\n", szFile);

  nBytesPerLine = pVar->nCols * 3;
  while (nBytesPerLine % 4)
    nBytesPerLine++;

  bmfh.bfType = 0x4d42;
  bmfh.bfSize = nBytesPerLine * pVar->nRows + 54;
  bmfh.bfRes1 = 0;
  bmfh.bfRes2 = 0;
  bmfh.bfOffBits = 54;

  bmih.biSize = sizeof(BITMAPINFOHEADER);
  bmih.biWidth = pVar->nCols;
  bmih.biHeight = pVar->nRows;
  bmih.biPlanes = 1;
  bmih.biBitCount = 24;
  bmih.biCompression = 0;
  bmih.biSizeImage = 0;
  bmih.biXPelsPerMeter = 0;
  bmih.biYPelsPerMeter = 0;
  bmih.biClrUsed = 0;
  bmih.biClrImportant = 0;

  fwrite(&bmfh.bfType, 2, 1, fp);
  fwrite(&bmfh.bfSize, 4, 1, fp);
  fwrite(&bmfh.bfRes1, 2, 1, fp);
  fwrite(&bmfh.bfRes2, 2, 1, fp);
  fwrite(&bmfh.bfOffBits, 4, 1, fp);
  fwrite(&bmih, sizeof(bmih), 1, fp);
  for (y=0; y<pVar->nRows; y++)
    fwrite((void *) ppbyData[y], nBytesPerLine, 1, fp);

  fclose(fp);
}

/* ================================================================= */
void WriteBinaryOutput(FILE * fp, float * pfData, emu_type tEmuType, int nElements)
{
  short int * psData;
  BYTE * pbyData;
  int i;
  char s[FILENAMELEN];

  if (tEmuType == EMU_FLOAT)
  {
    if (fwrite(pfData, sizeof(float), nElements, fp) != (size_t)(nElements))
      Warning("Couldn't write output data\n", "");
  }
  else if (tEmuType == EMU_SHORT)
  {
    psData = alloc1d_s(0,nElements-1);

    for (i=0; i < nElements; i++)
      psData[i] = (SHORT) pfData[i];

    if (fwrite(psData, sizeof(SHORT), nElements, fp) != (size_t)(nElements))
      Warning("Couldn't write output data\n", "");

    free1d_s(psData,0,nElements-1);
  }
  else if (tEmuType == EMU_BYTE)
  {
    pbyData = alloc1d_uc(0,nElements-1);

    for (i=0; i < nElements; i++)
      pbyData[i] = (BYTE) pfData[i];

    if (fwrite(pbyData, sizeof(BYTE), nElements, fp) != (size_t)(nElements))
      Warning("Couldn't write output data\n", "");

    free1d_uc(pbyData,0,nElements-1);
  }
  else
  {
    sprintf(s, "%d", tEmuType);
    Warning("Data type %s currently unsupported\n", s);
  }
}

/* ================================================================= */
/* Open output file for writing; return pointer to the output file  */

FILE *OUTPUTfile(char *filename, char *extension)
{
  FILE *fout;
  char writefile[FILENAMELEN];

  strcpy(writefile, pszOutputFilePath);
  strcat(writefile, "/");
  strcat(writefile, filename);
  /* strcat(writefile, "."); */
  strcat(writefile, extension);

  if ((fout = fopen(writefile, "w")) == NULL)
    nrerror("Couldn't open file %s for writing\n", writefile);

  return (fout);
}



/* ================================================================= */
/* ======================= LOCAL FUNCTIONS ========================= */

/* ================================================================= */

FILE * OpenFile(VAR var, int nYear, int nMonth, char * pszModeString, char * pszFile)
{
  FILE * fp;


  if (GetFileName(&var, pszFile))
  {
    ReplaceDynAlias(pszFile, "YEAR", (float)nYear, EMU_INT);
    ReplaceDynAlias(pszFile, "MONTH", (float)nMonth+1, EMU_INT);
  }
  else
    Warning("No Dynamic Alias in %s \n", pszFile);

  /* open file */
  if ((fp = fopen(pszFile, pszModeString)) == NULL)
    ExitError("Couldn't open file %s in OpenFile.\n", pszFile);

  return(fp);
}



/* ================================================================= */
/* ================================================================= */
/* ==== OBSOLESCENT FUNCTIONS, KEPT FOR BACKWARDS COMPATIBILITY ==== */

void ReadXYData(float * pfData, VAR * pVar)
{
  /* ReadBinaryData(fp, pfData, var.tEmuType, nRows*nCols, bSwap); */
  ReadXY(pfData, *pVar, NONE, NONE);
}
/* ================================================================= */
void WriteXYMap(float * pfData, VAR * pVar, int nYear, int nMonth)
{
  WriteXY(pfData, *pVar, nYear, nMonth);
}
/* ================================================================= */
void WriteXYData(float * pfData, VAR * pVar)
{
  WriteXY(pfData, *pVar, NONE, NONE);
}