#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <float.h>
#include "bds.h"
#include "gdstool.h"
#include "gds.h"
#include "gribw.h"
#include "angles.h"

/*
 * ggrib -- Shrinks GRIB-files to selected geographic subdomain
 *
 *  Compile: make -f ggrib.make
 *  Usage: type ggrib without arguments
 *  Requires: libgribw
 *
 *  Oyvind.Breivik@dnmi.no
 *
 */

/* Constants */
#define VERSION "ggrib v2.6, 2001-11-05, Oyvind.Breivik@dnmi.no"
#define USAGE "Usage: %s in.grb out.grb westlon southlat eastlon northlat\n"
#define MDEG 1000.0
#define N 100000
#define TOL 0.001

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

   unsigned char *pds, *gds, *bms, *bds;
   float *arr, *arr2; /* tmp work arrays */

   int nrec=0;     /* records */
   int nxny;       /* grid dimensions */
   int arr_size=N; /* default arr dimension */

   long int len_grib, pos=0;

   int nx, ny, nx2, ny2;                 /* grid dimensions, old and new */
   int i, j, k, l, m, n, i1, j1, i2, j2, ii, jj, mm, nn, scan; /* indices */
   int idir, jdir, iadj, cyclic, global; /* logical vars */

   float lon1, lon2, lat1, lat2;
   float dlon, dlat;                               /* from orig GDS */
   float firstlon2, lastlon2, firstlat2, lastlat2; /* for new GDS */
   float wlon, elon, slat, nlat;                   /* orig grid boundaries */
   float wlon2, elon2, slat2, nlat2;               /* new grid boundaries */

   unsigned char mode; /* GDS octet 17 */

   FILE *input, *output;

   /* Help line */
   if (argc == 1) {
      fprintf(stderr, "%s\n", VERSION);
      fprintf(stderr,"Shrinks lon/lat GRIB records to the domain "
      "(westlon southlat eastlon northlat)\n");
      fprintf(stderr, USAGE, argv[0]);
      fprintf(stderr,"  Bounding coords can be left out by typing a dot (.), ");
      fprintf(stderr,"for example:\n\n    ggrib in.grb out.grb -25 . . 78\n\n");
      fprintf(stderr,"  will produce a grid from 25 E and southern rim of orig "
      "grid to eastern rim of\n  orig grid and 78 N.\n");
      exit(8);
   }

   /* Open files */
   if (argc < 7) {
      fprintf(stderr, USAGE, argv[0]);
      exit(8);
   }
   if ((input = fopen(argv[1],"rb")) == NULL) {
      fprintf(stderr,"Could not open file: %s\n", argv[1]); 
      exit(7);                                              
   }                                                         
   if ((output = fopen(argv[2],"wb")) == NULL) {             
      fprintf(stderr,"Could not open file: %s\n", argv[2]); 
      exit(7);                                              
   }                                                         

   /* Preallocate BIG arrays for speed */
   arr = (float *) malloc(N * sizeof(float));
   arr2 = (float *) malloc(N * sizeof(float));
   
   /* Loop over input records */
   while ((len_grib = rd_grib_msg(input, pos, &pds, &gds, &bms, &bds)) > 0) {

      nrec++; /* count records */

      /* Reallocate arrays if nxny > N */
      nxny = get_nxny(pds, gds, bms, bds);
      if (nxny > arr_size) {
         free(arr);
         free(arr2);
         arr_size = nxny;
         arr =  (float *) malloc(arr_size * sizeof(float));
         arr2 = (float *) malloc(arr_size * sizeof(float));
         if (arr == NULL || arr2 == NULL) {
            fprintf(stderr,"Memory allocation failed\n");
         }
      } /* end if nxny */

      /* Unpack data, get_unpk_bds doesn't seem to work properly */
      unpk_bds(arr, pds, gds, bms, bds, nxny);

      /* Is grid in geographic coordinates? */
      if (!GDS_LatLon(gds)) {
         fprintf(stderr,\
         "Grid #%4d not in geographic coordinates (lon,lat)\n",nrec);
         exit(9);
      } 

      /* Original grid parameters */
      scan = GDS_LatLon_scan(gds);
      nx = GDS_LatLon_nx(gds);
      ny = GDS_LatLon_ny(gds);
      dlon = GDS_LatLon_dx(gds)/MDEG;
      dlat = GDS_LatLon_dy(gds)/MDEG;

      /* Grid orientation */
      idir = ((scan&128)==128);
      idir = 1-idir;          /* West to east is "true" */
      jdir = ((scan&64)==64); /* South to north is "true" */
      iadj = ((scan&32)==32);
      iadj = 1-iadj;          /* Zonal points adjacent is "true" */

      if (idir) {
         wlon = GDS_LatLon_Lo1(gds)/MDEG; /* grid running W to E (default)? */
         elon = GDS_LatLon_Lo2(gds)/MDEG;
      }
      else {
         elon = GDS_LatLon_Lo1(gds)/MDEG; /* or E to W? */
         wlon = GDS_LatLon_Lo2(gds)/MDEG;
      }
         
      if (jdir) {
         slat = GDS_LatLon_La1(gds)/MDEG; /* grid running S to N? */
         nlat = GDS_LatLon_La2(gds)/MDEG;
      }
      else {
         nlat = GDS_LatLon_La1(gds)/MDEG; /* or N to S (default)? */
         slat = GDS_LatLon_La2(gds)/MDEG;
      }

      if (nlat < slat) {
         fprintf(stderr,"Error: scan direction mismatches with south and "
         "north grid boundaries, scan=%d\n",scan);
         exit(5);
      }

      /* User input */
      if (argv[3][0]=='.') {
         lon1=wlon;
      }
      else {
         lon1 = (float) atof(argv[3]); /* western boundary */
      }
      if (argv[4][0]=='.') {
         lat1=slat;
      }
      else {
         lat1 = (float) atof(argv[4]); /* southern boundary */
      }
      if (argv[5][0]=='.') {
         lon2=elon;
      }
      else {
         lon2 = (float) atof(argv[5]); /* eastern boundary */
      }
      if (argv[6][0]=='.') {
         lat2=nlat;
      }
      else {
         lat2 = (float) atof(argv[6]); /* northern boundary */
      }

      /* Grid around the world in longitude? */
      cyclic = (abs(ang180(elon+dlon-wlon)) < TOL);

      /* Confine boundaries to grid domain */
      if (!cyclic) {
         /* Western lon outside grid? */
         if (ang360(lon1-wlon) > ang360(elon-wlon)) {
            lon1=wlon;
         }
         /* Eastern lon outside grid? */
         if (ang360(lon2-wlon) > ang360(elon-wlon)) {
            lon2=elon;
         }
         /* Eastern and western lons swapped round? */
         if (ang360(lon1-wlon) > ang360(lon2-wlon)) {
            lon2 = elon; /* Retain eastern part of grid and bitch */
            fprintf(stderr,"Error, boundary crossed in non-periodic grid "
            "record no %d ... retaining eastern part of grid\n",nrec);
         }
      } /* end if cyclic */
      wlon2 = lon1;
      elon2 = lon2;

      slat2 = lat1 < slat ? slat : lat1;
      nlat2 = lat2 > nlat ? nlat : lat2;

      /* Find bounding grid indices */
      i1 = (int) (ang360(wlon2-wlon)/dlon + 1.5); /* round index to nearest */
      i2 = (int) (ang360(elon2-wlon)/dlon + 1.5);
      j1 = (int) ((slat2-slat)/dlat + 1.5);
      j2 = (int) ((nlat2-slat)/dlat + 1.5);

      /* Find new grid dimensions */
      nx2 = (int) (ang360(elon2-wlon2)/dlon + 1.5);
      ny2 = j2-j1+1;

      /* Make new grid global in lon if lon1==lon2 */
      global = ((nx2==1) && cyclic);
      if (global) {
         nx2 = nx;
         i2 = i1>1 ? i1-1 : nx2; /* i2 = i1-1 with wrap-around */
      } 
      else {
         nx2 = i2-i1+1;
         if (nx2 < 1) {
            nx2 = nx+nx2; /* cyclic grids are tricky bastards */
         }
      }

      /* Adjust boundaries to nearest exact grid points */
      wlon2 = ang360(wlon+((float)(i1-1))*dlon); /* find exact grid longitude */
      elon2 = ang360(wlon+((float)(i2-1))*dlon); 
      slat2 = slat+((float)(j1-1))*dlat;
      nlat2 = slat+((float)(j2-1))*dlat;

      /* Shrink grid */
      for (m=1; m<=nx2; m++) {    /* lons */
         for (n=1; n<=ny2; n++) { /* lats */
            i = (i1+m-2)%nx + 1;  /* wrap around if grid is periodic in lon*/
            j = j1+n-1;           /* non-periodic in lat */

            /* Positive direction index i? */
            ii = i;
            mm = m;
            if (!idir) {
               ii = nx-i+1;
               mm = nx2-m+1;
            }

            /* Positive direction index j? */
            jj = j;
            nn = n;
            if (!jdir) {
               jj = ny-j+1;
               nn = ny2-n+1;
            }

            /* Index i adjacent (grid running east-west)? */
            if (iadj) {
               jj = (jj-1)*nx;
               nn = (nn-1)*nx2;
            }
            else {
               ii = (ii-1)*ny;  /* index j adjacent, running north-south */
               mm = (mm-1)*ny2;
            }

            k = ii+jj;
            l = mm+nn;
            arr2[l-1] = arr[k-1];
         } /* end for n */
      } /* end for m */

      /* Make sure longitudes lie in [-180,180) */
      wlon2 = ang180(wlon2);
      elon2 = ang180(elon2);

      /* Cosmetics: change e.g. range 150:-150 deg to read 150:210 deg */
      if ( (ang360(wlon2)<180.0) && (ang360(elon2)>=180.0) ) {
        elon2 = ang360(elon2);
      }

      /* Cosmetics: last longitude should not be -180 */
      if (elon2 == -180.0) {
         elon2 = 180.0;
      }

      /* Compute new grid boundaries in lon/lat */
      if (idir) {
         firstlon2 = wlon2;  /* grid starts on western rim (default) */
         lastlon2 = elon2;
      }
      else {
         firstlon2 = elon2;
         lastlon2 = wlon2;
      }
         
      if (jdir) {
         firstlat2 = slat2;  /* grid starts on southern rim */
         lastlat2 = nlat2;
      }
      else {
         firstlat2 = nlat2;  /* default */
         lastlat2 = slat2;
      }

      /* Remake grid definition section, GDS */
      mode = gds[16];
      gds = new_LatLon_GDS(pds,nx2,ny2,firstlon2,firstlat2,lastlon2,
                           lastlat2,dlon,dlat);
      gds[16] = mode; /* Restore the mode flags */

      /* Remake bitmap section, BMS */
      nxny = nx2*ny2;
      bms = mk_BMS(pds, arr2, &nxny, UNDEFINED_LOW, UNDEFINED_HIGH);

      /* Remake binary data section, BDS */
      if (BDS_BinScale(bds)) { /* ECMWF style? */
         set_BDSMinBits(BDS_NumBits(bds));
         set_BDSMaxBits(BDS_NumBits(bds));
      }
      bds = mk_BDS(pds, arr2, nxny);

      /* Write new GRIB record */
      wrt_grib_msg(output, pds, gds, bms, bds);

      /* Clean up */
      if (gds != NULL) {
        free(gds);
        gds = NULL;
      }
      if (bms != NULL) {
        free(bms);
        bms = NULL;
      }
      if (bds != NULL) {
        free(bds);
        bds = NULL;
      }

      pos += len_grib; /* file position */
   } /* end for nrec */

   /* Close files */
   fclose(input);
   fclose(output);

   /* Clean up */
   if (arr != NULL) {
     free(arr);
     arr = NULL;
   }
   if (arr2 != NULL) {
     free(arr2);
     arr2 = NULL;
   }

   /* Output */
   if (nrec > 1) {
      printf("Wrote %d records\n", nrec);
   }
   else if (nrec == 1) {
      printf("Wrote one record\n");
   }
   else {
      fprintf(stderr,"ERROR: no records written\n");
   } /* end if */

   return 0;
} /* end main */