/*
 * suppressriverrouting.c
 *
 * PIK/S. Petri/Jul 2023
 * $Id$
 *
 * Modify a river destination field so that runoff from Antarctica is suppressed.
 *
 * 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 output-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, *outfile;
  myname = argv[0];
  int *dest_i, *dest_j; // indices
  double *lons, *lats; // coordinates
  char input_format[80];

  if (argc != 3) usage(myname, "need 2 parameters");

  outfile = fopen(argv[2], "w");
  if (NULL == outfile) { perror(argv[2]); exit(1); }

  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));
  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");
  //}

  // South to North
  for (j = 1; j <= nlats; j++) {
    if (lats[j-1] <= -60) {  // south of Africa?
      for (i = 1; i <= nlons; i++) {
        if (i != dest_i[A(i,j)] || j != dest_j[A(i,j)]) { // not an ocean cell?
          //dest_i[A(i,j)] = i; dest_j[A(i,j)] = j; // suppress!
          dest_i[A(i,j)] = nlons/2; dest_j[A(i,j)] = 1; // all into south pole
        }
      }
    }
  }

  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?
    }
    printf("\n");
  }
  printf("\n");


  fprintf(outfile, "%3d, %3d, %6.3f\n", (int)nlons, (int)nlats, wb_degrees);
  fprintf(outfile, "(20i4)\n");
  for (j = 1; j <= nlats; j++) {
    int outcol = 0;
    for (i = 1; i <= nlons; i++) {
      fprintf(outfile, "%4d", dest_i[A(i,j)]);
      if (20 == ++outcol) { fprintf(outfile, "\n"); outcol = 0; }
    }
    fprintf(outfile, "\n");
  }
  for (j = 1; j <= nlats; j++) {
    int outcol = 0;
    for (i = 1; i <= nlons; i++) {
      fprintf(outfile, "%4d", dest_j[A(i,j)]); // destination rows run 1..nlats
      if (20 == ++outcol) { fprintf(outfile, "\n"); outcol = 0; }
    }
    fprintf(outfile, "\n");
  }

  return 0;
}