/*
 * createriverrouting.c
 *
 * PIK/S. Petri/May 2023
 * $Id$
 *
 * Create a river destination field in a format that can be used by
 * land_lad/soil/rivers.F90 for the idealised planet with one or two
 * circumpolar stripes of land.
 *
 * 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>

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 number-of-stripes equatorward-lat stripewidth output-file\n",
          myname);
  fprintf(stderr, "  number-of-stripes=1 : only in northern hemisphere");
  fprintf(stderr, "  number-of-stripes=2 : symmetric in both hemispheres");
  fprintf(stderr, "example: %s 2 10 20 rivers-10-20\nproduces a map for two stripes from 10N to 30N and 10S to 30S\n", myname);
  exit(1);
}



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

  int i, j;
  size_t nlons = 144, nlats = 90; /* hardcoded resolution here */
  double dlat=180.0/nlats;
  int nstripes = -1;
  double equatorlat, stripewidth, polelat, midlat;
  FILE *outfile;
  myname = argv[0];
  // In the StripedPlanet each cell drains into a cell of same longitude.
  // Thus we dont need longitude coordinates nor destination columns.
  int *destrow;
  int d1 = -1, d2 = -1, d3 = -1, d4 = -1; // latitude indices of the 4 coastlines, S to N
  double *lats; // latitude coordinates

  if (argc != 5) usage(myname, "need 4 parameters");

  nstripes = atoi(argv[1]);
  equatorlat = atof(argv[2]);
  stripewidth = atof(argv[3]);

  if (nstripes < 1 || nstripes > 2) usage(myname, "number-of-stripes must be 1 or 2");
  if (equatorlat < 0 || equatorlat > 90) usage(myname, "equatorlat must be between 0 and 90");
  if (stripewidth < 0) usage(myname, "stripewidth must be positive");
  polelat = equatorlat+stripewidth;
  if (polelat < 0 || polelat > 90)
    usage(myname, "equatorlat+stripewidth must be between 0 and 90");
  midlat = equatorlat + 0.5*stripewidth; // mid of land stripe

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

  lats = (double *)malloc(sizeof(double)*nlats);
  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);
    //printf("%3d %6.2f\n", j, lats[j]);
  }

  printf("equatorlat %f midlat %f polelat %f\n",
         equatorlat, midlat, polelat);
  for (j = 0; j < nlats-1; j++) {
    if (2 == nstripes) { // southern stripe
      if (lats[j] <= -polelat    && -polelat    <= lats[j+1]) { d1 = j; /*printf("d1 ");*/ }
      if (lats[j] <= -equatorlat && -equatorlat <= lats[j+1]) { d2 = j; /*printf("d2 ");*/ }
    }
    if (lats[j] <= equatorlat && equatorlat <= lats[j+1]) { d3 = j; /*printf("d3 ");*/ }
    if (lats[j] <= polelat    && polelat    <= lats[j+1]) { d4 = j; /*printf("d4 ");*/ }
  }
  if (2 == nstripes)
    printf("d1 %d %f d2 %d %f", d1, lats[d1], d2, lats[d2]);
  else
    printf("d1 -1 d2 -1");
  printf(" d3 %d %f d4 %d %f", d3, lats[d3], d4, lats[d4]);
  printf("\n");

  destrow = (int *)malloc(sizeof(int)*nlons*nlats);
  for (j = 0; j < nlats; j++) {
    int d = j; // ocean cell routes to itself
    if (-polelat <= lats[j] && lats[j] <= -midlat     && d1 > 0) d = d1;
    if (-midlat  <  lats[j] && lats[j] <= -equatorlat && d2 > 0) d = d2;
    if (equatorlat <= lats[j] && lats[j] <= midlat    && d3 > 0) d = d3;
    if (midlat     <  lats[j] && lats[j] <= polelat   && d4 > 0) d = d4;
    printf("j %3d lats[j] %6.2f d %3d\n", j, lats[j], d);
    for (i = 0; i < nlons; i++) {
      destrow[j*nlons+i] = d;
    }
  }

  fprintf(outfile, "%3d, %3d, %6.3f\n", (int)nlons, (int)nlats, -360.0/nlons/2.0);
  fprintf(outfile, "(20i4)\n");
  for (j = 0; j <nlats; j++) { // write out in N to S order
    int outcol = 0;
    for (i = 0; i < nlons; i++) {
      fprintf(outfile, "%4d", i+1); // destination columns run 1..nlons!
      if (20 == ++outcol) { fprintf(outfile, "\n"); outcol = 0; }
    }
    fprintf(outfile, "\n");
  }
  for (j = 0; j < nlats; j++) {
    int outcol = 0;
    for (i = 0; i < nlons; i++) {
      fprintf(outfile, "%4d", destrow[j*nlons+i]+1); // destination rows run 1..nlats
      if (20 == ++outcol) { fprintf(outfile, "\n"); outcol = 0; }
    }
    fprintf(outfile, "\n");
  }
  return 0;
}