/* vegdynvsci.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          12/7-8/98
 *                Wrote and debugged program.
 * COMPILE AS: gccemu vegdynvsci
 ******************************************************************************/

#include "emu.h"

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

/* local functions */
void Setup(void);
void ProcessAndOutput(void);

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


/* local variables */
int  binindex, bintot;
float binmin, binmax, bindel;
double *binmidpt;
float *x, *clindexMap, **out;
int errnc;


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Average Vegetation Dynamics vs. Climate Index");
  AssignString(&pszNChist, "file created by vegdynvsci");

  Initialize(&argc, argv, DOMASK);

  Setup();
  ProcessAndOutput();

  CleanUp();

  return(0);
}

/* ================================================================= */
void Setup(void)
{
  int b, tb, tvar;

  binmin = P[0];
  binmax = P[1];
  intot">bintot = (int)P[2];
  bindel = (binmax-binmin)/(float)bintot;

  /* write midpoint of each bin into array binmidpt */
  binmidpt = alloc1d_d(0,bintot-1);
  for (b=0; b < bintot; b++)
    binmidpt[b] = binmin + bindel*((double) b + 0.5);

  out = alloc2d_f(0,2-1,0,nTimeSteps*bintot-1);
  for (tvar=0; tvar < nTimeSteps; tvar++)
    for (b=0; b < bintot; b++)  {
      tb = bintot*tvar + b;
      out[SUM][tb]   = 0.0;
      out[COUNT][tb] = 0.0;
      }

  clindexMap = alloc1d_f(0,nRows*nCols-1);
  x = alloc1d_f(0,nRows*nCols-1);
}

/* ================================================================= */
void ProcessAndOutput(void)
{
  int v, yr, mo, i, b, tb, tvar, t1;    int bug=0;
  int ncID, tID;
  double *Tvalues;


  /* CREATE OUTPUT NETCDF FILE */
  Tvalues = alloc1d_d(0,nTimeSteps);
  t1 = (nFirstYear - (YEARORIGIN + 1))*12 + nFirstMonth;
  for (tvar=0; tvar < nTimeSteps; tvar++)
    Tvalues[tvar] = (double) MONTH2HOUR*(t1 + tvar + 0.001);

  FillAxisProp("z", "Rainfall Climate Index", "0-1", 1, binmin+bindel, binmax-bindel, bintot, binmidpt, FALSE);
  FillTAxisProp("time", "time", "default", 1, Tvalues[0], Tvalues[nTimeSteps-1], nTimeSteps, Tvalues, "default", FALSE);
  ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
  free1d_d(Tvalues, 0,nTimeSteps);


  /* read in climate index map (time-independent variable) */
  ReadXY(clindexMap, pstVars[0], NONE, NONE);

  for (v = 1; v < nVars; v++)
  {
    for (tvar = 0; tvar < nTimeSteps; tvar++)
    {
      getyrmonth(tvar, &yr, &mo);
      ReadXY(x, pstVars[v], yr, mo);

      for (b=0; b < bintot; b++)  {
        tb = bintot*tvar + b;
        out[SUM][tb]   = 0.0;
        out[COUNT][tb] = 0.0;
        }

      for (i = 0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          binindex = floor((double)(clindexMap[i] - binmin)/bindel);
          if (binindex >= 0 && binindex < bintot)
          {
            tb = bintot*tvar + binindex;
            out[SUM][tb]   += x[i];
            out[COUNT][tb] += 1.0;
          }
        }

      for (b=0; b < bintot; b++)  {
        tb = bintot*tvar + b;
        if (out[COUNT][tb] > 0)    out[MEAN][tb] /= out[COUNT][tb];
        }
    }

    /* Define and Write the output nc variable */
    writeNCvar(ncID, pstVars[v].param, "avci", "av long name", out[MEAN]);
    writeNCvar(ncID, pstVars[v].param, "count", "cnt long name", out[COUNT]);
  }

  errnc = nc_close(ncID);
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
}
/* ================================================================= */
void LocalCleanUp(void)
{
  free2d_f(out,0,2-1,0,nTimeSteps*bintot-1);
  free1d_d(binmidpt,0,bintot-1);
  free1d_f(clindexMap,0,nRows*nCols-1);
  free1d_f(x,0,nRows*nCols-1);
}

/* ============================================================ */
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, "zt", varLongName, "");
  errnc = nc_put_var_float(ncID, varID, array);
}