/* ncREADday7.c 
 *
 *  Revised:	 29jun99 - nanr
 *
 *  INPUT FILE:  geog3.nc 
 *		 <file pathnames>
 *      OUTPUT:  output file
 *	
 * USE TO READ:  TCLIMATE (daily)
 *		 TSCENARIO (daily)
 *		 TSPIN (daily)
 *
 *      TO RUN:  call the executable from within the
 *               directory containing the netCDF 
 *               file you wish to read.
 *
 *  This interactive C program reads 
 *  VEMAP2 netCDF files, and writes an output 
 *  file in your current directory
 *  containing the specified variable.
 *
 *  Note you must customize your pathname for geog3.nc.
 * ----------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* >>>>>> VEMAP USER's ----- Customize pathname here ---- */
#include "/contrib/netcdf-3.4/include/netcdf.h"

#define HISRECORD 36159
#define HADRECORD 38716
#define CCCRECORD 39081
#define TSPINRECORD 36500
#define ONEYEAR 365


void grid(int gridpt, int *row, int *col, float *olat, float *olon);
void background(int row,int col,int gridpt,int *startpos);
int leapyear(int);

int main(argc, argv)
int argc;
char *argv[];
{

   FILE *fp_out=NULL;

/* misc ids */
  char BUFFER[512];

/* Netcdf file id */
  int ncid;	

/* variable ids */
  int varid;
  char gp[5];
  char var_name[4];
  char name[512];
  char file_name[512];
  char units[20];

/* corners and edge lengths */
  long int count[2];
  long int start[2];

/* data variables */
  float var[CCCRECORD];
  float slice[ONEYEAR+1];

/* attribute vectors */
  char *vr_val;
  int vr_len;
  nc_type vr_type;


/* Grid variables */
  int gridpt;		/* For VEMAP US Grid 1 thru 5520 */
  int startpos;		/* Gridpt record position */
  int row;		/* Row number */
  int col;		/* Column number */
  int nc_err;		/* netCDF error */
  float olat;		/* cell center latitude */
  float olon;		/* cell center longitude */

  int i=0, lyr=0, nyrs=0, recs=0, y=0, fflag=0, flag=0, year=0, yr=0, ctr=0, ye=0;

/* --- Command line instructions --- */
  fprintf(stderr,"Enter <variable> (tx,tn,ppt,sr,vp,irr,rh)\n"); 
  scanf("%s",var_name);

  fprintf(stderr,"Enter <filename> (e.g., srdyT1.nc) \n"); 
  scanf("%s",file_name);

  *(file_name + strlen(file_name)) = '\0';

  fprintf(stderr,"Enter <grid cell> (1-5520)\n");
  scanf("%s",gp);

  fprintf(stderr,"Enter POR fflag:\n 1 = 1895-1993 (TCLIMATE)\n 2 = 1994-2099 (Hadley)\n 3 = 1994-2100 (CCC)\n 4 = 1-100 (TSPIN)\n 5 = 1895-1993 (HADLEY Historical)\n 6 = 1900-1993 (CCC Historical)\n");
    scanf("%d",&fflag);

  fprintf(stderr,"Do you want the full Monty? <1=yes,0=no> \n");
  scanf("%d",&flag);

  if(!flag) {
      /* TCLIMATE */
      if(fflag == 1 )
        fprintf(stderr,"Enter <year> (1895-1993) \n");
      /* Hadley */
      if(fflag == 2)
        fprintf(stderr,"Enter <year> (1994-2099) \n");
      /* CCC */
      if(fflag == 3)
        fprintf(stderr,"Enter <year> (1994-2100) \n");
      /* TSPIN */
      if(fflag == 4)
        fprintf(stderr,"Enter <year> (1-100) \n");

      scanf("%d",&year);
      if((fflag == 1  ) && (year > 1993 || year < 1895)) {
        fprintf(stderr,"Period of Record is (1895-1993) \n");
        fprintf(stderr,"Re-enter <year> \n");
      	scanf("%d",&year);
      }
      if(fflag == 2 && (year > 2099 || year < 1994)) {
        fprintf(stderr,"HADLEY Period of Record is (1994-2099) \n");
        fprintf(stderr,"Re-enter <year> \n");
      	scanf("%d",&year);
      }
      if(fflag == 3 && (year > 2100 || year < 1994)) {
        fprintf(stderr,"CCC Period of Record is (1994-2100) \n");
        fprintf(stderr,"Re-enter <year> \n");
      	scanf("%d",&year);
      }
  }


  /* -------- End Command Line ------------*/

  gridpt = atoi(gp);

/* --- Translate gridpoint --- */ 
/*    (1) latitude/longitude   */
/*    (2) row and column       */

  grid(gridpt,&row,&col,&olat,&olon);

/* ---------Screen for background point ------*/


  background(row,col,gridpt,&startpos);

  if(startpos == -9999) {
             printf("%5d is a background point %d \n",gridpt,startpos);
             exit(0);  
  }
  else {


	strcat(name,var_name); strcat(name,"DAY");
	fp_out = fopen((char *)strcat(strcpy(BUFFER,name),gp), "w");
	ncid = ncopen(file_name,NC_NOWRITE);
    		if (ncid == -1) printf("INPUT OPEN FAILED\n");

  	fprintf(fp_out,"          file = %s\n", file_name);

/* --- get variable ids from file --- */
	varid = ncvarid(ncid, var_name);
	if (varid == -1) printf("  VARID SCREWED\n");

	fprintf(fp_out," Variable name = %s\n", var_name);

/* --- retrieve variable units --- */
	ncattinq(ncid,varid,"units",&vr_type, &vr_len);
	vr_val = malloc((vr_len+1)*nctypelen(vr_type));
	ncattget(ncid,varid,"units",vr_val);
	vr_val[vr_len] = '\0';
	fprintf(fp_out,"         units = %s \n", vr_val);

/* --- retrieve variable scaling factor --- */
	ncattinq(ncid,varid,"scaled",&vr_type, &vr_len);
	vr_val = malloc((vr_len+1)*nctypelen(vr_type));
	ncattget(ncid,varid,"scaled",vr_val);
	vr_val[vr_len] = '\0';
	fprintf(fp_out,"scaling factor = %s \n", vr_val);

/* --- retrieve variable period_of_record --- */
        ncattinq(ncid,varid,"period_of_record",&vr_type, &vr_len);
        vr_val = malloc((vr_len+1)*nctypelen(vr_type));
        ncattget(ncid,varid,"period_of_record",vr_val);
        vr_val[vr_len] = '\0';
        fprintf(fp_out,"period_of_record = %s \n", vr_val);

/* --- set start position for gridpt --- */
/*     NOTE:  C and FORTRAN start positions differ */
	start[0] = startpos-1;
	start[1] = 0;
	count[0] = 1;

/* --- retrieve variable records --- */

        if(fflag == 1 ) {
                yr=1895; recs =  HISRECORD; nyrs =  99; lyr = 1993;
                count[1] = recs;
        }
        if(fflag == 2) {
                yr=1994; recs = HADRECORD; nyrs = 106; lyr = 2099;
                count[1] = recs;
        }
        if(fflag == 3) {
                yr=1994; recs = CCCRECORD; nyrs = 107; lyr = 2100;
                count[1] = recs;
        }
        if(fflag == 4) {
                yr=1; recs = TSPINRECORD; nyrs = 100; lyr = 100;
                count[1] = recs;
        }

	ncvarget(ncid, varid, start, count, var);

/* --- calculate desired POR from input year --- */

	if(!flag ) {
		/* calculate start index */
		if(fflag == 4) 
		{
		    y=0; 
		    for(ctr=yr;   ctr < year; ctr++) y += (ONEYEAR);
		    for(i=0; i<ONEYEAR; i++) { slice[i] = var[y]; y++; }
		}
		else 
		{
		    y=0; 
		    for(ctr=yr;   ctr < year; ctr++) y  += (ONEYEAR + leapyear(ctr)); 
		    printf("y = %d \n",y);
		    for(i=0; i<ONEYEAR+leapyear(year); i++) { slice[i] = var[y]; y++; }
		}
	}

/* --- Print gridpoint information --- */
	fprintf(fp_out,"      gridcell = %d\n", gridpt);
	fprintf(fp_out,"           row = %d\n", row);
	fprintf(fp_out,"        column = %d\n", col);
	fprintf(fp_out,"      latitude = %6.2f\n", olat);
	fprintf(fp_out,"     longitude = %6.2f\n", olon);
 	fprintf(fp_out," record number = %d\n\n", startpos);
        fprintf(fp_out," Records number from 1->3261.\n\n");
        fprintf(fp_out," Recall that record number does not \n");
        fprintf(fp_out," necessarily equal netCDF position index:\n\n");
        fprintf(fp_out," netCDF index: \n");
        fprintf(fp_out,"             C = %d \n",startpos-1);
        fprintf(fp_out,"       Fortran = %d\n\n",startpos);
        fprintf(fp_out,"-----------------------------------------\n");
	if(!(flag)) 
        fprintf(fp_out,"Year of Record = %d (daily) \n", year);
	else
        fprintf(fp_out,"Period of Record = %s (daily) \n",vr_val);
        fprintf(fp_out,"-----------------------------------------\n");

/* --- Print requested record --- */
	printf("Writing file:  %sDAY%d\n",var_name,gridpt);

	if(flag)  {
	    ctr=0;
	    for (i=0; i<recs;  i++) {
		   if(fflag == 4) { if(ctr == 365)  { ctr=0; yr++;} }
		else {
               	   if(leapyear(yr)) { if(ctr == 366)  { ctr=0; yr++;}}
               	   else             { if(ctr == 365)  { ctr=0; yr++;} }
		}
                fprintf(fp_out,"%4d %3d %8.2f\n",yr, ctr+1, var[i]);
                ctr++;
	    }
        }
	else    {
		if(fflag == 4) {
		   for (i=0; i<ONEYEAR; i++) 
 			fprintf(fp_out,"%4d %3d %8.2f\n",year,i+1,slice[i]);
		}
		else {
		   for (i=0; i<ONEYEAR+leapyear(year); i++) 
 			fprintf(fp_out,"%4d %3d %8.2f\n",year,i+1,slice[i]);
		}
	}

	nc_err = ncclose(ncid);
	exit(0);
   }
}


