#ifndef GLM_GRIDDED_INPUT_H_
#define GLM_GRIDDED_INPUT_H_

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <vector>
#include "netcdf.h"
#include "config.h"
#include "grid_info.h"

namespace glm {

inline
void check_err (const int stat, const int line, const char *file) {
   if (stat != NC_NOERR) {
      (void) fprintf(stderr, "line %d of %s: %s\n", line, file,
                     nc_strerror(stat));
      exit(1);
   }
}

#if 1
template<class T>
void aggregate (T orig[NY][NX], T agg[AGG_NY][AGG_NX], GridInfo &grid) {
   memset(agg, 0, sizeof(T)*AGG_NY*AGG_NX);
   int ratio = NX/AGG_NX;

   for (size_t n=0; n<NY; n++) {
      for (size_t m=0; m<NX; m++) {
         
         size_t i, j;
         
         if ((n+1) % ratio == 0) {
            i = (n+1)/ratio - 1;
         } else {
            i = (n+1 + (ratio - ( (n+1) % ratio ) ) ) / ratio - 1;
         }
         
         if ((m+1) % ratio == 0) {
            j = (m+1)/ratio - 1;
         } else {
            j = (m+1 + (ratio - ( (m+1) % ratio ) ) ) / ratio - 1;
         }
         
         if (orig[n][m] > 0 ){
            agg[i][j] = agg[i][j] + orig[n][m]* grid.cellArea[n][m] / grid.cellArea2Deg[i][j];
         }
         
      } // end of m-loop
   }// end of n-loop
}
#endif

template <class T>
class GriddedInput {
protected:
   int ncID_;
   int varID_;
   int firstYear_;

   void readVar (size_t start[], size_t count[], T* buffer);
   void readVar (int ncid, int varid, size_t start[], size_t count[], T* buffer);

public:
   T buffer_[NSTEP][NY][NX];

   GriddedInput (const char* fileName, const char* varName) {
      // Reads in netCDF file, gets first year from its time dimension
      // Checks for any errors in reading time dimension
      // Checks for any errors in reading varName
      int status = nc_open(fileName, NC_NOWRITE, &ncID_);
      check_err(status,__LINE__,__FILE__);

      int timeID;
      // nc_inq_varid returns the ID of a netCDF variable given its name. Its parameters are:
      // ncid: NetCDF ID, from a previous call to nc_open or nc_create
      // name: Variable name for which ID is desired
      // varidp: Pointer to location for returned variable ID
      status = nc_inq_varid(ncID_, "time", &timeID);
      check_err(status,__LINE__,__FILE__);

      size_t index[] = {0};
      // nc_get_var1_<type> get a single data value from a variable of an open netCDF dataset
      // ncid: NetCDF ID, from a previous call to nc_open or nc_create
      // varid: Variable ID
      // index[]: The index of the data value to be read
      // In this case, firstYear_ stores the value of the first year read in from the netCDF file
      status = nc_get_var1_int(ncID_, timeID, index, &firstYear_);
      check_err(status,__LINE__,__FILE__);

      status = nc_inq_varid(ncID_, varName, &varID_);
      check_err(status,__LINE__,__FILE__);
   }

   virtual void load () {
      load(firstYear_);
   }

   virtual void load (int year) {
      size_t startYear = year - firstYear_;
      size_t start[] = {startYear, 0, 0};
      size_t count[] = {1, NY, NX};
      
      readVar(start, count, &buffer_[0][0][0]);

   }

   virtual void load (const char* initFileName, const char* initVarName,
                      int firstIAMYear) {

   }
   
   virtual void postLoad () {
      // do nothing. only needed for harmonization
   }

   virtual T get(size_t y, size_t x) {
      return buffer_[0][y][x];
   }
   
