/**
 * \file drainage2ocean.c
 *
 * $Id$
 *
 * Reads an LPJ grid.bin file and a drainage map in ASCII format, e.g. DDM30.asc.
 * Produces a NetCDF file with information about river runoff routing into the ocean.
 *
 * Augmentation of the LPJ utility drainage.c .
 *
 * The data format of the original LPJ drainage maps, as generated by the utility
 * drainage.c, allows for storing river routing information only for runoff the goes
 * from a land cell to another land cell.
 * That drainage.c converts a drainage map in ASCII format (DDM30.asc)
 * together with a LPJ grid specification into a binary LPJ drainage map.
 * Since the LPJ is not at all interested in ocean cells,
 * there is no information in that binary map file where water flows out
 * into the ocean, even though that information is contained in the original
 * ASCII format map.
 * 
 * This program, in contrast, retains the information about runoff into the ocean.
 * It generates a NetCDF map for all cells land and ocean, and stores for each (land)
 * cell the destination coordinates of the runoff from that cell, regardless if the
 * destination is land or ocean.
 * To ease the task of the LPJ-FMS interface, also a mask of those coast cells is
 * computed, which actually drain into an ocean cell.
 *
 * NOTE: DDM30.asc data spans -180 to 180 deg, also LPJ uses that longitude range;
 *       but the land_lad grid, as defined by the FMS coupler, is set up as 0 to 360 deg.
 * Thus, the runoff information is encoded in lon/lat coordinates,
 * not in cell indices.
 *
 * The output file is organized in a global regular grid, so that it is easily compatible with
 * the FMS domain decomposition.
 *
 * This program also computes the drainage directions encoded as vector components in
 * x,y directions. This can be used e.g. with ferret to create overlay plots.
 *
 * Encoding of the original DDM30.asc :
 * https://www.uni-frankfurt.de/45217896/3_drainage_direction_map
 * The map is generated with a cellsize of 30 arc minutes and is
 * distributed as a zipped ASCII-grid that can be easily imported in
 * most GIS-software that support rasters or grids. The direction of
 * outflow from cell is according to the ArcInfo/ArcView convention as
 * follows: "0" is assigned to internal sinks and to cells draining
 * into the ocean while "1", "2", "4", "8", "16", "32", "64" and "128"
 * standing for the outflow directions E, SE, S, SW, W, NW, N or NE.  
 * -9 : missing value / ocean cell
 * -88: DDM30 ocean cell inside LPJ land area, as result from an earlier
 *     processing step with the utility asc2drainage2asc .
 *
 * Encoding of the variables in the generated NetCDF file:
 * dest_lon: longitude of destination cell; -180 <= dest_lon <= 180;
 * dest_lat: latitude of destination cell; -90 <= dest_lat <= 90;
 *    -999 = missing_value = ocean cell
 *    -333 = land cell that is ocean in the input routing file
 *    sink cells without outflow route into themselves
 * drain2ocean: mask of cells from which the runoff must be
 *    taken and forwarded to the ocean component of FMS
 *      1 for land (coast) cells which actually drain into ocean cells.
 *      0 == missing_data else
 *
 * written by Stefan Petri
 * inspired by code from Werner von Bloh and Stefanie Rost
 * Potsdam Institute for Climate Impact Research
 * PO Box 60 12 03
 * 14412 Potsdam/Germany
 */

#include "lpj.h"

#include <netcdf.h>
#include <time.h>
#include <assert.h>

static char *myname;

