//////////////////////////////////////////////////////////////////////////
//
//       (C) Copyright 2003 by GATS, Inc.
//           11864 Canon Blvd., STE 101
//           Newport News, VA 23606     
//
//   All Rights Reserved. No part of this software or publication may be
//   reproduced, stored in a retrieval system, or transmitted, in any form
//   or by any means, electronic, mechanical, photocopying, recording, or
//   otherwise without the prior written permission of GATS, Inc.
// 
//////////////////////////////////////////////////////////////////////////
//
// Class: NMC_FileReader
//
// Filename:  NMC_FileReader.cpp
//
// Description:  Implementation file for NMC_FileReader class.
//               Refer to header file for more details.
//
// Author:      Y.Wang
//              (757)873-5920
//              <ywang@gats-inc.com>
//////////////////////////////////////////////////////////////////////////

#include <exception>
#include <string>
#include <strstream>
#include <fstream>
#include <iostream>
#include <cstring>
//#include <stream.h>
//#include <netcdf-3/netcdf.h>
#include <netcdf.h>
#include <math.h>
#include "NMC_FileReader.h"
#include "GeneralConsts.h"
#define ABOVE10mbar
#ifndef GREF
#define GREF 9.8
#endif

 
//#undef ABOVE10mbar

using namespace std;

///////////////////////////////////////////////////////////////////////
//
// Method:  NMC_FileReader
//
// Description:  Default constructor.
//
//////////////////////////////////////////////////////////////////////////

NMC_FileReader::NMC_FileReader()
{
   BasicInitialize();
}

///////////////////////////////////////////////////////////////////////
//
// Method:  NMC_FileReader
//
// Description:  Class constructor.  Gives the class the name of the file
//               to be read.
//
// Inputs:  the_filename - The name of the file to be read.
//
//////////////////////////////////////////////////////////////////////////

NMC_FileReader::NMC_FileReader(const std::string &the_filename)
                                 throw(const char*, exception)
{
   //cout<<"Input File Name: "<<the_filename<<"\n";
   BasicInitialize();
   //cout<<"after BasicInitialize!!!!!!!\n";
   SetUp(the_filename);
   //cout<<"after SetUp!!!!!!!\n";
   
}

///////////////////////////////////////////////////////////////////////
//
// Method:  ~NMC_FileReader
//
// Description:  Destructor.
//
//////////////////////////////////////////////////////////////////////////

NMC_FileReader::~NMC_FileReader()
{
   if(file_is_open)
   {
      nc_close(fd);

      file_is_open = false;
   }
}


///////////////////////////////////////////////////////////////////////
//
// Method:  SetUp
//
// Description:  Does everything that needs to be done before calling 
//               any of the other non-constructor/destructor methods.
//
// Inputs:  the_filename - The name of the NMC file to be read.
//
//////////////////////////////////////////////////////////////////////////

void NMC_FileReader::SetUp(const std::string &the_filename)
                            throw(const char*, exception)
{
   //cout<<"filename in SetUp"<<the_filename<<"\n";

   /*if (setup_attempted == true)
   {
      cerr<<"ERROR:  Attempt to re-set-up the NMC file reader."<<endl;
      throw("ERROR:  Attempt to re-set-up the NMC file reader.");
   }
   */
   setup_attempted = true;

   filename = the_filename;
   //cout<<"filename in SetUp"<<filename<<"\n";
   //cout<<"filename.c_str() in SetUp"<<filename.c_str()<<"\n";

   // Open the file.

   if (NC_NOERR != nc_open(filename.c_str(),
                           NC_NOWRITE,
                           &fd))
   {
      cerr << "ERROR:  Unable to open file for reading " << endl
                    << "        " << filename.c_str() << endl;

      throw("Stop in NMC file reader");
   }

   file_is_open = true;

   setup_succeeded = true;
   return;
}

///////////////////////////////////////////////////////////////////////
//
// Method:  BasicInitialize
//
// Description:  Initializes some data members of the class that need
//               initialization.  Done in this method so that the same
//               initialization can be done for each of the constructors.
//
//////////////////////////////////////////////////////////////////////////

void NMC_FileReader::BasicInitialize()
{
   setup_attempted = false;
   setup_succeeded = false;
   file_is_open = false;
   setup_succeeded = false;
   setup_attempted = false;
   return;
}

///////////////////////////////////////////////////////////////////////
//
// Method:  getStatus
//
// Description:  Returns 1 if setup succeeded
//
//////////////////////////////////////////////////////////////////////////

int NMC_FileReader::getStatus()
{
   if(setup_succeeded) return 1;
   return 0;
}




///////////////////////////////////////////////////////////////////////////
// Method: Read_netCDF_NMC_data
// Reads the data from the given netCDF file into memory
/////////////////////////////////////////////////////////////////////////////