   friend class InputFiles;

};

template <class T>
class GriddedHarmonizedInput : public GriddedInput<T> {
   size_t thisYear_;
   GridInfo* grid_;
   int firstYear_;
   
public:
   // This will store current and next timestep from IAM
   // Currently assuming same resolution as GLM, but it wouldn't have to be (mostly)
   T iamBuffer_[2][NY][NX];
   // This will store those two layers aggregated to 2.0 degree (or some such)
   // as well as the current state of GLM aggreaged to same resolution
   T aggBuffer_[3][AGG_NY][AGG_NX];
   // This will store the harmonized 2.0 degree data before it is disagregated and interpolated
   T harmBuffer_[2][AGG_NY][AGG_NX];

   size_t now_;
   size_t next_;

   GriddedHarmonizedInput (const char* fileName, const char* varName,
                           GridInfo* grid)
   : GriddedInput<T> (fileName, varName) {
      thisYear_ = NSTEP-1; // this will make sure data is loaded the first year
      now_ = 0;
      next_ = 1;
      grid_ = grid;
      
   }
   
   void load (const char* initFileName, const char* initVarName, int firstIAMYear) {
     // This is the first load. Put the last year from the historic data into the last
      // layer of the buffer.
      int ncID, varID, timeID, status;
      size_t timeLength;
      
      status = nc_open(initFileName, NC_NOWRITE, &ncID);
      check_err(status,__LINE__,__FILE__);
      
      // Get the time dimension and determine its length. We want to load the last record
      status = nc_inq_dimid (ncID, "time", &timeID);
      check_err(status,__LINE__,__FILE__);
      
      status = nc_inq_dimlen(ncID, timeID, &timeLength);
      check_err(status,__LINE__,__FILE__);

      status = nc_inq_varid (ncID, initVarName, &varID);
      check_err(status,__LINE__,__FILE__);
      
      size_t index[] = {0};
      status = nc_get_var1_int(ncID, timeID, index, &firstYear_);
      check_err(status,__LINE__,__FILE__);
      
      size_t start[] = {firstIAMYear-firstYear_, 0, 0};
     
      size_t count[] = {1, NY, NX};

      GriddedInput<T>::readVar(ncID, varID, start, count, &GriddedInput<T>::buffer_[NSTEP-1][0][0]);

      ncclose(ncID);
      
      // Then load the first year of the IAM data into the iamBuffer "now"
      start[0] = static_cast<size_t>(firstIAMYear - GriddedInput<T>::firstYear_);
      // For now, assume count has same dimensions for iam data
      
      GriddedInput<T>::readVar(start, count, &iamBuffer_[now_][0][0]);
      // Aggregate it to 2.0 degrees
      aggregate(iamBuffer_[now_], aggBuffer_[now_], *grid_);
   }
   
   void load (int year) {
      // increment thisYear_
      thisYear_++;
      
      if (thisYear_ == NSTEP) { // we need to load more data
         // aggregate the current state (will be in the last time step of 'buffer_')
         // aggregated results will be placed into the last slot of aggBuffer_ (first two slots
         // are for 'now_' and 'next_' of the IAM data
         aggregate(GriddedInput<T>::buffer_[NSTEP-1], aggBuffer_[2], *grid_);

         // Load the next IAM year into "next" iamBuffer
         size_t start[] = {static_cast<size_t>(year - GriddedInput<T>::firstYear_), 0, 0};
         size_t count[] = {1, NY, NX};
  
         GriddedInput<T>::readVar(start, count, &iamBuffer_[next_][0][0]);
	
         // aggregate the incoming IAM data to 2.0 degree
         aggregate(iamBuffer_[next_], aggBuffer_[next_], *grid_);

      }
   }
   
   void postLoad() {
      // swap 'now_' and 'next_' in preperation for next load call
      next_ = now_;
      now_ = (now_ == 0) ? 1 : 0;
      
      thisYear_ = 0;
   }
   
   T get(size_t y, size_t x) {
      return GriddedInput<T>::buffer_[thisYear_][y][x];
   }
   
   
   

}; // END GriddedHarmonizedInput


} // END glm namespace
#endif // GLM_GRIDDED_INPUT_H_