#include <stdio.h>
#include "hrdi.h"
#include "netcdf.h"
#include "timesgcm.h"

extern char *malloc();
extern int cdfid;
extern int locint();
extern int htint();

/*
 * global number of recs in data file, number of altitudes, and
 * altitude profile
 */
int nrec,nalt;
float *alt;
int interpslt = INTERPSLT; /* if = 1, interp to slt rather than slon */

char *cmalloc(n,s)
int n; char *s;
{
	char *ptr;
	if (n <= 0) {
	  (void)fprintf(stderr,"cmalloc: bad n=%d\n",n);
	  (void)fprintf(stderr,"%s\n",s);
	  exit(1);
	}
	if ((ptr = malloc(n)) == NULL) {
	  (void)fprintf(stderr,"cmalloc: malloc failed for n=%d\n",n);
	  (void)fprintf(stderr,"%s\n",s);
	  exit(1);
	}
	return(ptr);
}

/*
 * Return longitude corresponding to given ut and local time slt:
 */
float getsltlon(slt,ut)
float slt,ut;
{
	float xlon = (slt-ut)*15.;
	if (xlon >= 180.) xlon -= 360.;
	if (xlon < -180.) xlon += 360.;
	return(xlon);
}

void defdat(itime,rlat,rlon,rsza,rslt,un0,vn0,dat)
int itime;
float rlat,rlon,rsza,rslt,*un0,*vn0;
Dat_struct *dat;
{
	int i;
	dat->ut = (float)itime / 3.6e6;
	dat->glat = rlat;
	dat->glon = rlon;
	if (dat->glon >= 180.) dat->glon = dat->glon - 360.;
	dat->sza = rsza;
	dat->slt = rslt;
	for (i=0; i < nalt; i++) {
	  dat->un[i] = *(un0+i);
	  dat->vn[i] = *(vn0+i);
	}
}

/*
 * Interpolate gcm to uars lat,lon, ht: 
 * Note pnt contains un, vn, and ht in that order at global gcm grid.
 * zpcol[KMX] pressure column field interpolated to slat,lon
 * htzpcol[KMX] = heights at pressure col, interpolated to slat,slon
 * gcm->un, gcm->vn = final gcm interpolated fields to plot
 */
void defgcm(dat,gcm)
Dat_struct *dat;
Gcm_struct *gcm;
{
	float zpcol[KMX], htzpcol[KMX];
/*
 * Save gcm pressure column heights, interpolated to slat,slon, in htzpcol:
 */
	locint(dat->glat,dat->glon,dat->slt,&pnt[2][0][0][0],tgcmhdr->ut,
	  interpslt,IMX,KMX,JMX,gcmlat,gcmlon,htzpcol);
/*
 * Interpolate un to slat,slon (output of locint in zpcol), then to height 
 * scale (output of htint returned in gcm->un,vn, the final plot fields.
 */
	locint(dat->glat,dat->glon,dat->slt,pnt,tgcmhdr->ut,interpslt,
	  IMX,KMX,JMX,gcmlat,gcmlon,zpcol);
	htint(zpcol,htzpcol,0,alt,nalt,KMX,gcm->un,SPVAL);
/*
 * Same for vn:
 */
	locint(dat->glat,dat->glon,dat->slt,&pnt[1][0][0][0],tgcmhdr->ut,
	  interpslt,IMX,KMX,JMX,gcmlat,gcmlon,zpcol);
	htint(zpcol,htzpcol,0,alt,nalt,KMX,gcm->vn,SPVAL);
}

void printdat(ir,dat)
int ir;
Dat_struct *dat;
{
	int i;
	printf("\nrec=%d ut=%6.2f lat,lon=%6.2f %7.2f slt=%5.2f\n",
  	ir,dat->ut,dat->glat,dat->glon,dat->slt);
	printf("un=");
	for (i=0; i < nalt; i++) printf("%8.1f",*(dat->un+i));
	printf("\n");
	printf("vn=");
	for (i=0; i < nalt; i++) printf("%8.1f",*(dat->vn+i));
	printf("\n");
}

