/***************************************************************************/
/**                                                                       **/
/**             n  e  t  c  d  f  2  g  r  i  d  .  c                     **/
/**                                                                       **/
/**     Program converts a land mask in netcdf format to grid clm         **/
/**     file format.                                                      **/
/**     The land_mask fiel that is generated by the FMS grid generation   **/
/**     tools does not write axis values. It is hard assumed that the     **/
/**     coordinates run from 0 to 360 and -90 to 90.                      **/
/**     For LPJ grid file, longitudes are  converted to -180 to 180       **/
/**                                                                       **/
/**                                                                       **/
/**                                                                       **/
/**     written by Stefan Petri                                           **/
/**     inspired by code from Werner von Bloh                             **/
/**                                                                       **/
/**     Potsdam Institute for Climate Impact Research                     **/
/**     PO Box 60 12 03                                                   **/
/**     14412 Potsdam/Germany                                             **/
/**                                                                       **/
/**     Last change: $date$ **/
/**                                                                       **/
/***************************************************************************/

/* TODO: surround the land with a thin rim of ocean cells,
 * so that lateron the river runoff can always find a target cell
 * to finally send the water into the ocean.
 */

#include "lpj.h"
#include "netcdf.h"

static char *mask_variable_name = "mask";

static char *myname;

static void myusage(char *m) {
  if (m) fprintf(stderr, "%s", m);
  fprintf(stderr, "usage: %s netcdf-mask-file lpj-grid-file [mask-variable-name]\n"
          "reads a land mask file in netcdf format and creates an LPJ grid file\n",
          myname);
  exit(EXIT_FAILURE);
}

int main(int argc, char **argv)
{
  FILE *gridfile;
  Coord grid;
  int i, j, dim, err;
  Header header;
  
  int ncid, maskid, ndims, dimids[2];
  size_t dimlens[2];
  double *mask = NULL;
  char dimname[NC_MAX_NAME+1];
  double delta_lat, delta_lon; /* size of cells in degrees */
  double off_lat, off_lon;     /* offset of cell centers from edge in degrees */

  myname = argv[0];

  if (argc < 3) myusage("need two file names");
  if (argc == 4) mask_variable_name = argv[3];
  if (argc > 4) myusage("found too many arguments");

  err = nc_open(argv[1], NC_NOWRITE, &ncid);
  if (err) { fprintf(stderr, "nc_open(%s): %s\n", argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }

  err = nc_inq_varid(ncid, mask_variable_name, &maskid);
  if (err) { fprintf(stderr, "cannot find variable %s in %s: %s\n", mask_variable_name, argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }
  /* TODO: make sure that the data type is correct */

  err = nc_inq_varndims(ncid, maskid, &ndims);
  if (err) { fprintf(stderr, "cannot get number of dimensions for variable %s in %s: %s\n", mask_variable_name, argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }
  if (2 != ndims) { fprintf(stderr, "variable %s must have 2 dimensions but has %d\n", mask_variable_name, ndims); exit(EXIT_FAILURE); }
  err = nc_inq_vardimid(ncid, maskid, dimids);
  if (err) { fprintf(stderr, "cannot get dimension ids for variable %s in %s: %s\n", mask_variable_name, argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }

  for (dim = 0; dim < ndims; dim++) {
    err = nc_inq_dim(ncid, dimids[dim], dimname, &dimlens[dim]);
    if (err) { fprintf(stderr, "cannot get length of dimension %d for variable %s in %s: %s\n", dim, mask_variable_name, argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }
    printf("dimension %d name %s len %lu\n", dim, dimname, dimlens[dim]);
  }
  mask = (double *)malloc(dimlens[0] * dimlens[1] * sizeof(double));
  if (NULL == mask) { perror("cannot allocate storage for mask input"); exit(EXIT_FAILURE); }

  err = nc_get_var_double(ncid, maskid, mask);
  if (err) { fprintf(stderr, "cannot read data for variable %s from %s: %s\n", mask_variable_name, argv[1], nc_strerror(err)); exit(EXIT_FAILURE); }

  header.ncell=0;
  header.nbands=2;
  header.order=CELLYEAR;
  header.firstcell=0;
  header.firstyear=0;
  header.nyear=1;
  gridfile=fopen(argv[2],"wb");
  if(gridfile==NULL)
  {
    fprintf(stderr,"Error creating '%s': %s\n",argv[2],strerror(errno));
    myusage(NULL);
  }
  fwriteheader(gridfile,&header,LPJGRID_HEADER,LPJGRID_VERSION);

  delta_lat = (180.0 / dimlens[0]);
  delta_lon = (360.0 / dimlens[1]);
#if LPJGRID_VERSION >= 2
  header.scalar = 0.01;
#endif
#if LPJGRID_VERSION == 3
  header.datatype = LPJ_SHORT;
  header.cellsize_lon = delta_lon;
  header.cellsize_lat = delta_lat;
#else
  header.cellsize = delta_lon;
#endif

  off_lat = delta_lat / 2.0;
  off_lon = delta_lon / 2.0;
  printf("delta_lat %g off_lat %g delta_lon %g off_lon %g\n", delta_lat, off_lat, delta_lon, off_lon);
  for (j = 0; j < dimlens[0]; j++) {
    double lat = delta_lat * j + off_lat - 90.0;
    for (i = 0; i < dimlens[1]; i++) {
      double lon = delta_lon * i + off_lon;
      if (lon > 180) lon -= 360.0;
      if (mask[j*dimlens[1] + i] > 0.0) {
        grid.lon=lon;
        grid.lat=lat;
        writecoord(gridfile,&grid);
        header.ncell++;
      }
    }
  }
  printf("Number of grid cells: %d\n",header.ncell);

  /* write header again with correct number of grid cells */
  rewind(gridfile);
  fwriteheader(gridfile,&header,LPJGRID_HEADER,LPJGRID_VERSION);
  printf("header size %u LPJGRID_VERSION %d\n", sizeof header, LPJGRID_VERSION);
  fclose(gridfile);
  return EXIT_SUCCESS;
} /* of 'main' */