//
//  $Id: CO2ProfileSimulation.cpp,v 1.7 2008/01/08 02:28:50 deaver Exp deaver $
//-----------------------------------------------------------------------
//
//            (c) Copyright 2006 by GATS, Inc.,
//     11864 Canon Blvd, Suite 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.
//
//-----------------------------------------------------------------------
//
// Module:	CO2ProfileSimulation.cpp
//
// Author:	John Burton	<john@gats-inc.com>
//
// Date:  	Thu May  4 17:01:31 2006
//
//-----------------------------------------------------------------------
//
// Modification History:
//
//     $Log: CO2ProfileSimulation.cpp,v $
//     Revision 1.7  2008/01/08 02:28:50  deaver
//     Added code so that it will update the CO2 Profile with a modified WACCM profile
//     as is done for SABER.  Added a UseWaccmCO2 flag and a WaccmDir entry to the control file
//
//     Revision 1.6  2007/11/27 20:56:19  deaver
//     Renamed F77Linepak.h to F77Radtran.h
//
//     Revision 1.5  2007/11/27 15:51:50  deaver
//     renamed GeneralUtilities.h to F77Functions.h
//     made atmfile a static string
//
//     Revision 1.4  2007/07/30 15:36:21  deaver
//     Uses readfasatm to get the climatology VMR
//     added an eventvar to ini file for the calculated refraction angle
//
//     Revision 1.3  2007/07/25 01:12:58  deaver
//     added a fixed CO2 Profile  TODO need to output refraction angle in eventVAR
//
//     Revision 1.2  2007/07/20 15:15:31  deaver
//     Modifications now needs testing
//
//     Revision 1.1  2007/07/10 03:46:55  deaver
//     Initial Checkin
//
//
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

#include "CMarkup.h"
#include "ConfigFile.h"
#include "Event.h"
#include "CO2ProfileSimulation.h"
#include "LoadLinesFromDB.h"
#include "F77Radtran.h"
#include "F77Functions.h"
#include "FORTRANarray.hpp"
#include "GATS_Utilities.hpp"
#include "DatabaseHandles.h"

#include <cmath>
#include <ctime>
//#include<iostream>
//#include <numeric>

//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//
#define THIS_MODULE	"CO2ProfileSimulation"
extern "C" {

void read_waccm_interp ( char *INFile, int *DATE1, int *DATE2, int *date_out, float *z_out,
                  float tk_out[][46][66],  float h2o_out[][46][66],    float co2_out[][46][66],    float o3_out[][46][66],
                  float n2o_out[][46][66], float ch4_out[][46][66],    float o_out[][46][66],      float o2_out[][46][66],
                  float no_out[][46][66],  float no2_out[][46][66],    float oh_out[][46][66],     float n2_out[][46][66],
                  float ho2_out[][46][66], float cfcl3_out[][46][66],  float cf2cl2_out[][46][66], float h_out[][46][66],
                  float o1d_out[][46][66], float *latitude_out, float *p_out,        float *wght);

};


//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//

// there are several global routines and structures defined in FORTRAN code


//
//-----------------------------------------------------------------------
// Utility Routines:
//-----------------------------------------------------------------------
//

using namespace GATS_Utilities;
using namespace gatsDBpp;