FILE *getdatfile(flnmdat)
char *flnmdat;
{
	FILE *fpdat;

	printf("Enter hrdi data file name (d = %s): ",DEFDATFLNM);
  	if ((fscanf(stdin,"%s",flnmdat)) == NULL)
	  printf("getdatfile: error reading flnmdat from stdin\n");
	else if (!strcmp("d",flnmdat)) {
	  strcpy(flnmdat,DEFDATFLNM);
	  printf("Using default dat file name %s\n",flnmdat);
	} else  
	  printf("getdatfile: read from stdin flnmdat = %s\n",flnmdat);
  	if ((fpdat = fopen(flnmdat,"r")) == NULL) {
	  printf("getdatfile: error from fopen opening file %s\n",flnmdat);
	  return(NULL);
	} else
	  return(fpdat);
}

int getgcmfile(flnmgcm)
char *flnmgcm;
{
	int icdf;
	printf("Enter gcm (netcdf) file name (d = %s): ",DEFGCMFLNM);
  	if ((fscanf(stdin,"%s",flnmgcm)) == NULL)
	  printf("getdatfile: error reading flnmgcm from stdin\n");
	else if (!strcmp("d",flnmgcm)) {
	  strcpy(flnmgcm,DEFGCMFLNM);
	} else  
	  printf("getgcmfile: read from stdin flnmgcm = %s\n",flnmgcm);
          icdf = ncopen(flnmgcm,NC_NOWRITE);
	  if (icdf < 0) {
	    printf("getgcmfile: error opening file %s\n",flnmgcm);
	    return(-1);
	  } else
	    return(icdf);
}