void NMC_FileReader::Read_netCDF_NMC_data(NMCdata &nmcdata)
                        throw(const char*, exception)
{

   long  start[MAX_NC_VARS];
   long  end[MAX_NC_VARS];
   // Initialize start so all elements are zero.
   bzero(start, MAX_NC_VARS * sizeof(long));

   // Read in the variables individually.
   // Read: lat.
   end[0] = NMCLATdim;
   ncvarget(fd, 0, start, end, nmcdata.lat);
   //cout<<"nmcdata.lat[0]="<<nmcdata.lat[0]<<"\n";
 
   // Read: lon.
   end[0] = NMCLONdim;
   ncvarget(fd, 1, start, end, nmcdata.lon);
   //cout<<"nmcdata.lon[0]="<<nmcdata.lon[0]<<"\n";
   
   // Read: level.
   end[0] = NMCLEVELdim;
   ncvarget(fd, 2, start, end, nmcdata.level);
   //cout<<"nmcdata.level[0]="<<nmcdata.level[0]<<"\n";
#ifdef ABOVE10mbar

   // Read: hgt.
   end[0] = NMCLEVELdim;   end[1] = NMCLATdim;   end[2] = NMCLONdim;
   ncvarget(fd, 3, start, end, nmcdata.hgt);
   
   // Read: temp.
   end[0] = NMCLEVELdim;   end[1] = NMCLATdim;   end[2] = NMCLONdim;
   ncvarget(fd, 4, start, end, nmcdata.temp);
   
#else

   // Read: spres.
   end[0] = NMCLATdim;   end[1] = NMCLONdim;
   ncvarget(fd, 3, start, end, nmcdata.spres);
   //cout<<"nmcdata.spres[0]="<<nmcdata.spres[0][0]<<"\n";

   // Read: shgt.
   end[0] = NMCLATdim;   end[1] = NMCLONdim;
   ncvarget(fd, 4, start, end, nmcdata.shgt);
   //cout<<"nmcdata.shgt[0][0]="<<nmcdata.shgt[0][0]<<"\n";

   // Read: hgt.
   end[0] = NMCLEVELdim;   end[1] = NMCLATdim;   end[2] = NMCLONdim;
   ncvarget(fd, 5, start, end, nmcdata.hgt);
   
   // Read: temp.
   end[0] = NMCLEVELdim;   end[1] = NMCLATdim;   end[2] = NMCLONdim;
   ncvarget(fd, 6, start, end, nmcdata.temp);

#endif

   // Close the netCDF file.
   ncclose(fd);
   file_is_open = false;
}



void NMC_FileReader::Read_netCDF_NMC_data(float lat,float lon,
                             float *press_nmc,float *temp_nmc, float *gheight, double *altz_tp)
			     
                        throw(const char*, exception)
{
  long  start[MAX_NC_VARS];
  long  count[MAX_NC_VARS];
  bzero(start, MAX_NC_VARS * sizeof(long));
  int latReq= int(lat+0.5);  //find the closest integer
  int lonReq = int(lon+0.5);
  if (lonReq == 360) lonReq =0;
  int idxLat = (90-latReq);
  int idxLon = lonReq;
   
   // Read: level.
   count[0] = NMCLEVELdim;
   start[0]=0;
   ncvarget(fd, 2, start, count, press_nmc);
  
   // Read: geopotential height at requested latitude and longitude.
   
   start[0]=0; start[1]=idxLat; start[2]=idxLon;
   count[0] = NMCLEVELdim;   count[1] = 1;   count[2] = 1;
#ifdef ABOVE10mbar
   if(NC_NOERR != ncvarget(fd, 3, start, count, gheight)){
     cout<<"Error in reading gheight from NMC at lat "<<latReq<<" and lon at "
     <<lonReq<<endl;
#else
   if(NC_NOERR != ncvarget(fd, 5, start, count, gheight)){
     cout<<"Error in reading gheight from NMC at lat "<<latReq<<" and lon at "
     <<lonReq<<endl;
#endif     
   }  
   
   // Read: temp.
#ifdef ABOVE10mbar
   ncvarget(fd, 4, start, count, temp_nmc);
#else   
   ncvarget(fd, 6, start, count, temp_nmc);
#endif   

   convert_hgt2ztp(gheight, altz_tp, lat);
} 


/*void NMC_FileReader::convert_hgt2ztp(float* nmcGht, float* nmcZTP){
        for(int i=0; i<NMCLEVELdim; i++)
                nmcZTP[i] = nmcGht[i];
        return;
}*/
void NMC_FileReader::convert_hgt2ztp(float* hgt, double* ztp){
        for(int i=0; i<NMCLEVELdim; i++)
		ztp[i] = hgt[i];
        return;
}

void NMC_FileReader::convert_hgt2ztp(float* hgt, double* ztp, float lat){
	double g0, gEnd, gMean;
	g0 = calculate_gravity(0.0, lat);
        for(int i=0; i<NMCLEVELdim; i++){
                gEnd = calculate_gravity(hgt[i], lat);
		gMean = (g0 + gEnd) / 2.0;
		ztp[i] = hgt[i]*GREF/gMean;
	}
        return;
}

double NMC_FileReader::calculate_gravity(float alt, float lat){
	double latRadx2 = 2*lat*RADperDEGREE;
	double gSeaLevel = 9.80160 * (1.0 - 2.6373e-3 * cos(latRadx2) + 5.9e-6 * cos(latRadx2) * cos(latRadx2));		
	double gDiff = (3.085462e-4 + 2.27e-7*cos(latRadx2))*alt/100.0;
	return gSeaLevel - gDiff;

}