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

#include <algorithm>
#include <iostream>

using namespace GATS_Utilities;

extern "C" void sofie_mega_model (
	int* iblk,	int* nch,
	int* idg,  	int* ng,	
	int* isl,	int* iel,	int* nl, 	int* ml,
	int **idg_b,	int **irad_b,	int **icor_b,
	int* itab,	int* mgs,	
	int* ntb,	int *ntt,	int *nst,
	float*  pout,	float ***rlt,	float ***rut,
	float* re,	float* tearth,	float* albedo,
	float **qmix,	float ***qg,	
	float ***tbl,	float *tout,	float **sout,
	double ***emr,  double ***rmr,  double ***esr,
	double ***rsr,	float ***bsr,
	double  **tau, float *zt, float* tapr)
//, float* zt, float* pt, float* tt,  float* xlat, float* za)
{
	static int ic1 = -1;		//beginning cell for integration
	static int ic2 = 1;		//ending cell for integration

//	static int old_nl = -1;
	static FORTRANarray<float> qmixB, qgB, rltB, rutB;

//std::cout << *nch << "  " << idg[0] << "  " << idg[1] << "  " << *ng << std::endl; 

// comment below

//if(iblk[0] == 13 ) {

/*
if(iblk[0] == 3 ) {
  std::cout << "# " << *ng << "   " << idg[0] << "  " << idg[1]  << "   " << iblk[0] <<  std::endl;
  for(int i= *isl-1; i< *iel ; ++i)
    std::cout << qmix[0][i] << "  " <<  zt[i] << std::endl;
//  std::cout << "&" << std::endl;
}

*/

// comment

	if(qmixB.size() == 0) {
//std::cout << " building arras " << std::endl;
		qmixB = FORTRANarray<float>(*ml, *mgs);
		qgB = FORTRANarray<float>(2, *ml, *mgs);
	    rltB =FORTRANarray<float>(*ml, *mgs);
	    rutB = FORTRANarray<float>(*ml, *mgs);
	}

	for(int ich = 0; ich< *nch; ++ich) {   // loop over the bands
		int idch = iblk[ich]-1;
//std::cout << *isl << " " << *iel << std::endl;
		int ngs = 0;
		std::vector<int>  iradB, icorB;	
		for(int ig = 0; ig< *mgs; ++ig) {
			if(idg_b[idch][ig] == 0 )  break;
//std::cout <<" idch " <<  idch << std::endl;
			++ngs;
			int* j = std::find(idg, idg+ *ng, idg_b[idch][ig] );
			int ik = (int)std::distance(idg, j);
			if( ik < *ng) {
//std::cout << "# pre " << ik << "  " << ig << "  "  << *ng << std::endl;
				iradB.push_back( irad_b[idch][ig] );
				icorB.push_back( icor_b[idch][ig] );
				std::copy( &qmix[ik][0], &qmix[ik][*nl], qmixB.getArray(1,ig+1) );
				std::copy( &rlt[ich][ik][0], &rlt[ich][ik][*nl], rltB.getArray(1,ig+1) );
                std::copy( &rut[ich][ik][0], &rut[ich][ik][*nl], rutB.getArray(1,ig+1) );
				std::copy( &qg[ik][0][0], &qg[ik][0][0]+ 2* *nl, qgB.getArray(1,1,ig+1) ); 
// Channel 2 DV
if(*nch == 2 && iblk[ich] == 3  && idg_b[idch][ig] == 9000  ) {
//std::cout << "mod extinction " << za[*isl-1] << " " << qmixB(*isl,ig+1) ;
   for(int i=1; i<= *nl; ++i) qmixB.setVal(qmixB(i,ig+1)*2.0,i,ig+1);   // make extinction 2 * ( band 4) 
//   std::cout << "   "   <<  qmixB(*isl,ig+1) << std::endl;
}

			} else {
				iradB.push_back( 0 );
				icorB.push_back( 0 );
			}
		}
		if( ngs == 0) continue;

//   comment below

/*
if(iblk[ich] == 16) {
  for(int i=1; i<= ngs; i++) {
  std::cout << "# " << idg_b[idch][i-1] << std::endl; 
      for(int j= 0; j< *nl; ++j ) 
          std::cout <<  qmixB(j+1,i) << "  " << zt[j]   << std::endl;
     std::cout << "&" << std::endl;
  }


//std::cout << iradB[0] << "  " << icorB[0] << std::endl;
exit(23);
}
*/

// comment above

//std::cout << "#before mega " << *isl << "  " << *iel << std::endl; 

//for(int i = 1; i<= *nl ; i++) {

 //std::cout <<    qmixB(i, 1) << "   " << qmixB(i,2) << std::endl;
//}

		::mp_mega__( ntb,  nst,  ntt,  mgs,    ml,    nl,   
			isl, iel,   &ic1,  &ic2,   itab,  &ngs,  &idg_b[idch][0],
			qmixB.getArray(), qgB.getArray(), 
			rltB.getArray(), rutB.getArray(),
			albedo, tearth,   &iradB[0],  &icorB[0],   &emr[0][0][0],
			&rmr[0][0][0],  &esr[0][0][0],  &rsr[0][0][0],
			&bsr[0][0][0],  &tbl[ich][0][0],   pout,
			&sout[0][0], tout  );

//if( *iexo ) {
//		for(int ir= *isl-1; ir <= *iel-1; ++ir) {
//			tau[ich][ir] =  emr[0][0][ir] ; // emr[0][0][ir] ;  // 
  //      }
//} else {
        for(int ir= *isl-1; ir <= *iel-1; ++ir) {
            tau[ich][ir] = (double)1.0 - emr[0][0][ir] ;

//std::cout << emr[0][1][ir] << "  " << zt[ir] << std::endl;
//std::cout <<  tau[ich][ir] << "  " << tapr[ir] << "  " << zt[ir] << std::endl;
        }
//}
//if(*iel - *isl > 10) exit(23);

//if(*iel - *isl > 10) exit(23);
	
/*
for(int ir= *isl-1; ir <= *iel-1; ++ir) 
    std::cout <<   ConvertToStringPrec(tau[0][ir] )   <<  std::endl;
    exit(23);
*/
   

	} // loop over bands

}