static void myusage(char *m) {
  if (m) fprintf(stderr, "%s\n", m);
  fprintf(stderr, "usage: %s lpj_grid_file DDM30.asc netcdf_file\n"
          "  reads an LPJ grid description and a drainage map in ASCII format, e.g. DDM30.asc,\n"
          "  and writes a NetCDF file with information about river routing into the ocean.\n"
          "  (CF-1.0, regular lat-lon grid).\n"
          "  dest_lon: longitude of destination cell; -180 <= dest_lon <= 180;\n"
          "  dest_lat: latitude of destination cell; -90 <= dest_lat <= 90;\n"
          "    -999 = missing_value = ocean cell;\n"
          "    -333 = land cell that is ocean in the input routing file;\n"
          "  sink cells without outflow route into themselves.\n"
          "  drain2ocean: 1 for coast cells wich actually drain into an ocean cell;\n"
          "               0 == missing_value else.\n"
          "  The program also computes the drainage directions encoded as vector\n"
          "  components in x/y directions.\n"
          "  vectors give x/y directions as pairs of -1/0/1 values, suitable for plotting with ferret:\n",
          myname);
  fprintf(stderr,
          "let xx=x[gx=dirx]; let xmax=`xx[x=@max]`; let xmin=`xx[x=@min]`;\n"
          "let yy=y[gy=dirx]; let ymax=`yy[y=@max]`; let ymin=`yy[y=@min]`;\n"
          "let ii=i[gx=dirx]; let nlons=`ii[i=@max]`; let imin=`ii[i=@min]`; let imax=`ii[i=@max]`\n"
          "let jj=j[gy=dirx]; let nlats=`jj[j=@max]`; let jmin=`jj[j=@min]`; let jmax=`jj[j=@max]`\n"
          "let resx=`(xmax-xmin)/nlons`; let resy=`(ymax-ymin)/nlats`\n"
          "let dirxres=dirx*resx; let diryres=diry*resy\n"
          "set region/x=15:40/y=20:35; go basemap; go land; can region\n"
          "cancel mode verify\n"
          "repeat/j=`jmin`:`jmax` ( \\\n"
          "  repeat/i=`imin`:`imax` ( \\\n"
          "    if `dirx gt -2` then ( \\\n"
          "      say `i`,`j` {`xx`,`xx+dirxres`},{`yy`,`yy+diryres`} ; \\\n"
          "      plot/vs/over/color=black/line/nolabel/nokey {`xx`,`xx+dirxres`},{`yy`,`yy+diryres`} \\\n"
          "    )\\\n"
          "  )\\\n"
          ")\n");
  exit(EXIT_FAILURE);
}

#define ncerror(e) { if(e) {fprintf(stderr,"Error: %s\n", nc_strerror(e)); exit(1);}}

#define OCEAN_CELL -999


typedef struct
{
  const Coord_array *index;
  int ncid;
  int vardest_iid, vardest_jid, vardirxid, vardiryid;  // i/j index of dest cell var, x-component var, y-component var
  int vardrain2ocean_id;
  //int varroutexid, varrouteyid; // x and y components of route markers
  //int dcount_var_id; // count of cells which drain (indirectly) into this one
  //int basin_var_id;
  //int basinboundx_var_id, basinboundy_var_id; // basin boundaries
} Cdf;


static int ///< returns NetCDF id of created file
create_cdf_file(const char *name
#ifdef USE_NETCDF4
                , int compress
#endif
                )
{
  int rc;
  int ncid = -1;
#ifdef USE_NETCDF4
  rc = nc_create(name, compress ? NC_CLOBBER|NC_NETCDF4 : NC_CLOBBER, &ncid);
#else
  rc = nc_create(name, NC_CLOBBER, &ncid);
#endif
  if (rc) {
    fprintf(stderr,"ERROR426: Cannot create file '%s': %s.\n",
            name,nc_strerror(rc));
    exit(-1);
  }
  return ncid;
}

static int ///< returns NetCDF id of created dimension
create_cdf_dim(int ncid, ///< file
               const char *name,
               int len) {
  int rc;
  int dimid = -1;
  rc = nc_def_dim(ncid, name, len, &dimid);
  ncerror(rc);
  return dimid;
}

