/*
 * printrivermap.c
 *
 * PIK/S. Petri/May 2023
 * $Id$
 *
 * Print a land_lad/soil/river.F90 river_destination_field
 *
 * The file river_destination_field contains two times 96x80
 * values (2*96*80=2*7680).
 * (802 = 2+2*400 lines, 15364 = 3+1+7680*2 words)
 *
 * --------
 96,  80,  -1.876
   (20i4)
  85  85  85  85  85  85  85  85  85  85  84  84  84  84  84  84  84  84  84  84
  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84  84
 * --------
 *
 * Apparently, there is one pair of values for each surface cell.
 * The first value denotes column, the second value denotes row index of the
 * destination cell.
 * Cells which point to themselfes are ocean cells.
 * The input data is regridded dynamically to the actually used model
 * grid resolution.
 * After regridding, for each land cell it is made sure that its runoff
 * goes into an ocean cell.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

static char *myname = NULL;

double eps = 1e-4;

static void
usage(const char *myname, const char *msg) {
  if (msg) fprintf(stderr, "%s\n", msg);
  fprintf(stderr, "usage: %s input_file\n", myname);
  exit(1);
}

// A macro to convert from 2-D array indices in Fortran convention to
// 1-D C-array index
#define A(_i, _j) ((_i-1)+nlons*(_j-1))


int main(int argc, char **argv) {

  int i, j, I;
  size_t nlons = -1, nlats = -1;
  double wb_degrees, dlon, dlat;
  FILE *infile;
  myname = argv[0];
  int *dest_i, *dest_j; // indices
  int *dcount, maxdcount, *dpoint, *basin, nbasins;
  double *lons, *lats; // coordinates
  char input_format[80];

  if (argc != 2) usage(myname, "need 1 parameter");

  infile = fopen(argv[1], "r");
  if (NULL == infile) { perror(argv[1]); exit(1); }
  fscanf(infile, "%ld,%ld,%lg\n", &nlons, &nlats, &wb_degrees);
  printf("nlons %lu nlats %lu western boundary at %lg\n", nlons, nlats, wb_degrees);
  fscanf(infile, "%s\n", input_format);
  // test array access macro
  for (j = 1; j <= nlats; j++) {
    for (i = 1; i <= nlons; i++) {
      I = A(i,j);
      assert(I>=0);
      assert(I<nlons*nlats);
    }
  }
  

  dest_i = (int *)malloc(nlons*nlats*sizeof(int));
  dest_j = (int *)malloc(nlons*nlats*sizeof(int));
  dcount = (int *)malloc(nlons*nlats*sizeof(int));
  dpoint = (int *)malloc(nlons*nlats*sizeof(int));
  basin = (int *)malloc(nlons*nlats*sizeof(int));
  lons = (double *)malloc(sizeof(double)*nlons);
  lats = (double *)malloc(sizeof(double)*nlats);
  dlon = 360.0/nlons;
  dlat = 180.0/nlats;
  for (i = 0; i < nlons; i++) {
    lons[i] = wb_degrees + i*dlon;
  }
  for (j = 0; j < nlats; j++) {
    //j=0 -> North Pole j=nlats-1 -> South Pole
    //lats[j] = 90-dlat/2.0-dlat*j;
    //j=0 -> South Pole j=nlats-1 -> North Pole
    lats[j] = dlat*j-90.0+(dlat/2.0);
  }
  // use loop indices in Fortran order, 1..n
  // subtract 1 for actual access to C arrays
  for (j = 1; j <= nlats; j++) {
    for (i = 1; i <= nlons; i++) {
      int n = fscanf(infile, "%4d", &dest_i[A(i,j)]);
      if (n < 1) { printf("Error: could not read dest_i(%d, %d)\n", i, j); }
    }
  }
  for (j = 1; j <= nlats; j++) {
    for (i = 1; i <= nlons; i++) {
      int n = fscanf(infile, "%4d", &dest_j[A(i,j)]);
      if (n < 1) { printf("Error: could not read dest_j(%d, %d)\n", i, j); }
    }
  }

  //for (j = 1; j <= nlats; j++) {
  //  for (i = 1; i <= nlons; i++) {
  //    printf("%4d/%4d, ", dest_i[A(i,j)],dest_j[A(i,j)]);
  //  }
  //  printf("\n");
  //}

  maxdcount = 0;
  nbasins = 0;

  for (j = nlats; j >= 1; j--) {
    printf("%6.2f ",lats[j-1]);
    for (i = 1; i <= nlons; i++) {
      if (i == dest_i[A(i,j)] && j == dest_j[A(i,j)])
        printf("."); // ocean cell routes into itself
      else
        printf("@"); // a land cell?
      dcount[A(i,j)]  = 0;
      dpoint[A(i,j)]  = 0;
      basin[A(i,j)]   = 0;
    }
    printf("\n");
  }
  printf("\n");

  for (j = nlats; j >= 1; j--) {
    for (i = 1; i <= nlons; i++) {
      I = A(i,j);
      assert(I >= 0);
      assert(I <= nlons*nlats);
      int id = dest_i[A(i,j)];
      int jd = dest_j[A(i,j)];
      if (id != i || jd !=j) { // point does not discharge to itself => it is not ocean
        I = A(id,jd);
        assert(I >= 0);
        assert(I <= nlons*nlats);
        if (dpoint[A(id,jd)] == 0) {
          // we have not encountered this destination point before,
          // add a basin
          nbasins       = nbasins+1;
          dpoint[A(id,jd)] = nbasins;
        }
        //printf("from basin %d to %d\n", basin[A(i,j)], dpoint[A(id,jd)]);
        basin[A(i,j)] = dpoint[A(id,jd)];
        dcount[A(id,jd)]++;
        if (maxdcount < dcount[A(id,jd)]) maxdcount = dcount[A(id,jd)];
      }
    }
  }
  printf("\nfound %d basins, max dcount %d\n", nbasins, maxdcount);
  if (0 == nbasins) {
    printf("No Basins???\n");
    exit(0);
  }

  char charmap[] = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";
  int nchars = sizeof(charmap)-3;

  printf("\nBasins:\n");
  for (j = nlats; j >= 1; j--) {
    printf("%6.2f ",lats[j-1]);
    for (i = 1; i <= nlons; i++) {
      //if (0 == basin[A(i,j)]) printf("."); else printf("B");
      int mapindex = (int)(nchars * (double)(basin[A(i,j)]) / (double)nbasins);
      if (basin[A(i,j)] > 0) mapindex++; // make 0 uniq
      //printf("basin %d mapindex %d char %c\n", basin[A(i,j)], mapindex, charmap[mapindex]);
      printf("%c", charmap[mapindex]);
    }
    printf("\n");
  }

  printf("\nDestination count:\n");
  for (j = nlats; j >= 1; j--) {
    printf("%6.2f ",lats[j-1]);
    for (i = 1; i <= nlons; i++) {
      //if (0 == dcount[A(i,j)]) printf("."); else printf("D");
      int mapindex = (int)(nchars * (double)(dcount[A(i,j)]) / (double)maxdcount);
      if (dcount[A(i,j)] > 0) mapindex++; // make 0 uniq
      //printf("dcount %d mapindex %d char %c\n", dcount[A(i,j)], mapindex, charmap[mapindex]);
      printf("%c", charmap[mapindex]);
    }
    printf("\n");
  }

  return 0;
}