#include "Level1Data.h"
#include "GATS_Utilities.hpp"
#include "F77Radtran.h"

#include "fortranfunctions.h"

#include <iostream>
#include <cmath>
#include <cstring>


using namespace GATS_Utilities;


extern "C" void get_waccm_climatology ( int* usewaccm, char* waccm_dir, float* pr,float* zr,  float **q,
	int* igs,  int* nr,   int*  ng)
{

	static int c__1 = 1;
	if( *usewaccm == 0) return;
	
	char waccm_file[90];
	int ntimes = 1, nlatitudes=46, nlevels=66;
        int mr_temp = nlevels+1;
	int mgm_temp =16; 
        int idg_temp[] = {700,200,300,100,400,600,800,1000,1300,2200,3300,3400,3700,5100,5200,6100} ;
	float qmix_temp[mgm_temp][mr_temp],zr_temp[mr_temp],pr_temp[mr_temp],
		tr_temp[mr_temp],co2_temp[mr_temp] ;
	

	int	idate1,idate2;
	int     date_out[1];                         // time
	float   latitude_out[46];                    // lat
	float   pres_out[66];                        // lev
	float   co2_out[1][46][66];                  // time, lat, lev
	float   o_out[1][46][66];                    // time, lat, lev
	float   o1d_out[1][46][66];                  // time, lat, lev
	float   o2_out[1][46][66];                   // time, lat, lev
	float   n2_out[1][46][66];                   // time, lat, lev
	float   o3_out[1][46][66];                   // time, lat, lev
	float   t_out[1][46][66];                    // time, lat, lev
	float   z_out[66];                    // lev
	float   h2o_out[1][46][66];
	float   n2o_out[1][46][66];
	float   ch4_out[1][46][66];
	float   no_out[1][46][66];
	float   no2_out[1][46][66];
	float   oh_out[1][46][66];
	float   ho2_out[1][46][66];
	float   cfcl3_out[1][46][66];
	float   cf2cl2_out[1][46][66];
	float   h_out[1][46][66];

	float 	wght;

	std::string dt = gLevel1Data.StartTime();
	struct tm timestruct={0};
	
	timestruct.tm_year = ConvertFromString<int>( dt.substr(0,4), std::dec) - 1900;
	timestruct.tm_mon = ConvertFromString<int>( dt.substr(5,2), std::dec) - 1;
	timestruct.tm_mday = ConvertFromString<int>( dt.substr(8,2), std::dec) ;

	mktime(&timestruct) ;

	int idate = (1900+ timestruct.tm_year)*1000 + 1 + timestruct.tm_yday;
	std::string ds_time = ConvertToString( idate );


	::get_waccmfile_interp__( ds_time.c_str() , waccm_dir, waccm_file, &idate1,  
		&idate2, &wght,
		ds_time.size(), strlen(waccm_dir), sizeof(waccm_file) );


	::read_waccm_interp ( waccm_file, &idate1, &idate2, date_out, z_out, 
	      	    t_out,          h2o_out,     co2_out,      o3_out, 
		    n2o_out,        ch4_out,     o_out,        o2_out,
                  no_out,         no2_out,     oh_out,       n2_out,
                  ho2_out,        cfcl3_out,   cf2cl2_out,   h_out,
                  o1d_out,        latitude_out, pres_out,    &wght );


	float co2_fix;
//	::get_co2_trend__(&idate, &co2_fix);   // no longer using this... Its all WACCM CO2

	float  
		re_temp, 
		lat_temp = gLevel1Data.Lat83();


	::waccm_atmstate_interp__( &ntimes, &nlatitudes, &nlevels, z_out,  
           pres_out,  &t_out[0][0][0], &h2o_out[0][0][0],  &co2_out[0][0][0], &o3_out[0][0][0], &n2o_out[0][0][0], &ch4_out[0][0][0],
	  &o_out[0][0][0], &o2_out[0][0][0], &no_out[0][0][0], &no2_out[0][0][0], &oh_out[0][0][0], &n2_out[0][0][0], &ho2_out[0][0][0],
          &cfcl3_out[0][0][0], &cf2cl2_out[0][0][0], &h_out[0][0][0], &o1d_out[0][0][0], &mr_temp, &mgm_temp, idg_temp , &re_temp, &lat_temp, 
		zr_temp,  pr_temp, tr_temp,&qmix_temp[0][0] , co2_temp , &co2_fix);

	
/*
        for(int i = 0; i < *nr; ++i) {
                std::cout << q[j][i] << "  " << zr[i] << std::endl;
        }
       std::cout << "&" << std::endl;
*/

// we interpolate the waccm to the input grid on pressure coordinate
//
	std::vector<float> LogPR(pr, pr+ *nr);
	std::transform(LogPR.begin(), LogPR.end(), LogPR.begin(), log );
	std::vector<float> LogPOUT(pr_temp, pr_temp+mr_temp);
	std::transform(LogPOUT.begin(), LogPOUT.end(), LogPOUT.begin(), log );

	reverse(LogPR.begin(), LogPR.end() );
	reverse(LogPOUT.begin(),  LogPOUT.end() );

	std::vector<float> A_( *nr) ;
	for(int i = 0; i< *ng; ++i) {

		int* ik = std::find(idg_temp,  idg_temp+mgm_temp , igs[i] );
		int j = std::distance( idg_temp,ik);

		if( j <  mgm_temp ) {
			std::vector<float> B_(&qmix_temp[j][0], &qmix_temp[j][0]+ mr_temp); 
			std::reverse(B_.begin(), B_.end() );		
			::kp_intprof_noxtrp__( nr, &mr_temp, &c__1, nr, &LogPOUT[0], &LogPR[0], 
			&B_[0], &A_[0]  );
     
           std::reverse(A_.begin(), A_.end());
		   std::copy(A_.begin(), A_.end(), &q[i][0] );

		}
     }

}