static int ///< returns NetCDF id of created variable
create_cdf_var(int ncid, ///< file
               const char *name,
               const char *units,
               const char *long_name,
               int nc_type,
               const void *fill_value,
               int ndims,
               const int *dimids,
               int compress ) {
  int rc;
  int varid = -1;
  rc = nc_def_var(ncid, name, nc_type, ndims, dimids, &varid);
  ncerror(rc);
#ifdef USE_NETCDF4
  if (compress) {
    rc = nc_def_var_deflate(ncid, varid, 0, 1, compress);
    ncerror(rc);
  }
#endif
  if (units) {
    rc = nc_put_att_text(ncid, varid, "units", strlen(units), units);
    ncerror(rc);
  }
  if (long_name) {
    rc = nc_put_att_text(ncid, varid, "long_name", strlen(long_name), long_name);
    ncerror(rc);
  }
  if (fill_value) {
    switch(nc_type) {
    case NC_SHORT:
      rc = nc_put_att_short(ncid, varid, "missing_value", NC_SHORT, 1, (short *)fill_value);
      ncerror(rc);
      rc=nc_put_att_short(ncid, varid, "_FillValue", NC_SHORT, 1, (short *)fill_value);
      ncerror(rc);
      break;
    case NC_INT:
      rc = nc_put_att_int(ncid, varid, "missing_value", NC_INT, 1, (int *)fill_value);
      ncerror(rc);
      rc=nc_put_att_int(ncid, varid, "_FillValue", NC_INT, 1, (int *)fill_value);
      ncerror(rc);
      break;
    case NC_FLOAT:
      rc = nc_put_att_float(ncid, varid, "missing_value", NC_FLOAT, 1, (float *)fill_value);
      ncerror(rc);
      rc=nc_put_att_float(ncid, varid, "_FillValue", NC_FLOAT, 1, (float *)fill_value);
      ncerror(rc);
      break;
    case NC_DOUBLE:
      rc = nc_put_att_double(ncid, varid, "missing_value", NC_DOUBLE, 1, (double *)fill_value);
      ncerror(rc);
      rc=nc_put_att_double(ncid, varid, "_FillValue", NC_DOUBLE, 1, (double *)fill_value);
      ncerror(rc);
      break;
    default:
      fprintf(stderr, "Error: dont know how to put missing value for variable %s\n",
              name);
      exit(1);
      break;
    }
  }
  return varid;
}
                          