int CO2ProfileSimulation(Event& L0, Event& L1, Event& Tmp, Event& SD, ConfigFile& cf)
{
	const static int msiz = 30000000;
	const static int mfilt = 1000;
	const static int mf = 1;
	const static int mgm = 14;

	char nam[mgm*8],atmos[2], waccm_file[100];
	float filt[mfilt];
	int idf[mf];
	FORTRANarray<double> crss(msiz);
		
	int mcel, ncel, ner, mgas, ngas, mray, 
		nray,  nr,  mr, ngm, mpout, mcout,
		mcntrl,nhvy, nl;

	static int debug = 0;   		// 0 is a false (this is the debug flag)
	static int nfilt=1;			// Number of filters to read in
	static int irf =  1;			// 1 is true always do refraction
	static int formatted = 1;			// file is opened formattted
	static int ihydro = 0;			// Do not call hydro in readfasatm
	static int ppmv = 1;			// convert vmrs from ppmv to mixing ratio in readfasatm
	static int readzpt=0;			// Do not read in a z,p,t profile in readfasatm
	static int iplat = 0;  		// observer in space
	static int nsr = 1 ;  		// always start at layer 1;
	static int nset = 1;			//only one set of cross sections
	static int imr = 2;          	// multi Ray
	static int segment=0;			// don't segment here
	static int icalc = 1;  		// transmission only
	static int iopt = 1;			// limtran method of integration
	static int ic1 = -1;		//beginning cell for integration
	static int ic2 = 1;		//ending cell for integration
	static int nfpmx = 750;		// max number points per filter
	static int mcmix = 1;			// max number of parameters for line mixing
	static int mlmix = 5;         // Max number of lines allowd to have mixing
	static int mbnd = 5830;		// Max number of points for band model
	static int mcof = 6;			// Max number of coefficients for each band model
	static int mset = 1;          // only one cross section set per gas;
	static int lfilt=10;			// Lun of filter file.
	static int lhvy=11;			// LUN of the cross section file
	static int lout = 12;			// LUN of the debug file
	static int latmos = 20;		// Lun used by the readfasatm subroutine
	static float solar=0.0;		// No solar weighting applied in Lp_filter
	static float albedo=0.0;
	static float dvel=0.0;		// no doppler shift
	static int ifl=1;			// tells cfilter to convolve first filter only;
	static double resmax =0.;		// Linepak picks its own resolution
	static float odmin = -1.0;	// if negative Linepak will use its default value 
	static float resfc = 0.3;  	// if Negative Linepak will use its default value   
	static  bool firstCall = true;
	static bool useWaccm= true;
	static ConfigFile cfLocal;
	static std::vector<int> SOFIE_bands;
	static float wnmax;
	static std::string atmfile;
	static std::string waccm_dir;
	static int c__1 = 1;

		
	FORTRANarray<float> rlt, rut, rl, ru, zt, pt, 
		tt, gr, qmix, qg, pb, tb, tbs, tedge, gasb, pmass,
		q, qgr;
	std::vector<float> zr, pr, tr, tgr, pout, coutA, cntrl,
		hvy,	afgl;
	std::vector<double> transmittance, refractionAng;
	std::vector<int> igas, igm;
//	float lat83;
	float  re, tearth, wave,erad,elat;
	int ir, isc, done, initCross, initIntegrate;
	int bandNo;
	double vs1, vs2, resfilt, avtran, avrad, filtti;
	std::string sectionName, TranEventVarName, AltEventVarName,
		RefracEventVarName,
		lineTableName;

	DatabaseHandles* handles;
	GATS_DB  *dbConnection;



//  executable code starts here  
// initialization stuff done once
	handles = GetDatabaseHandles();
	if( firstCall  ) {
		std::string Filename = cf.GetStr(THIS_MODULE, "IniFile");
		cfLocal.open(Filename);
		sectionName="BandSimulations";

		SOFIE_bands = cfLocal.GetIntValues(sectionName,"BANDS");

		Filename = cfLocal.GetStr(sectionName,"FilterFile");
		::openfile_( &formatted,&lfilt, Filename.c_str(),  
			(int)Filename.size() );
			
		Filename = cfLocal.GetStr(sectionName,"BandModelFile");
		::openfile_( &formatted,&lhvy, Filename.c_str(),  
			(int)Filename.size() );
			
		Filename = cfLocal.GetStr(sectionName,"DebugFile");
		::openfile_( &formatted,&lout, Filename.c_str(),  
			(int)Filename.size() );

		atmfile = cfLocal.GetStr(sectionName,"AtmosDir");
		atmfile += " ";   // important that it has a blank at the end

		if( cfLocal.ValidEntry(sectionName , "UseWaccmCO2")  && 
			cfLocal.GetInt(sectionName , "UseWaccmCO2") != 0 ) {
			waccm_dir = cfLocal.GetStr(sectionName,"WaccmDir") + " ";
			useWaccm = true;
		} else {
			useWaccm = false;
		}

		Filename = cfLocal.GetStr(sectionName,"DatabaseInfoFile");
		CMarkup* xml = CMarkup::LoadXMLFile(Filename);
		xml->InsideRoot();
		xml->FindElem("databaseConnections");
		handles->Init(xml);
		delete xml;

		wnmax = (cfLocal.ValidEntry(sectionName , "WNMAX") )? 
			cfLocal.GetReal(sectionName , "WNMAX" ) :  50.0;
			
// Call  LP_inita  so we can PROCEED
		mcel=mray=2;
		mgas=1;
		rlt= FORTRANarray<float>(mcel,mgas);      	
		rut = rlt;
		ru = FORTRANarray<float>(mcel,mray,mgas); 		
		rl = ru;
		pout.resize(mcel);
		::lp_inita__(&debug, &mcel, &mray, &mgas, rlt.getArray(),
			rut.getArray(), rl.getArray(), ru.getArray(),
			(int*)&pout[0]  )  ;

		firstCall = false;
	}


// Pull data from the L1 Event object.
  	EventVar<double> Z, P, T,  radCurvature ;
	L1.getEventVar("Z_Merged", Z);  //  merged NCEP and MSIS z,p,t profile in hydrostatics
	L1.getEventVar("P_Merged", P);
	L1.getEventVar("T_Merged", T);
	L1.getEventVar("L1_CurvatureRadius", radCurvature);
	re = (float)radCurvature[0];
//	lat83 = (float)L1.getLatitude83km() ;
	time_t unixTime = (time_t)L0.getEventStartTime() ;
	struct tm* timeinfo = gmtime(&unixTime) ;
//	isc = L1.getEventNumber();
		
// make floats from doubles
	zr = std::vector<float>(&Z[0], &Z[0]+Z.size() );
	pr = std::vector<float>(&P[0], &P[0]+P.size() );
	tr = std::vector<float>(&T[0], &T[0]+T.size() );	
	nr = mr=(int)zr.size();
// no temperature Gradients	
	tgr = std::vector<float>(2*mr, 0.);

// Load in the mixing ratios and extinctions from climatology file	
	igm = std::vector<int>(mgm);
	q=FORTRANarray<float>(mr,mgm);
	qgr=FORTRANarray<float>(2,mr,mgm);

	::readfasatm2_( &latmos,  atmos ,  &erad , &elat ,
        	&mr,    &nr,    &ngm,      nam,     
        	&igm[0],  &zr[0], &pr[0],  &tr[0] , 
        	&tgr[0],  q.getArray(), qgr.getArray() , &ihydro,
       		atmfile.c_str(),  &ppmv, &readzpt, sizeof(atmos), 8, 
		atmfile.size() );
// We need to replace the CO2 with the WACCM  profile as per SABER

if(useWaccm) {
	int ntimes = 1, nlatitudes=46, nlevels=66;
	float wght;
	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];


	int idate1 = (1900+ timeinfo->tm_year)*1000 + 1 + timeinfo->tm_yday, idate2;
	std::string ds_time = ConvertToString( idate1 );

	::get_waccmfile_interp__( ds_time.c_str() , waccm_dir.c_str(), waccm_file,
		&idate1, &idate2, &wght, ds_time.size(), waccm_dir.size(), 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);

	float 
		re_temp, 
		lat_temp =  (float)L1.getLatitude83km();

        int mgm_temp =16;
        int idg_temp[] = {700,200,300,100,400,600,800,1000,1300,2200,3300,3400,3700,5100,5200,6100} ;

        int mr_temp = nlevels+1;//, mgm_temp =16, idg_temp[mgm_temp];
	float qmix_temp[mgm_temp][mr_temp],zr_temp[mr_temp],pr_temp[mr_temp],
		tr_temp[mr_temp],co2_temp[mr_temp] ;

	::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);


	std::vector<float> LogPR(pr.begin(), pr.end() );
	for(unsigned j =0; j< LogPR.size(); ++j) {
		LogPR[j] = log(LogPR[j]);
//std::transform(LogPR.begin(), LogPR.end(), LogPR.begin(), log);
	}
	std::vector<float> LogPOUT(pr_temp, pr_temp+mr_temp);
	for( unsigned j=0; j< LogPOUT.size(); ++j) {
		LogPOUT[j] = log(LogPOUT[j] );
//	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< ngm; ++i) {

		int* ik = std::find(idg_temp,  idg_temp+mgm_temp , igm[i] );
		int j = std::distance( idg_temp,ik);
//std::cout << "this is j " << j << "  " << *ik << std::endl;
		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.getArray(1,i+1) );
		}
	}

}
// wasn't that icky???