/*------------------------------------------------------------*/
void grid(int gridpt, int *row, int *col, float *olat, float *olon)

{
/* Program converts gridpoint to 2D matrix (row and col) */
/*  and returns gridcell center latitude and longitude.  */
/*  nanr 12-4-95 	           		         */

      int z;

      for (z=0; z<48; z++) {
            if ((gridpt > z*115) && (gridpt <= (z+1)*115)) {
                   *row = z+1;
                   *col = gridpt-z*115;
            }
      }

      *olat =  49.0 + (0.5 * (1 - *row)) - 0.25;
      *olon = 124.5 + (0.5 * (1 - *col)) - 0.25;

      return;
}
/*------------------------------------------------------------*/
void background(int row,int col,int gridpt,int *startpos)

{
         int backgrpt=0, datapt=0, sum=0;

	         int varid, ncid;
         int nrows=48, ncols=115;
         int irow=0, i=0;
         long int count[2] = {1, 115};
         long int start[2];
         short var[115];
         FILE *fopen(), *fp_in;
 
         if(!(ncid = ncopen("/fs/cgd/home0/nanr/tdaily/source/geog3.nc",NC_NOWRITE))) {
                printf("ERROR:  Can't find 'geog3.nc' file. \n");
                exit(0);
         }
 
        varid = ncvarid(ncid, "mask");
        while(irow < row) {
                start[0] = irow;
                start[1] = 0;
                count[0] = 1;
                count[1] = ncols;
                ncvarget(ncid, varid, start, count, var);
                if(irow<row-1) {
                        for(i=0; i<ncols; i++) {
                                if(var[i]) datapt++;
                                else       backgrpt++;
                        }
                }
                else if(irow==row-1) {
                        for(i=0; i<col; i++) {
                                if(var[i]) datapt++;
                                else       backgrpt++;
                        }
                }
                irow++;
         }

         sum = backgrpt + datapt;
         if(sum != gridpt) {
                printf("bad count %d\n",sum);
                exit(1);
         }
         if(!(var[col-1])) *startpos = -9999;
         else              *startpos = datapt;
         return;
}
/* Returns 1 if year is a leap year, 0 otherwise */
int leapyear(int year)
{
        return year%4?0:(year%100?1:(year%400?0:1));
}