// LEB.cpp : Defines the entry point for the console application.
//

#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>

//#define pkdur 0.25
#define npulse 3
#define nyrs 200

int main(int argc, char* argv[])
{
	char dumhd[256];
	float soildep, por, drainpor, ks, fwr, maxlstor, soilstor, keff, drainstor;
	float prcp, dur, tp, ip, tmax, tmin, rad, wspd, wdir, tdew;
	float pkdur, p_int[npulse], p_dur[npulse], q[npulse],q_dur[npulse];
	float stor, lstor;
	float pvol,q_rate,re_dur,re_vol,qmax,qrmax,q_tot,q_dur_tot,re_tot;
	float fcap,evap,kcond,conductance;
	float maxq_tot[nyrs],temp;
	int da, mo, yr, ii, iii;
	
	ifstream par_file("LEBParam.dat");

	par_file >> soildep;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << soildep << "  " << dumhd << endl;
	
	par_file >> por;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << por << "  " << dumhd << endl;

	par_file >> drainpor;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << drainpor << "  " << dumhd << endl;

	par_file >> fwr;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << fwr << "  " << dumhd << endl;
	
	par_file >> ks;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << ks << "  " << dumhd << endl;

	par_file >> maxlstor;
	par_file.getline (dumhd, sizeof(dumhd));
	cout << maxlstor << "  " << dumhd << endl;

	par_file.close();

	
	keff = ks * (1-fwr);
	soilstor = soildep * por;
	drainstor = soildep * drainpor;
	fcap = soilstor - drainstor;
	
	for (ii=0; ii<=nyrs-1; ii++) {
		maxq_tot[ii]=0;
		}


	// begin simulation

	ifstream cli_file(argv[1]);
	ofstream out_file("lebout.txt");

	for (ii=0; ii<=13; ii++) {
		cli_file.getline(dumhd,sizeof(dumhd));
		}
	cout << dumhd << endl;
	cli_file.getline(dumhd,sizeof(dumhd));
	cout << dumhd << endl;
	
	while (! cli_file.eof()){
		cli_file >> da >> mo >> yr >> prcp >> dur >> tp >> ip >> tmax
			>> tmin >> rad >> wspd >> wdir >> tdew;
		
		if ((mo<6) || (mo>9)) {
			stor = soilstor - drainstor;
			}
		else {                      // summer months 

			// Handle Additions to storage and runoff generation
			
			// 1 Build hyetograph
		
			if (prcp > 0){
				pkdur = 0.25;
				p_int[1] = ip * prcp/dur;
				if (dur > pkdur) {p_int[0] = (prcp - p_int[1] * pkdur)/(dur - pkdur);}
				else{ p_int[0]=0;}
				p_int[2] = p_int[0];
				p_dur[0] = (dur - pkdur) * tp;
				p_dur[1] = pkdur;
				p_dur[2] = (dur - pkdur) - p_dur[0];

		
			// 2 Partition into runoff and infiltration by pulse
			
			qrmax = 0;
			qmax = 0;
			q_tot = 0;
			q_dur_tot = 0;
			re_tot = 0;

			re_vol = 0;
			for (ii=0; ii<=npulse-1; ii++) {
					q[ii] = 0;
					q_dur[ii]=0;
			}
						
			for (ii=0; ii<=npulse-1; ii++) {
				pvol = p_dur[ii] * p_int[ii];  // precip in pulse
				
				if  (stor + pvol < soilstor) {   // soil above wr layer unsaturated 
					stor = stor + pvol;
					q[ii] = 0;
					q_dur[ii]=0;
					q_rate = 0;
					re_vol = 0;
				}
				else if (p_int[ii] < keff)    {  //soil above wr sat but all infiltrates
					stor = soilstor;
					q[ii] = 0;
					q_dur[ii] = 0;
					q_rate = 0;
					re_vol = 0;
				}
				else  {                             //rainfall excess
					
					re_dur = p_dur[ii] - (soilstor-stor)/p_int[ii];
					re_vol = re_dur * (p_int[ii] - keff);

					if (lstor + re_vol < maxlstor) {   //stored behind logs
						stor = soilstor;
						lstor = lstor + re_vol;
						q[ii] = 0;
						q_dur[ii] = 0;
						q_rate = 0;
					}                                  //end lstor if
					else {                             //runoff
						stor = soilstor;
						q_dur[ii] = re_dur - (maxlstor - lstor) / (p_int[ii] - keff);
						q[ii] = q_dur[ii] * (p_int[ii] - keff);
						q_rate = q[ii] / q_dur[ii];
						lstor = maxlstor;
					}                                    // end lstor else
				}                                 // end partitioning if
				//now summarize pulse information
				
				q_tot = q_tot + q[ii];
				q_dur_tot = q_dur_tot + q_dur[ii];
				re_tot = re_tot + re_vol;
				
				if (q[ii] > qmax) qmax = q[ii];
				if (q_rate > qrmax) qrmax = q_rate;

			
			} // end partitioning "for"
			// summarize daily information

			if (q[1] > maxq_tot[yr]) maxq_tot[yr] = q[1];
				
			
			} // end precip if
			else {                 //no precip
				qmax = 0;
				qrmax = 0;
				for (ii=0; ii<=npulse-1; ii++) {
					p_int[ii] = 0;
					p_dur[ii]=0;
					q[ii]=0;
					q_dur[ii]=0;
				}
				q_tot=0;
				q_dur_tot=0;
				re_vol=0;
				re_tot=0;
			}	    //end no-precip else

			//Drain and evaporate

			lstor = 0;   //logs drain in a day?
			evap = 5;
			kcond = .5;
			conductance = exp(-kcond*(soilstor-stor));

			stor = stor - conductance * evap;
			if (stor < 0) stor = 0;
			if (stor > fcap) stor = fcap;



			//Output Results

			cout << da << "  " << mo << "  " << yr <<" " << prcp << "  " << p_int[1] << "  "
				<< qmax << "  " << qrmax << endl;


			out_file << da << "  " << mo << "  " << yr <<" " << prcp << " " << dur <<"  " << p_int[0] <<" " 
				<< p_int[1] << "  " << p_int[2] << "  " << p_dur[0] <<" " << p_dur[1] << "  " 
				<< p_dur[2] << " " <<q_tot << "  " << qrmax << "  " << re_tot << "  " << stor << endl;


			
			}   //end summer month if
		

	
		}       //end while for file input


		ofstream ann_file("lebannual.txt");
		 
		//sort annual max
		for (ii=1; ii<=yr-1; ii++){
			for (iii=1; iii<=yr-ii; iii++){
				if (maxq_tot[iii] < maxq_tot[iii+1]){
					temp = maxq_tot[iii+1];
					maxq_tot[iii+1] = maxq_tot[iii];
					maxq_tot[iii]=temp;
				} //end switch
			} // end embedded for
		}// end sort


		for (ii=1; ii<=yr; ii++) ann_file << maxq_tot[ii] << endl;
		
		ann_file.close();

		
		
		cli_file.close();
		out_file.close();



	return 0;
}