main()
{
	FILE *fpdat;
	char *flnmdat, *flnmgcm, *ans, *histname;
	int ir,i,itime, norb,nporb, ire,irb, sizef=sizeof(float);
	int iborb=IBORB, ieorb=IEORB, dogcm;
	float rlat,rlon,rsza,rslt, *un, *vn, sltp;
	Dat_struct *dat;
	Gcm_struct *gcm;
/*
 * Get and open data file under file descriptor fpdat:
 */
	flnmdat = cmalloc(STRLEN,"cb main for flnmdat");
	if ((fpdat = getdatfile(flnmdat)) == NULL)
	  exit(1);
	else
	  printf("main: opened file %s\n",flnmdat);
/*
 * Read number of records (nrec), number of altitudes (nalt),
 * and allocate space for winds profiles:
 */
	if ((fscanf(fpdat,"%d%d",&nrec,&nalt)) == NULL) 
	  printf("main: error reading nrec,nalt from file %s\n",flnmdat);
	else
	  printf("nrec=%d nalt=%d\n",nrec,nalt);
	un = (float*)cmalloc(sizef*nalt,"cb main for un");
	vn = (float*)cmalloc(sizef*nalt,"cb main for vn");
/*
 * Ask about doing gcm, and get netcdf file if we are:
 */
	dogcm = 0;
	histname = "                ";
	ans = cmalloc(1,"cb main for ans");
	printf("\nGet timesgcm data? (y/n): ");
  	if ((fscanf(stdin,"%s",ans)) == NULL)
	  printf("main: error reading answer from stdin\n");
	else if (!strcmp("y",ans)) {
	  printf("Will get timesgcm history\n");
	  dogcm = 1;
	  flnmgcm = cmalloc(STRLEN,"cb main for flnmgcm");
  	  if ((cdfid = getgcmfile(flnmgcm)) < 0) {
	    printf("main: error from getgcmfile\n");
	    exit(1);
	  } else {
	    tail_(flnmgcm,histname,strlen(flnmgcm),strlen(histname));
	    for (i=0; i < strlen(histname); i++) 
	      if (*(histname+i) == '.') *(histname+i) = '\0';
	    printf("Reading netCDF file %s (histname=%s)...\n",flnmgcm,
	      histname);
	    rdcdf(cdfid);
	    printf("Completed reading netcdf file\n");
	  }
	  gcm = (Gcm_struct*)cmalloc(sizeof(Gcm_struct)*MXPORB,
		"cb main for gcm");
	  for (i=0; i < MXPORB; i++) {
	    (gcm+i)->un = (float*)cmalloc(sizef*nalt,"cb main for gcm->un");
	    (gcm+i)->vn = (float*)cmalloc(sizef*nalt,"cb main for gcm->vn");
	  }
	}
/*
 * Allocate space for data structures:
 */
	dat = (Dat_struct*)cmalloc(sizeof(Dat_struct)*MXPORB,"cb main for dat");
	for (i=0; i < MXPORB; i++) {
	  (dat+i)->un = (float*)cmalloc(sizef*nalt,"cb main for dat->un");
	  (dat+i)->vn = (float*)cmalloc(sizef*nalt,"cb main for dat->vn");
	}
/* 
 * Read altitude profile:
 */
	alt = (float*)cmalloc(sizef*nalt,"cb main for alt");
	printf("alt = ");
  	for (i=0; i < nalt; i++) { 
	  if ((fscanf(fpdat,"%f",alt+i)) == NULL)
	    printf("main: error reading record %i\n",i);
	  printf("%8.1f",*(alt+i));
	}
	printf("\n");
/*
 * Record loop:
 */
	norb = nporb = 0; sltp = 0.;
	for (ir=0; ir < nrec; ir++) {
/*
 * Read data for this record:
 */
	  if ((fscanf(fpdat,"%d%f%f%f%f",&itime,&rlat,&rlon,&rsza,&rslt)) == NULL)
	    printf("main: error reading record %d\n",ir);
  	  for (i=0; i < nalt; i++) { 
	    if ((fscanf(fpdat,"%f",un+i)) == NULL) 
	      printf("main: error reading un for record %d\n");
	  }
  	  for (i=0; i < nalt; i++) { 
	    if ((fscanf(fpdat,"%f",vn+i)) == NULL) 
	      printf("main: error reading vn for record %d\n");
	  }
/*
 * Plot if end of current orbit or end of file:
 */
	  if (sltp-rslt >= 1. || ir == nrec-1) {
	    norb++;
	    if (norb >= iborb && norb <= ieorb) {
	      irb = ir-nporb+1;
	      ire = ir;
	      if (ir == nrec-1) {
	        nporb++;
	        if (nporb > MXPORB) {
		  printf("Too many points this orbit: nporb=%d MXPORB=%d\n",
		    nporb,MXPORB);
		  clsgks_();
		  exit(0);
	        }
	        defdat(itime,rlat,rlon,rsza,rslt,un,vn,dat+nporb-1);
	        if (dogcm) defgcm(dat+nporb-1, gcm+nporb-1);
	        irb = ir-nporb+2;
	        ire = ir+1;
	      }
	      printf("\n");
	      printf("Plot orbit %2d nporb=%3d ir=%3d to %3d\n",
	        norb,nporb,irb,ire); 
	      pltorb(nporb,dat,dogcm,histname,gcm);
	      if (norb == ieorb) {
		printf("Exit: norb = ieorb = %d\n",norb);
		exit(0);
	      }
	    } else printf("Skipping orbit %2d (iborb,ieorb=%2d %2d)\n",
	             norb,iborb,ieorb);
	    nporb = 0; sltp = 0.;
/*
 * Define data for 1st point of next orbit (just read):

   9/18/92: binned data files rcvd from dwu have 1st slt of each orbit
   apparently out of order, so skip 1st point. e.g., last point of 
   1st orbit (file hrdi154.dat) is lt=16.25; then next 2 point are 
   lt=14.16 and 7.61, thence increasing from 7.61. So we skip 14.16.

	    if (ir != nrec-1) {
	      nporb++;
	      defdat(itime,rlat,rlon,rsza,rslt,un,vn,dat+nporb-1);
	      if (dogcm) defgcm(dat+nporb-1, gcm+nporb-1);
	    }
 */

/*
 * If not done with orbit or file, update arrays for current orbit:
 */
	  } else {
	    nporb++;
	    defdat(itime,rlat,rlon,rsza,rslt,un,vn,dat+(nporb-1));
	    if (dogcm) defgcm(dat+nporb-1, gcm+nporb-1);
	    sltp = rslt;
	  }
/*	  if ((ir+1) % 100 == 0) printdat(ir+1,dat+ir); */
	}
	clsgks_();
}