static Cdf *create_cdf(const char *filename,
                       const char *clm_filename,
                       Coord res,
                       const Coord_array *array,
                       int compress)
{
  Cdf *cdf;
  float *lon, *lat;
  int miss = OCEAN_CELL;
  short missdir = NC_FILL_SHORT, zero = 0;
  //int missint = NC_FILL_INT;
  //float missfloat = NC_FILL_FLOAT;
  int i,rc,dim[2];
  String s;
  time_t t;
  int lat_var_id,lon_var_id,lat_dim_id,lon_dim_id;
  //int linelen_dim_id, routes_dim_id; // record dimension is number of routes
  //int basinbounds_dim_id;

  cdf=new(Cdf);
  lon=newvec(float,array->nlon);
  if (lon==NULL) { printallocerr("lon"); return NULL; }
  lat=newvec(float,array->nlat);
  if (lat==NULL) { printallocerr("lat"); return NULL; }
  lon[0] = (float)array->lon_min;
  for (i=1; i<array->nlon; i++)
    lon[i] = lon[i-1]+(float)res.lon;
  lat[0] = (float)array->lat_min;
  for (i=1; i<array->nlat; i++)
    lat[i] = lat[i-1]+(float)res.lat;

  cdf->index=array;
#ifdef USE_NETCDF4
  rc=nc_create(filename,(compress) ? NC_CLOBBER|NC_NETCDF4 : NC_CLOBBER,&cdf->ncid);
#else
  rc=nc_create(filename,NC_CLOBBER,&cdf->ncid);
#endif
  if(rc)
  {
    fprintf(stderr,"ERROR426: Cannot create file '%s': %s.\n",
            filename,nc_strerror(rc));
    free(lon);
    free(lat);
    free(cdf);
    return NULL;
  }
  rc=nc_def_dim(cdf->ncid,LAT_DIM_NAME,array->nlat,&lat_dim_id);
  ncerror(rc);
  rc=nc_def_dim(cdf->ncid,LON_DIM_NAME,array->nlon,&lon_dim_id);
  ncerror(rc);

  snprintf(s,STRING_LEN,"%s %s",myname, clm_filename);
  rc=nc_put_att_text(cdf->ncid,NC_GLOBAL,"source",strlen(s),s);
  ncerror(rc);
  time(&t);
  snprintf(s,STRING_LEN,"Created for user %s on %s at %s",getuser(),gethost(),
         ctime(&t));
  s[strlen(s)-1]='\0';
  rc=nc_put_att_text(cdf->ncid,NC_GLOBAL,"history",strlen(s),s);
  ncerror(rc);

  lon_var_id = create_cdf_var(cdf->ncid, LON_NAME, "degrees_east", "longitude",
                              NC_FLOAT, NULL, 1, &lon_dim_id, 0);
  lat_var_id = create_cdf_var(cdf->ncid, LAT_NAME, "degrees_north", "latitude",
                              NC_FLOAT, NULL, 1, &lat_dim_id, 0);
  rc=nc_put_att_text(cdf->ncid, lon_var_id,"standard_name",strlen("longitude"),"longitude");
  ncerror(rc);
  rc=nc_put_att_text(cdf->ncid, lon_var_id,"axis",strlen("X"),"X");
  ncerror(rc);
  rc=nc_put_att_text(cdf->ncid, lat_var_id,"standard_name",strlen("latitude"),"latitude");
  ncerror(rc);
  rc=nc_put_att_text(cdf->ncid, lat_var_id,"axis",strlen("Y"),"Y");
  ncerror(rc);
  dim[0]=lat_dim_id;
  dim[1]=lon_dim_id;

  cdf->vardest_iid = create_cdf_var(cdf->ncid, "dest_lon", NULL, "longitude of destination cell",
                              NC_INT, &miss, 2, dim, compress);
  cdf->vardest_jid = create_cdf_var(cdf->ncid, "dest_lat", NULL, "latitude of destination cell",
                              NC_INT, &miss, 2, dim, compress);

  cdf->vardrain2ocean_id = create_cdf_var(cdf->ncid, "drain2ocean", NULL, NULL,
                                          NC_SHORT, &zero, 2, dim, compress);
  cdf->vardirxid = create_cdf_var(cdf->ncid, "dirx", NULL, NULL,
                                  NC_SHORT, &missdir, 2, dim, compress);
  cdf->vardiryid = create_cdf_var(cdf->ncid, "diry", NULL, NULL,
                                  NC_SHORT, &missdir, 2, dim, compress);
  rc=nc_put_att_text(cdf->ncid, cdf->vardirxid,"ferret_command", strlen("vector/over/xskip=1/yskip=1 dirx,diry"), "vector/over/xskip=1/yskip=1 dirx,diry");
  ncerror(rc);

  rc=nc_enddef(cdf->ncid);
  ncerror(rc);
  rc=nc_put_var_float(cdf->ncid,lat_var_id,lat);
  ncerror(rc);
  rc=nc_put_var_float(cdf->ncid,lon_var_id,lon);
  ncerror(rc);
  nc_sync(cdf->ncid);
  free(lat);
  free(lon);
  return cdf;
}




static void close_cdf(Cdf *cdf)
{
  nc_close(cdf->ncid);
  free(cdf);
}


