/* runoffbasins.c
 *                *
 *              Emilio          4/1/99
 *                Updated program.
 *              Emilio          7/21/98
 *                *
 * COMPILE AS: gccemu runoffbasins
 * RUN AS: runoffbasins <output file>
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: runoffbasins \n\
Given a mask file with unique integer ids (values), runoffbasins\n\
creates a table containing time series of total runoff in m^3/sec.\n\
\n\
Command-line switches specific to this program are:\n\
\t-nh    Don't print series headers (1st two rows) \n\
\n\
"};

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

/* local variables */
float *x, ***series;
double *sum_mo, mean_mo;
BOOL bNoHeaders;
int errnc;


int main(int argc, char *argv[])
{
  Initialize(&argc, argv, DOMASK);

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

  CleanUp();

  return(0);
}

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

  series = alloc3d_f(0,FlMskSitesN-1,0,nVars-1,0,nTimeSteps-1);
  for (sindex=0; sindex < FlMskSitesN; sindex++)
    for (v=0; v < nVars;v++)
      for (i=0; i < nTimeSteps;i++)
        series[sindex][v][i] = 0.0;

  x = alloc1d_f(0,nRows*nCols-1);
  sum_mo = alloc1d_d(0,FlMskSitesN-1);
}

/* ================================================================= */
void Process(void)
{
  int v, yr, mo, i;    int bug=0;
  int tvar = 0, sindex, n;

  for (v=0; v < nVars; v++)
  {
    for (tvar=0; tvar < nTimeSteps; tvar++)
    {
      getyrmonth(tvar, &yr, &mo);
      if (bDebug)   printf("\n pstVars[%d]=%s: tvar=%d, year=%d, month=%d  *  ",
                           v,pstVars[v].param,tvar,yr,mo+1);
      ReadXY(x, pstVars[v], yr, mo);

      /* initialize/clear cumulative temp variables */
      mean_mo = 0.0;
      for (sindex=0; sindex < FlMskSitesN; sindex++)
        sum_mo[sindex]  = 0.0;

      for (i=0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          sindex = siteindex(FlMskSitesIDs, pMask[i]);

          /* catch floating-point underflow errors with tolerance "TOL" */
          if (x[i] > TOL || x[i] < -TOL)
            sum_mo[sindex] += cellarea(i) * x[i];
        }

      for (n=0; n < FlMskSitesN; n++)
      {
        mean_mo = sum_mo[n] / (double) FlMskSitesArea[n];
        series[n][v][tvar] = mean_mo;
      }
    }   /* END tvar FOR-LOOP (loop thru each time step) */
  }  /* END v FOR-LOOP (loop thru each variable) */
}

/* ================================================================= */
/* ========== WRITE OUT COLUMN (SERIES) HEADERS PER SITE =========== */
void Output(void)
{
  FILE  *fout;
  double UnitConversion;
  int yr, mo, n, tvar = 0;    int bug=0;
  int m, sindex;


  fout = OUTPUTfile(pszOutputFile, ".mts");

  /* write series headers */
  if (!bNoHeaders)
    fprintf(fout,"RUNOFF\n");

  for (n=0; n < MskSitesN; n++)
  {
    if (!bQuiet) printf("\n  sites: MskSitesIDs[%d] = %d, .cellcount = %d, .area = %.1f ",
                        n,MskSitesIDs[n],Masks[n].cellcount,Masks[n].area/1000000.0);
    /* fprintf(fout,", %d", MskSitesIDs[n]); */
    fprintf(fout,"%d,", MskSitesIDs[n]);
  }


  /* write series data (series[][][]) per site  */
  for (tvar=0; tvar < nTimeSteps; tvar++)
  {
    getyrmonth(tvar, &yr, &mo);

    if (!bNoHeaders)
      if (pstVars[0].tTimeType == MONMEAN)       /* if monmean, don't print series a year */
        fprintf(fout,"\n%d", mo+1);
      else
        fprintf(fout,"\n%d/%d", mo+1, yr);

    fprintf(fout,"\n");

    for (n=0; n < MskSitesN; n++)
    {
        mean_mo = 0.0;
        for (m=0; m < Masks[n].nMaskids; m++)
        {
          sindex = siteindex(FlMskSitesIDs, Masks[n].maskids[m]);
          /* cellcountWt = (double) FlMskSitesCellcount[sindex]  /
                           (double) Masks[n].cellcount; */
          UnitConversion = FlMskSitesArea[sindex] * MM2M / MONTH2SEC;
          mean_mo  += series[sindex][0][tvar] * UnitConversion;
          /* runoff = (runofflow+runoffup [in mm])*(10^-3 m/mm) * (total area of pixels)
                    * (10^6 m2/km2) / [(30.417 days/month)*(86400 secs/day)] */
        }
        if (!bNoHeaders)
          fprintf(fout,", %f", mean_mo);
        else
          fprintf(fout,"%f,", mean_mo);
    }
  } /* END tvar FOR-LOOP */

  fclose(fout);
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
  int i;

  for (i=1; i<*pArgc; i++)
    if (strcmp(argv[i], "-nh") == 0)
      bNoHeaders = TRUE;
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free3d_f(series,0,FlMskSitesN-1,0,nVars-1,0,nTimeSteps-1);
  free1d_f(x,0,nRows*nCols-1);
  free1d_d(sum_mo,0,MskSitesN-1);
}