# include "/sw/include/gmt.h"

/* Fortran callable routine to read a grd file in pixel registration */
/* June 23, 1995 - David Sandwell */
/* Revised for GMT3.4 December 28, 2002 - David Sandwell */

int readgrd_(rdat,nx,ny,rlt0,rln0,dlt,dln,rdum,title,filein)

  float *rdat;            /* real array for output */
  int *nx;                /* number of x points */
  int *ny;                /* number of y points */
  double *rlt0;            /* starting latitude */
  double *rln0;            /* starting longitude */
  double *dlt;             /* latitude spacing */
  double *dln;             /* longitude spacing */
  double *rdum;            /* dummy value */
  char  *title;           /* title */
  char  *filein;          /* filename of input file */
  
  {
   int i, j;
   int argc2 = 2;
   char *argv2[2];
   struct GRD_HEADER grd;

/* put some dummy data in argv2  to fool GMT_begin*/

   argv2[0]="string1";
   argv2[1]="string2";

/* execute GMT_begin */

   argc2 = GMT_begin (argc2, argv2);

/*  read the header */

   if (GMT_read_grd_info (filein, &grd)) {
      fprintf (stderr, "Error opening file %s\n", filein);
      exit (EXIT_FAILURE);
   }

/*  read the grid */

   if(GMT_read_grd(filein, &grd, rdat, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE )) {
     fprintf (stderr, "Error reading file %s\n", filein);
     exit (EXIT_FAILURE);
   }

   *nx = grd.nx;
   *ny = grd.ny;
   *dlt  = -grd.y_inc;
   *dln  =  grd.x_inc;
   if(grd.node_offset == 1) {
     *rlt0 = grd.y_max;
     *rln0 = grd.x_min; 
   }
   else {
     *rlt0 = grd.y_max + 0.5*grd.y_inc;
     *rln0 = grd.x_min - 0.5*grd.x_inc; 
   }
   *rdum = floor(grd.z_max + 1.);
   strncpy(title,grd.title,80); 

/*  Calculate rdum and reset the dummy values. */

   for (i = 0; i < *nx * *ny; i++) {
       if( isnan (rdat[i]) ) rdat[i] = *rdum; 
   }
  }