/*
	int* ik = std::find(&igm[0],&igm[0]+ngm, 200);  //locate co2 id 200
	int j = std::distance(&igm[0],ik);
for(int i = 0; i < nr; ++i) {
	std::cout << q(i+1,j+1) << "  " << zr[i] << std::endl;
}
exit(23);
*/

// loop over the requested band numbers and calculate transmittance for each

	for(unsigned iband = 0; iband < SOFIE_bands.size() ; ++iband) {

		bandNo = SOFIE_bands[iband];
		sectionName = "Band_" + ConvertToString(bandNo);
		std::vector<float> za = cfLocal.GetRealValues(sectionName,"ImpactZ");
		TranEventVarName = cfLocal.GetStr(sectionName, "TranEventVarName");
		AltEventVarName = cfLocal.GetStr(sectionName, "AltEventVarName");
		RefracEventVarName = cfLocal.GetStr(sectionName, "RefracEventVarName");
		mcel=ncel=ner=nray=mray= (int)za.size() ;
    		int it2 = mcel * mcel;
    		int it1 = it2 + mcel;
    		mpout = mcel * 11 + 3 + it1 * 12 + it2 * 3;	
   		mcout = mcel * 11 + it1 * 6 + it2 * 3; 			
		transmittance.clear();
		refractionAng.clear();

// Allocate some arrays
		zt=FORTRANarray<float>(mcel);
		pt=FORTRANarray<float>(mcel);
		tt=FORTRANarray<float>(mcel);
		gr=FORTRANarray<float>(2,mcel);
		tb=FORTRANarray<float>(mcel);
		pb=FORTRANarray<float>(mcel);
		pout.resize(mpout);
		coutA.resize(mcout);
				
// read the filter data file for this band.
		idf[0] = bandNo;

		::lp_filter__(&lfilt, &mf, &mfilt, 
			&nfilt, &nfpmx, &solar, &vs1, 
			&vs2, idf, &resfilt,  filt);

//apply a temperature shift spread here when ready

// get the filter wavenumber limits  return as vs1 and vs2.
		::filterlimits_( &vs1, &vs2, &idf[0], &nfilt, filt);
// refraction wavenumber  may just hard code to something
		wave = (float)(vs2 + vs1)/2.0;

// Get the absorbers used in this band		
		igas = cfLocal.GetIntValues(sectionName,"GasIds");
		ngas=mgas=(int)igas.size();

// Allocate/initialize more arrays
		rlt = FORTRANarray<float>(mcel,mgas);
		std::fill(rlt.getArray(), rlt.getArray()+rlt.size(), 1.0) ;
		rut = rlt;
		ru = FORTRANarray<float>(mcel,mray,mgas); 
		std::fill(ru.getArray(), ru.getArray()+ru.size(), 1.0)	;	
		rl = ru;
		qmix = FORTRANarray<float>(mcel,mgas);
		qg = FORTRANarray<float>(2,mcel,mgas);
		tbs = FORTRANarray<float>(mcel,mray,2,mgas);
		tedge=FORTRANarray<float>(mcel,mray,2,mgas);
		gasb = FORTRANarray<float>(mgas,mcel);
		pmass = FORTRANarray<float>(mgas,mcel,mray,2);
		
		FORTRANarray<float> amass = pmass;
		std::fill(amass.getArray(), amass.getArray()+amass.size(),0.);

		mcntrl = (mgas << 1) + 6 + mcel * (mgas * 13 + 5 + mgas * mray);
		cntrl.resize(mcntrl);		
		
// Perform the Mass Calculations
		isc = bandNo;
		::lp_mass__( &mr, &mcel, &mgas,
			&mgm, &mpout, 
			&ncel, &nr, &irf, 
			&nsr, &ner, &ngas, &igas[0],
			&ngm, &igm[0], &isc, &re, &wave,
			&iplat, &za[0], &zr[0], &pr[0], &tr[0], &tgr[0],
			q.getArray(), qgr.getArray(), zt.getArray(), pt.getArray(),
			tt.getArray(), gr.getArray() , 
			&tearth, qmix.getArray(), qg.getArray(), &pout[0]   );

// Calculate the total Refraction angle along the limb path for each for the ray. (radians)

		for(ir=nsr; ir<= ner; ++ir ) {
			refractionAng.push_back( ::refraction_( &ir, &ic1, &ic2, &pout[0] ) );	
		}

// mass apportioning
		::lp_celltb__( &mcel, &mray, &mgas, 
			&mpout, &mcout, &ncel, &nray,  
			&nset, &imr, &nsr, &ner, &ngas, 
			&igas[0], &re, pt.getArray(), tt.getArray(), rlt.getArray(), 
			rut.getArray(), qmix.getArray(), qg.getArray(), tb.getArray(), 
			tbs.getArray(), tedge.getArray(), rl.getArray(), 
			ru.getArray(), pb.getArray(), gasb.getArray(), 
			pmass.getArray(), &pout[0], &coutA[0]  );

// Get the line parameters and band model data

 		nhvy=4+ngas*(2+ 11 * mset+ mcof * mbnd * mset);
		hvy.resize(nhvy);
		::lp_hvyset__( &lout, &lhvy, &debug,
	 		&mcof, &mbnd, &mset, &ngas, &igas[0], 
	 		&vs1, &vs2,  &hvy[0]  );
// NOTE: afgl is allocated inside this function	 
		dbConnection = handles->Get("Linelist").get();	
		lineTableName = "hitran2004" ;	
		nl = LoadLinesFromDB(  igas, vs1, vs2, wnmax, mcmix, mlmix, 
			afgl,  dbConnection ,lineTableName);

// calculate the cross sections

		initCross = 1;

		::lp_crossm__(  
			&lout,  &iopt, &icalc,
	 		&mcel, &mray, &mgas, &mcntrl, &msiz, 
	 		&nl, &mcmix, &ncel, &nray, 
			&ngas, &igas[0], &vs1,  &vs2, 
			&wnmax, &resfc, &resmax, &odmin, &nsr, 
			&ner, pb.getArray(), tb.getArray(), gasb.getArray(), 
			amass.getArray(), &segment, &cntrl[0], &afgl[0], &hvy[0], 
			crss.getArray(), &initCross, &done, &nset, &imr, 
			&mpout,  &pout[0], rl.getArray(), ru.getArray()  );

// 	   Integrate the rays and compute the band averaged transmission

		for( ir=nsr; ir <= ner; ++ir) {
			initIntegrate=1;

			::lp_cingrate__( &lout, &icalc, &iopt, &mcel, 
				&mray, &mgas, &ncel, &nray, &initIntegrate, 
				&ir, &ic1, &ic2,
				&ngas, &igas[0], &vs1, &vs2, &albedo, 
				&tearth, tbs.getArray(), 
				tedge.getArray(), pmass.getArray(), &pout[0],
				&cntrl[0], crss.getArray() )	;

			::lp_cfilter__( &icalc, &nfilt,
				&ifl,  &mf, &nfpmx, &dvel, &avtran, 
				&avrad, &filtti, &cntrl[0], crss.getArray(), filt );
			
			transmittance.push_back( avtran / filtti );
		}
// make the EventVar's for this Band
	
		EventVar<double>  ev1( TranEventVarName , 
			sectionName+" Simulated Transmittance", 
			&transmittance[0], transmittance.size() );
// altregistration needs EventVar<double> for everything 
		std::vector<double> zdub(za.begin(), za.end() );
		EventVar<double>  ev2( AltEventVarName , 
			sectionName+" Impact Altitudes", "km",
			&zdub[0], zdub.size() );
	
		L1.addEventVar( ev1 ); 
		L1.addEventVar( ev2 );

	EventVar<double> ev3(RefracEventVarName,
		sectionName+" Simulated Refraction Angle","radians",
		&refractionAng[0], refractionAng.size() );
		L1.addEventVar(ev3);

	}  // go to the next band 

  return 0;
}