/* annmeannc.c
 *                *
 *              Emilio          3/31/99
 *                Updated program.
 *              Emilio          3/7/98
 *                Upgraded to new toolbox changes.
 *              Emilio          4/22/98
 *                *
 * COMPILE AS: gccemu annmeannc
 * RUN AS: annmeannc -o <Name of output NC file> -tit  -hist 
 ******************************************************************************/

#include "emu.h"

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

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

/* local variables */
float *x, *varf, **ansum;
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, "Annual Means");
  AssignString(&pszNChist, "annual means file created by annmeannc");

  Initialize(&argc, argv, DOMASK);

  Setup();
  ProcessAndOutput();

  CleanUp();

  return(0);
}

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

  ansum = alloc2d_f(0,nTotalYears-1,0,nRows*nCols-1);
  for (jyr = 0; jyr < nTotalYears; jyr++)
    for (i = 0; i < nRows*nCols; i++)
      ansum[jyr][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, FALSE);
  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 & WRITE TO ANNUAL MEANS ARRAY ansum === */
void ProcessAndOutput(void)
{
  int v, yr, mo, i, jyr, 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)
            ansum[yr-nFirstYear][i] += x[i]/12.0;
          }
    }  /* END tvar FOR-LOOP */

    Output(v);

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

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


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

  for (jyr = 0; jyr < nTotalYears; jyr++) {
    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) ansum[jyr][i])) ?
         pstVars[v].nFillValue : (short)floor((ansum[jyr][i]/pstVars[v].scale_factor) + 0.5);
      else if (pstVars[v].tEmuType == EMU_FLOAT)
        varf[iout] = (!masktest(pMask[i]) || IsBad((double) ansum[jyr][i])) ?
         pstVars[v].fFillValue : ansum[jyr][i];
      }  /* END i FOR-LOOP */

    /* write current time value to time-axis */
    start[0] = jyr;
    tvalue = (double) MONTH2HOUR*(YEAR2MONTH*((nFirstYear - YEARORIGIN) + jyr) - 5.999);
    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]; Year: yr=%d  ",v, nFirstYear+jyr);
    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 jyr FOR-LOOP */

}

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

/* ================================================================= */
void LocalCleanUp(void)
{
  free2d_f(ansum,0,nTotalYears-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);
}