/* monmeannc.c
 *                *
 *              Emilio          3/27/98
 *                Upgraded to new toolbox changes.
 *                *
 * COMPILE AS: gcc -o monmeannc monmeannc.c -DUNIX -lm -lnetcdf -L./ -R./ -lemu
 * RUN AS: monmeannc -o <Name of output NC file> -tit  -hist 
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: monmeannc \n\
Using time series of maps, calculate and create a new netcdf\n\
dataset of monthly means (climatology).\n\
\n\
"};

/* local functions */
void Setup(void);
void ProcessAndOutput(void);
void Output(int v);

/* local variables */
float *x, *varf, **monsum;
short *vars;
int ncID;     /* netCDF file IDs  */
int varID;    /* variable ID    */
int tID;      /* time axis ID   */
static size_t count[3];
int errnc;


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Climatologies: Monthly Means");
  AssignString(&pszNChist, "monthly means file created by monmeannc");

  Initialize(&argc, argv, DOMASK);

  Setup();
  ProcessAndOutput();

  CleanUp();

  return(0);
}


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

  monsum = alloc2d_f(0,12-1,0,nRows*nCols-1);
  for (mo = 0; mo < 12; mo++)
    for (i = 0; i < nRows*nCols; i++)
      monsum[mo][i] = 0.0;

  x = alloc1d_f(0,nRows*nCols-1);
  varf = alloc1d_f(0,nRows*nCols-1);
  vars = alloc1d_s(0,nRows*nCols-1);

  /* setup new NC file from NC headers, define variables */
  SetupDefaultAxes(TRUE, TRUE);
  count[0] = 1;
  count[1] = stAxes[1].length;
  count[2] = stAxes[0].length;
  ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);
  DefNCvarsFromNCs(ncID, "none", 0);
}

/* ================================================================= */
/* ============= READ DATA, CALCULATE MONTHLY MEANS ================ */
void ProcessAndOutput(void)
{
  int v, yr, mo, i, tvar = 0;     int bug=0;


  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);

      for (i = 0; i < nRows*nCols; i++)
        if (masktest(pMask[i]))
        {
          /* catch floating-point underflow errors with tolerance "TOL" */
          if (x[i] > TOL || x[i] < -TOL)
            monsum[mo][i] += x[i]/(float)nTotalMonths[mo];
        }
    } /* END tvar FOR-LOOP */

    Output(v);

    /* clear monsum array for current month mo, to set up for the next variable */
    for (mo = 0; mo < 12; mo++)
      for (i = 0; i < nRows*nCols; i++)
        monsum[mo][i] = 0.0;
  } /* END v FOR-LOOP */
}

/* ================================================================= */
/* ===============  WRITE DATA TO FILE(S)  ========================= */
void Output(int v)
{
  int mo, i, iout, col, row;
  double tvalue = 0;
  static size_t start[] = {0,0,0};


  errnc = nc_inq_varid(ncID, pstVars[v].param, &varID);

  for (mo = 0; mo < 12; mo++)
  {
    for (i = 0; i < nRows*nCols; i++)
    {
      iout = getRowRevCol(i, &col, &row);

      /* if outside pMask or IsBad, use _FillValue */
      /*
      if (pstVars[v].tEmuType == EMU_SHORT)
        vars[iout] = (!masktest(pMask[i]) || IsBad((double) monsum[mo][i])) ?
         pstVars[v].nFillValue : (short)floor((monsum[mo][i]/pstVars[v].scale_factor) + 0.5);
      else if (pstVars[v].tEmuType == EMU_FLOAT)
        varf[iout] = (!masktest(pMask[i]) || IsBad((double) monsum[mo][i])) ?
         pstVars[v].fFillValue : monsum[mo][i];
      */
      /* THERE IS A PROBLEM WITH THE USE OF IsBad(), BUT I DON'T UNDERSTAND IT */


      if (pstVars[v].tEmuType == EMU_SHORT)
        vars[iout] = (!masktest(pMask[i])) ?
         pstVars[v].nFillValue : (short)floor((monsum[mo][i]/pstVars[v].scale_factor) + 0.5);
      else if (pstVars[v].tEmuType == EMU_FLOAT)
        varf[iout] = (!masktest(pMask[i])) ?
         pstVars[v].fFillValue : monsum[mo][i];
    }  /* END i FOR-LOOP */

    /* write current time value to time-axis */
    start[0] = mo;
    tvalue = (double) MONTH2HOUR*(mo + 0.501);
    errnc = nc_put_var1_double(ncID, tID, &start[0], &tvalue);

    /* write array for this month into netCDF file */
    if (bDebug)  printf("\n   Write variable pstVars[%d]; Month: mo=%d  ",v, mo);
    if (pstVars[v].tEmuType == EMU_SHORT)
      errnc = nc_put_vara_short(ncID, varID, start, count, vars);
    else if (pstVars[v].tEmuType == EMU_FLOAT)
      errnc = nc_put_vara_float(ncID, varID, start, count, varf);
  } /* END mo FOR-LOOP */

}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{}
/* ================================================================= */
void LocalCleanUp(void)
{
  free2d_f(monsum,0,12-1,0,nRows*nCols-1);
  free1d_f(x,0,nRows*nCols-1);
  free1d_f(varf,0,nRows*nCols-1);
  free1d_s(vars,0,nRows*nCols-1);

  errnc = nc_close(ncID);
}