/* climateindex.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          12/3/98
 *                Wrote program. Debugged completely.
 * COMPILE AS: gccemu climateindex
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: climateindex \n\
\n\
.\n\
"};

/* local functions */
void Setup(void);
void Process(void);
void Output(void);

void writeBinaryVar(char *filename, char *Type, float *array);
void writeNCvar(int ncID, char *parameter, char *Type, char *LongType, float *array);

/* local variables */
float **TXYb_1, **TXYb0, **TXYb1, ***out;
int errnc;

#define  NUMPARAM  3
#define  CLINDX    0
#define  CLINDXMO  1
#define  SUMYR     0
#define  SUM3MO    1


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Climate Index: 3-month minimum / annual sum");
  AssignString(&pszNChist, "Climate index file created by climateindex");

  Initialize(&argc, argv, DOMASK);

  Setup();
  Process();
  Output();

  CleanUp();

  return(0);
}

/* ================================================================= */
void Setup(void)
{
  int v, i;

  /* allocate memory and initialize out[][][] TO 0.0  */
  out = alloc3d_f(0,2-1,0,nVars-1,0,nRows*nCols-1);
  for (v=0; v < nVars;v++)
    for (i=0; i < nRows*nCols;i++)
    {
      out[CLINDX][v][i]   = 0.0;    /* climate index                 */
      out[CLINDXMO][v][i] = 0.0;    /* center month of climate index */
    }

  TXYb_1 = alloc2d_f(0,12-1,0,nBlockRows*nCols-1);
  TXYb0 = alloc2d_f(0,12-1,0,nBlockRows*nCols-1);
  TXYb1 = alloc2d_f(0,12-1,0,nBlockRows*nCols-1);
}

/* ================================================================= */
/* ======== READ DATA AND CALCULATE CLIMATE INDEX ========= */
void Process(void)
{
  int v, i, iout, col, row, i0, nyblock, mo, mo_min;
  float txyval_1, txyval0, txyval1, fill_val;
  double stats_i[NUMPARAM];


  for (v = 0; v < nVars; v++)
  {
    if (!bQuiet)  printf("v=%d *** \n",v);
    for (nyblock = 0; nyblock < nYBlocks; nyblock++)
    {
      ReadTXYbl(TXYb_1, pstVars[v], nyblock);
      ReadTXYbl(TXYb0, pstVars[v], nyblock);
      ReadTXYbl(TXYb1, pstVars[v], nyblock);
      i0 = nyblock*nBlockRows*nCols;

      for (i = i0; i < i0+nBlockRows*nCols; i++)
      {
        iout = (bOutBinary) ? i : getRowRevCol(i, &col, &row);
        if (masktest(pMask[i]))
        {
          stats_i[SUMYR]  = 0.0;
          stats_i[SUM3MO] = 0.0;
          stats_i[MIN]    = 100000000.0;

          for (mo = 0; mo < 12;mo++)
          {
            txyval_1 = (mo == 0) ? TXYb_1[11][i-i0]:TXYb_1[mo][i-i0];
            txyval0  = TXYb0[mo][i-i0];
            txyval1  = (mo == 11) ? TXYb1[0][i-i0]:TXYb1[mo][i-i0];

            /* catch floating-point underflow errors w/ tolerance "TOL" */
            /* if (txyval > TOL || txyval < -TOL) */

            stats_i[SUMYR] += txyval0;
            stats_i[SUM3MO] = txyval_1 + txyval0 + txyval1;

            if (stats_i[SUM3MO] < stats_i[MIN])
            {
              stats_i[MIN] = stats_i[SUM3MO];
              mo_min = mo+1;
            }
          } /* END mo FOR-LOOP (loop thru each time step) */

          out[CLINDX][v][iout]   = stats_i[MIN] / stats_i[SUMYR];
          out[CLINDXMO][v][iout] = mo_min;

        } /* END masktest() IF-BLOCK */
        else
        {                /* outside the pMask: FillValue */
          if (pstVars[v].tEmuType == EMU_SHORT)
            fill_val = (float) pstVars[v].nFillValue;
          else if (pstVars[v].tEmuType == EMU_FLOAT)
            fill_val = pstVars[v].fFillValue;

          out[CLINDX][v][iout]   = fill_val;
          out[CLINDXMO][v][iout] = fill_val;
        }
      } /* END i-loop (read each pixel within nyblock) */
    } /* END nyblock-loop (loop thru each nyblock) */
  } /* END v FOR-LOOP (loop thru each variable) */
}

/* ================================================================= */
/* =============  WRITE CLIMATE INDEX TO FILE(S)  ================== */
void Output(void)
{
  int v, ncID, tID;
  char FileName[NAMELEN];         /* Base File Name for binary files output */

  if (bOutBinary)    /* write output to binary files */
    for (v = 0; v < nVars; v++)  {
      sprintf(FileName,"%s.%s", pszOutputFile, pstVars[v].param);
      writeBinaryVar(FileName, ".minindx", out[CLINDX][v]);
      writeBinaryVar(FileName, ".momin", out[CLINDXMO][v]);
      }
  else  {                              /* write output to a single netCDF file */
    /* setup new NC file from NC headers */
    SetupDefaultAxes(FALSE, FALSE);
    ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
    for (v = 0; v < nVars; v++)   {
      writeNCvar(ncID, pstVars[v].param, "minindx", "3-Month Min Index", out[CLINDX][v]);
      writeNCvar(ncID, pstVars[v].param, "momin", "Center month of 3-Month Min", out[CLINDXMO][v]);
    }
    errnc = nc_close(ncID);
    }
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(out,0,2-1,0,nVars-1,0,nRows*nCols-1);
  free2d_f(TXYb_1,0,12-1,0,nBlockRows*nCols-1);
  free2d_f(TXYb0,0,12-1,0,nBlockRows*nCols-1);
  free2d_f(TXYb1,0,12-1,0,nBlockRows*nCols-1);
}

/* ============================================================ */
void writeBinaryVar(char *filename, char *Type, float *array)
{
  FILE *fout;

  fout = OUTPUTfile(filename, Type);
  fwrite(array, sizeof(float), nRows*nCols, fout);
  fclose(fout);
}

/* ============================================================ */
void writeNCvar(int ncID, char *parameter, char *Type, char *LongType, float *array)
{
  int varID;                  /* variable ID    */
  char varName[NAMELEN];      /* Name of NC variable(s) to be written      */
  char varLongName[NAMELEN];  /* Long Name of NC variable(s) to be written */

  sprintf(varName,"%s_%s", parameter,Type);
  sprintf(varLongName,"%s %s", parameter,LongType);
  varID = DefNCvar(ncID, varName, EMU_FLOAT, "xy", varLongName, "0-1");
  errnc = nc_put_var_float(ncID, varID, array);
}