int main(int argc, char **argv) {
  FILE *ascfp;
  Coordfile coordfile;
  Coord_array *index;
  Coord *grid,res;
  Cdf *cdf;

  // Metadata from .asc file
  int ncols, nrows;
  float xllcorner, yllcorner, cellsize; // degrees
  int nodata;
  // Data from .asc file
  int *inputmap;
  //short *data;
  float *dest_lon, *dest_lat; // lon/lat of destination cell
  //int *dcount;   // count of cells which drain into this cell
  //int *basin;    // basin number
  //int maxbasin = 0;
  short *drain2ocean; // mask of coast cells which actually darin into the ocean, and thus need to be considered by the FMS-LPJ coupler
  short *dirx, *diry;
  int cell, i,j, ngrid,version,compress,rc;
  //Bool swap, iserror;
  Filename filename;
  size_t routeindex[2] = { 0, 0 };
  //int basinboundscount = 0;

  myname = argv[0];
  if (argc != 4) myusage("wrong number of arguments");

  filename.fmt = CLM;
  filename.name = argv[1];
  coordfile = opencoord(&filename);
  if (NULL == coordfile) {
    fprintf(stderr, "Error opening grid file '%s'.\n", filename.name);
    return EXIT_FAILURE;
  }
  ngrid = numcoord(coordfile);
  res.lon = res.lat = getcellsizecoord(coordfile);
  printf("grid has %d cells, resolution %gx%g\n", ngrid, res.lon, res.lat);
  grid = newvec(Coord,numcoord(coordfile));
  if (NULL == grid) {
    printallocerr("grid");
    return EXIT_FAILURE;
  }
  for (cell = 0; cell < numcoord(coordfile); cell++)
    readcoord(coordfile, grid+cell, res);
  closecoord(coordfile);

  index=createindex(grid,ngrid,res);
  if (NULL == index)
    return EXIT_FAILURE;

  printf("created index from grid file %s: lon_min %g lat_min %g nlon %d nlat %d\n",
         argv[1],index->lon_min, index->lat_min, index->nlon, index->nlat);

  if ((ascfp = fopen(argv[2], "r")) == NULL) {
    perror(argv[2]);
    exit(1);
  } 
  /* read metadata lines from .asc file */
  fscanf(ascfp,"%*s %d",&ncols);
  fscanf(ascfp,"%*s %d",&nrows);
  fscanf(ascfp,"%*s %f",&xllcorner);
  fscanf(ascfp,"%*s %f",&yllcorner);
  fscanf(ascfp,"%*s %f",&cellsize);
  fscanf(ascfp,"%*s %d",&nodata);

  inputmap = newvec(int, ncols*nrows);
  // DDM30.asc is organized N->S , but for NetCDF output we expect S->N order
  // thus read bottom-to-top here.
  for (j = nrows-1 ; j >= 0; j--) {
    for (i = 0; i < ncols; i++) {
      if (1 != fscanf(ascfp, "%d", inputmap+i+j*ncols)) {
        fprintf(stderr, "cannot read code for cell %d,%d from %s\n", i, j, argv[2]);
        perror("fscanf failed");
        exit(1);
      }
    }
  }
  fclose(ascfp);

  if (ncols != (int)(360.0/res.lon)) {
    fprintf(stderr, "Error: number of columns in %s does not match resolution of %s: %d != %d\n",
            argv[1], argv[2], ncols, (int)(360.0/res.lon));
    exit(1);
  }
  if (nrows != (int)(180.0/res.lat)) {
    fprintf(stderr, "Error: number of rows in %s does not match resolution of %s: %d != %d\n",
            argv[1], argv[2], nrows, (int)(180.0/res.lat));
    exit(1);
  }

  /* sigh. here we want a map of the entire world.
   * but createindex() considers only the bounding box of the
   * actually populated land cells.
   * Thus we need to override the results of createindex()
   */
  index->lon_min = -180.0+(res.lon/2.0);
  index->lat_min =  -90.0+(res.lat/2.0);
  index->nlon = ncols; // == (int)(360.0/res.lon);
  index->nlat = nrows; // == (int)(180.0/res.lat);
  printf("patched index lon_min %g lat_min %g nlon %d nlat %d\n",
         index->lon_min, index->lat_min, index->nlon, index->nlat);
  free(index->index); index->index = NULL; // should not be used subsequently. Else let it crash.
  // some lines copied from lpj_trunk/src/netcdf/createindex.c
  index->index=newvec(int,ngrid);
  if (NULL == index->index) { printallocerr("index"); free(index); return EXIT_FAILURE; }
  for(cell=0; cell<ngrid;  cell++) {
    index->index[cell] = (int)((grid[cell].lon-index->lon_min)/res.lon+0.5)+
                      (int)((grid[cell].lat-index->lat_min)/res.lat+0.5)*index->nlon;
#ifndef NDEBUG
    //printf("cell %d lon %g lat %g -> %d,%d = index %d\n",
    //       cell, grid[cell].lon, grid[cell].lat, 
    //       (int)((grid[cell].lon-index->lon_min)/res.lon+0.5),
    //       (int)((grid[cell].lat-index->lat_min)/res.lat+0.5),
    //       index->index[cell]);
    if (index->index[cell]<0 || index->index[cell] >= index->nlon*index->nlat)
      fprintf(stderr,"Invalid index %d.\n", index->index[cell]);
#endif
  }

  //free(grid);

  dest_lon = newvec(float, ncols*nrows);
  if (dest_lon==NULL) { printallocerr("dest_lon"); return EXIT_FAILURE; }
  dest_lat = newvec(float, ncols*nrows);
  if (dest_lat==NULL) { printallocerr("dest_lat"); return EXIT_FAILURE; }
  drain2ocean = newvec(short, ncols*nrows);
  if (drain2ocean==NULL) { printallocerr("drain2ocean"); return EXIT_FAILURE; }
  dirx = newvec(short, ncols*nrows);
  if (dirx==NULL) { printallocerr("dirx"); return EXIT_FAILURE; }
  diry = newvec(short, ncols*nrows);
  if (diry==NULL) { printallocerr("diry"); return EXIT_FAILURE; }
  for (i=0; i<ncols*nrows; i++) {
    dest_lon[i] = dest_lat[i] = OCEAN_CELL; 
    drain2ocean[i] = 0;
    dirx[i] = diry[i] = NC_FILL_SHORT;
  }

  // loop over land cells and re-encode their routing information from DDM30 format into
  // dest_lon/dest_lat destination cell indices.


  // read drainage info for each cell
  // and count number of outflow routes
  for (cell=0; cell<ngrid; cell++) {
    float lon = (float)grid[cell].lon;
    float lat = (float)grid[cell].lat;
    int mapindex = index->index[cell]; // position in output arrays
    if (mapindex < 0 || mapindex >= ncols*nrows) {
      printf("invalid mapindex at %d : %d\n", cell, mapindex);
      exit(EXIT_FAILURE);
    }
    printf("cell %5d %7.2f %6.2f inputmap %d ", cell, lon, lat, inputmap[mapindex]);
    if (-9 == inputmap[mapindex]) {
      // Sigh. DDM30 says this is ocean, but the LPJ grid file says it is land.
      printf("DDM30 ocean cell (-9) inside LPJ land area");
      dest_lon[mapindex] = dest_lat[mapindex] = -333;
      //dirx[mapindex] = diry[mapindex] = NC_FILL_SHORT; // already pre-set above
    } else if (-88 == inputmap[mapindex]) {
      // Same as above, but already detected by a previous run of DDM30 through asc2drainage2asc
      printf("DDM30 ocean cell (-88) inside LPJ land area");
      dest_lon[mapindex] = dest_lat[mapindex] = -333;
      //dirx[mapindex] = diry[mapindex] = NC_FILL_SHORT; // already pre-set above
    } else if (0 == inputmap[mapindex]) {
      printf("DDM30 sink");
      dest_lon[mapindex] = lon; dest_lat[mapindex] = lat; // route to self
      dirx[mapindex] = diry[mapindex] = 0;
    } else if (1 == inputmap[mapindex]) { // East
      dest_lon[mapindex] = (float)(lon + res.lon); if (dest_lon[mapindex] > 180.0) dest_lon[mapindex] -= 360.0;
      dest_lat[mapindex] = lat;
      dirx[mapindex] = 1; diry[mapindex] = 0;
      routeindex[0]++;
    } else if (2 == inputmap[mapindex]) { // SouthEast
      dest_lon[mapindex] = (float)(lon + res.lon); if (dest_lon[mapindex] > 180.0) dest_lon[mapindex] -= 360.0;
      dest_lat[mapindex] = (float)(lat - res.lat);
      dirx[mapindex] = 1; diry[mapindex] = -1;
      routeindex[0]++;
    } else if (4 == inputmap[mapindex]) { // South
      dest_lon[mapindex] = lon;
      dest_lat[mapindex] = (float)(lat - res.lat);
      dirx[mapindex] = 0; diry[mapindex] = -1;
      routeindex[0]++;
    } else if (8 == inputmap[mapindex]) { // SouthWest
      dest_lon[mapindex] = (float)(lon - res.lon); if (dest_lon[mapindex] < -180.0) dest_lon[mapindex] += 360.0;
      dest_lat[mapindex] = (float)(lat - res.lat);
      dirx[mapindex] = -1; diry[mapindex] = -1;
      routeindex[0]++;
    } else if (16 == inputmap[mapindex]) { // West
      dest_lon[mapindex] = (float)(lon - res.lon); if (dest_lon[mapindex] < -180.0) dest_lon[mapindex] += 360.0;
      dest_lat[mapindex] = lat;
      dirx[mapindex] = -1; diry[mapindex] = 0;
      routeindex[0]++;
    } else if (32 == inputmap[mapindex]) { // NorthWest
      dest_lon[mapindex] = (float)(lon - res.lon); if (dest_lon[mapindex] < -180.0) dest_lon[mapindex] += 360.0;
      dest_lat[mapindex] = (float)(lat + res.lat);
      dirx[mapindex] = -1; diry[mapindex] = 1;
      routeindex[0]++;
    } else if (64 == inputmap[mapindex]) { // North
      dest_lon[mapindex] = lon;
      dest_lat[mapindex] = (float)(lat + res.lat);
      dirx[mapindex] = 0; diry[mapindex] = 1;
      routeindex[0]++;
    } else if (128 == inputmap[mapindex]) { // NorthEast
      dest_lon[mapindex] = (float)(lon + res.lon); if (dest_lon[mapindex] > 180.0) dest_lon[mapindex] -= 360.0;
      dest_lat[mapindex] = (float)(lat + res.lat);
      dirx[mapindex] = 1; diry[mapindex] = 1;
      routeindex[0]++;
    } else printf("ERROR invalid inputmap %d!\n", inputmap[mapindex]);
    printf(" -> %7.2f %6.2f\n", dest_lon[mapindex], dest_lat[mapindex]);
  } // loop over cells
  printf("found %d outflow routes\n", (int)(routeindex[0]));


  // now we have all the information to create NetCDF file
  cdf = create_cdf(argv[3],argv[2],res,index, 0);
  if (NULL == cdf) return EXIT_FAILURE;

  // Now all cells outside the gird-file land mask should be marked
  // as ocean cells -999, and land cells with either -333 or -180...180 in dest_lon.
  // We can now search for those coast cells which actually drain into the ocean.
  for (cell=0; cell<ngrid; cell++) {
    int mapindex = index->index[cell]; // position in output arrays
    i = (int)((grid[cell].lon-index->lon_min)/res.lon+0.5);
    j = (int)((grid[cell].lat-index->lat_min)/res.lat+0.5);
    if (mapindex < 0 || mapindex >= ncols*nrows) {
      printf("invalid mapindex at %d : %d\n", cell, mapindex);
      exit(EXIT_FAILURE);
    }
    assert(i+j*ncols == mapindex);
    // only look at cells with actual outflow
    if (inputmap[mapindex] > 0) {
      int dest_i = i+dirx[mapindex];
      int dest_j = j+diry[mapindex];
      int destindex;
      if (dest_i >= ncols) dest_i -= ncols;
      if (dest_i < 0) dest_i += ncols;
      destindex = dest_i+dest_j*ncols;
      if (destindex < 0 || destindex >= ncols*nrows) {
        printf("invalid destindex at %d : %d\n", cell, destindex);
        exit(EXIT_FAILURE);
      }
      if (OCEAN_CELL == dest_lon[destindex]) {
        drain2ocean[mapindex] = 1;
      }
    }
  }


  rc = nc_put_var_float(cdf->ncid, cdf->vardest_iid, dest_lon);
  ncerror(rc);
  rc = nc_put_var_float(cdf->ncid, cdf->vardest_jid, dest_lat);
  ncerror(rc);
  rc = nc_put_var_short(cdf->ncid, cdf->vardrain2ocean_id, drain2ocean);
  ncerror(rc);
  rc = nc_put_var_short(cdf->ncid, cdf->vardirxid, dirx);
  ncerror(rc);
  rc = nc_put_var_short(cdf->ncid, cdf->vardiryid, diry);
  ncerror(rc);
#if 0
  rc = nc_put_var_int(cdf->ncid, cdf->dcount_var_id, dcount);
  ncerror(rc);
  rc = nc_put_var_int(cdf->ncid, cdf->basin_var_id, basin);
  ncerror(rc);
#endif

  close_cdf(cdf);
  // the main purpose of these free() invocations is to catch out-of-bound errors in the preceding code
  free(grid);
  free(dest_lon);
  free(dest_lat);
  free(drain2ocean);
  free(dirx);
  free(diry);
  free(index->index);
  free(index);
  return EXIT_SUCCESS;
} /* of 'main' */