#include "fortranfunctions.h"
#include "F77Radtran.h"

#include <vector>
#include <numeric>
#include <iostream>
#include <cassert>
#include <fstream>
#include <algorithm>

std::ofstream tempFile;
std::ofstream h2oFile;
std::ofstream co2File;

extern "C" void merge_into_clim( float* za,	float* ta,  float* pa,
	float**  qmix,	 int*   nl,	int*  igas,	int* ngas,	
	int* igdo,	int* ngdo,
	int* nch,	float *ztop,	float *zbot,	float* top_dz,	float* bot_dz,
	float* zr,	float* pr,	float* tr,	float** qmix_r,
	int* nr,	int* igs,	int* ngg,
	float* xlat,	float *zreg,	int* iscale_top,	int* iscale_bot,
	float* earth_rad )
{
	static int c__1=1;
/*
	static int once = 1;
	if(once == 1) {
     std::cout << " found true " << std::endl;
		tempFile.open("./tempRet.out");
		co2File.open("./co2Ret.out");
		h2oFile.open("./h2oRet.out");
		once = 0;
	} 
*/
//for(int il=0; il< *nr; ++il)
//	std::cout << tr[il] << "  " << zr[il] << std::endl;
//std::cout << "&" << std::endl;
//for(int il=0; il< *nr; ++il)
 //   std::cout << pr[il]*0.28438/tr[il] << "  " << zr[il] << std::endl;
//std::cout << "&" << std::endl;

//exit(23);

	float w, re;
	std::vector<float> profile( *nr);
	for(int ig = 0; ig< *ngdo; ++ig) {

		assert( top_dz[ig] >= 0.0);
		assert(bot_dz[ig] >= 0.0);
/*
		float ztop=0.0 ,zbot=900.0;
		for(int ich=0; ich< *nch; ++ich) {
			ztop = std::max(zstart[ich][ig], ztop);
			zbot = std::min(zend[ich][ig], zbot);
		}
*/
//std::cout << "ztop " << *ztop << " zbot " << *zbot << std::endl; 
		float* z1 = std::find_if(za, za+ *nl, std::bind2nd(std::less_equal<float>(),  ztop[ig] ) ) ;
		float* z2 = std::find_if(za, za+ *nl, std::bind2nd(std::less<float>(),  zbot[ig] ) ) -1 ;
		int istart = (int)std::distance(za, z1);
		int istop = (int)std::distance(za, z2);
//std::cout << "#z1 " << *z1  << " z2 " << *z2  << " zreg " << zreg[0] << " dz top " << top_dz[0] << std::endl;
//exit(23);
		if(igdo[ig] < 0)  {    // temperature merge

/*
for(int il = 0; il < *nr; ++il )
    std::cout << tr[il] << " " << pr[il] << std::endl;
std::cout << "&" << std::endl;
for(int jt=istart; jt<=istop ; ++jt)
	std::cout << ta[jt] << "  " << pa[jt] << std::endl;
std::cout << "&" << std::endl;
exit(23);
*/
			::liner_(za, ta, zr, &profile[0], nl, nr);

//std::cout << iscale_top[ig] << " " << top_dz[ig] << std::endl; exit(23);
			if(iscale_top[ig])  {
				float tt1, tt2;
				float zz = za[istart] - (0.5* top_dz[ig]);
				::liner_(za, ta, &zz, &tt1, nl, &c__1); 
				::liner_(zr, tr, &zz, &tt2, nr, &c__1);	
				int il=0;
				float delta = tt2-tt1;
//std::cerr << "**delta " << delta << " " << tt2 << " " << tt1 << " " << zz << " " << 
//ta[istart] << std::endl;
//exit(23);
				do {
					tr[il] -= delta;
					++il;
				} while(zr[il] >= za[istart]-top_dz[ig]);
			}
//std::cout << " scale bot " << iscale_bot[ig] << std::endl;
			if(iscale_bot[ig]) {
                                float tt1, tt2;
				float zz = za[istop] + (0.5* bot_dz[ig]);
				::liner_(za, ta, &zz, &tt1, nl, &c__1);
                                ::liner_(zr, tr, &zz, &tt2, nr, &c__1);
                                float delta = tt2-tt1;
                                for(int il=0; il< *nr; ++il)  {
					if(zr[il] > za[istop]+bot_dz[ig]) continue;
                                        tr[il] -= delta;
                                } 
			}

//for(int il = 0; il < *nr; ++il )
 //   std::cout << tr[il] << " " << pr[il] << std::endl;
//exit(23);

			std::replace_if(tr,tr+ *nr,std::bind2nd(std::less<float>(),110.), 110.) ;
  			std::replace_if(tr,tr+ *nr,std::bind2nd(std::greater<float>(),990.), 990.) ;

			for(int il = 0; il < *nr; ++il ) {
				if(zr[il] >= za[istart] || zr[il] <= za[istop]) continue;
				
				if(zr[il] >  za[istart]-top_dz[ig] ) {
					w = (zr[il]- za[istart])/ -top_dz[ig] ;
					tr[il] = w*profile[il] +  (1.0-w)*tr[il];
				} else if (zr[il] >= za[istop]+bot_dz[ig] ) {
					tr[il] = profile[il] ;
				} else {
					w = (zr[il]- (za[istop]+bot_dz[ig] ))/ -bot_dz[ig] ;
					tr[il] = (1.0- w)*profile[il] + w * tr[il];
				}
			}

//			float zreg_local = *zreg;
//			if( *zreg < za[istop]) zreg_local = za[istop];
			
			float* zhydro = std::find_if(zr, zr+ *nr, std::bind2nd(std::less<float>(), *zreg ) )-1;
			
//std::cout << "this is zhyrdro " << *zhydro << std::endl;  exit(23);
			int iptie = (int)std::distance(zr, zhydro)+1;
if( *zhydro >= za[istop] && *zhydro <= za[istart] )  {
// rebuild  pressures based on retrieved values 
			float wave=1000.,zt,  tt,  tearth;
			int mpout=200;
			std::vector<float> tgr(*nl, 0.0), gr(2),pout(mpout);
			int c__0 = 0;
			::kp_mpath__(nl, &c__1, &mpout, &c__0, &c__1, &c__1, &c__1, earth_rad,
				&c__1, nl, &wave, &c__0, &zr[iptie-1], 
				za, pa, ta, &tgr[0], &zt, &pr[iptie-1], 
			 	&tt, &gr[0], &tearth, &pout[0] );
}
//std::cout << "# p at iptie " << pr[iptie-1] << " and z " << zr[iptie-1] <<" this is iptie  " << iptie << std::endl;
//exit(23);
			::hydro_(nr, &c__1, nr, &iptie,  &re, xlat, zr, pr, tr);
//for(int il = 0; il < *nr; ++il )
//	std::cout << tr[il] << " " << zr[il] << std::endl;
// std::cout <<"&" << std::endl;
//std::cout << "# zreg " << *zreg << " " << zr[iptie-1] << "  " << pr[iptie-1] << std::endl;
//exit(23);
//   for(int i = 0; i< *nr; ++i )
//     tempFile << tr[i] << "   " << zr[i] << std::endl;
//   tempFile << "&" << std::endl;

		} else {       // VMR or extinction

			int* ik = std::find( igs,  igs+ *ngg, igdo[ig] );
			int iatm = (int)std::distance(igs, ik);
// take these two out as they are below
            ik = std::find( igas,  igas+ *ngas, igdo[ig] );
            int iret = (int)std::distance(igas, ik);

//std::cout << "#time for VMRS " << std::endl;
//for(int i = 0; i< *nr; ++i )
 //   std::cout << qmix_r[iatm][i] << "   " << zr[i] << std::endl;
//std::cout << "&" << std::endl;
//for(int i = 0; i< *nr; ++i ) {
// float dummy =  qmix_r[iatm][i] * 1.e-6 * 6.022e23*100.*pr[i]/(8.3144*tr[i]);
//    std::cout << dummy << "   " << zr[i] << std::endl;
//}
//std::cout << "&" << std::endl;

//std::cout << "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" << std::endl;

//for(int i = 0; i< *nl; ++i )
  // std::cout << qmix[iret][i] << "   " << za[i] << std::endl;
//exit(23);

			if(igdo[ig] >= 9000 || ik == igs+ *ngg ) continue;    // go to next product

//			ik = std::find( igas,  igas+ *ngas, igdo[ig] );			
//			int iret = (int)std::distance(igas, ik);

/*
for(int i = 0; i< *nr; ++i )
	std::cout << qmix_r[iatm][i] << "   " << zr[i] << std::endl;

std::cout << "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" << std::endl;

for(int i = 0; i< *nl; ++i )
	std::cout << qmix[iret][i] << "   " << za[i] << std::endl;
exit(23);

std::cout << "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&" << std::endl;
std::cout << za[istart] << "  " << za[istop] << std::endl;


exit(32);

*/

			::liner_(za, &qmix[iret][0], zr, &profile[0], nl, nr);
//std::cout << iscale_top[ig] << " top  " << top_dz[ig] << " bottom " << iscale_bot[ig] << "  "  << bot_dz[ig] << std::endl;
			if(iscale_top[ig])  {
				float Q1, Q2;
				float zz = za[istart] - (0.5* top_dz[ig]);
				::liner_(za, &qmix[iret][0], &zz, &Q1, nl, &c__1); 
				::liner_(zr, &qmix_r[iatm][0], &zz, &Q2, nr, &c__1);
				int il=0;
				float ratio = Q1/Q2;
				do {
					qmix_r[iatm][il] *= ratio;
					++il;
				} while(zr[il] >= za[istart]-top_dz[ig]);
			}

			if(iscale_bot[ig] ) {
				float Q1, Q2;
				float zz = za[istop] + (0.5* bot_dz[ig]);
				::liner_(za, &qmix[iret][0], &zz, &Q1, nl, &c__1); 
				::liner_(zr, &qmix_r[iatm][0], &zz, &Q2, nr, &c__1);
				float ratio = Q1/Q2;
                                for(int il=0; il< *nr; ++il)  {
					if(zr[il] > za[istop]+bot_dz[ig]) continue;
                                        qmix_r[iatm][il] *= ratio;
                                } 
			}
//exit(23);
			for(int il = 0; il < *nr; ++il ) {
				if(zr[il] >= za[istart] || zr[il] <= za[istop]) continue;
				
				if(zr[il] >  za[istart]-top_dz[ig] ) {
					w = (zr[il]- za[istart])/ -top_dz[ig] ;
					qmix_r[iatm][il] = w*profile[il] +  (1.0-w)*qmix_r[iatm][il];
				} else if (zr[il] >= za[istop]+bot_dz[ig] ) {
					qmix_r[iatm][il] = profile[il] ;
				} else {
					w = (zr[il]- (za[istop]+bot_dz[ig] ))/ -bot_dz[ig] ;
					qmix_r[iatm][il] = (1.0- w)*profile[il] + w * qmix_r[iatm][il];
				}
			}
//std::cout <<igdo[ig] << std::endl;
//exit(23);
//if(igdo[ ig ] == 200) {

//   for(int i = 0; i< *nr; ++i )
//	 std::cout << qmix_r[iatm][i] << "   " << zr[i] << std::endl;
//   co2File << "&" << std::endl;
 
//   for(int i = 0; i< *nr; ++i )
//	 tempFile << tr[i] << "   " << zr[i] << std::endl;
 //  tempFile << "&" << std::endl;
//exit(23);
//} else if(igdo [ig ] == 100)  {
 //  for(int i = 0; i< *nr; ++i )
//	 h2oFile << qmix_r[iatm][i] << "   " << zr[i] << std::endl;
 //  h2oFile << "&" << std::endl;
//std::cout << "done with water " << std::endl;
//exit(23);
//}
//for(int i = 0; i< *nr; ++i ) {
// float dummy =  qmix_r[iatm][i] * 1.e-6 * 6.022e23*100.*pr[i]/(8.3144*tr[i]);
 //   std::cout << dummy << "   " << zr[i] << std::endl;
//}

//std::cout << istart << "   " << istop << std::endl;  exit(34);
//for(int i = istart; i<= istop ; ++i )
//	std::cout << qmix[iret ][i] << "   " << za[i] << std::endl;
//std::cout << "&" << std::endl;

//std::cout << "# end of VMRS " << std::endl;
//exit(23);

		}


	}
std::cout << "# Finished merge into clim " << std::endl;

}