// Code for reading in the fetch and depth data created by pre_proc.cpp OR avg_modify.cpp
// and computing wave heights for each of the 4 formulations (WEMO,SPM,TMA,and TMA). These are calculated 
// on the wind forcing provided (ADCIRC .22 file) in addition to the physical processes of friction, breaking,
// and shoaling as computed from each formulation's period.
// The final results are 5 global elevation files (ADCIRC .63 files), one for each parametric formulation, and
// one for the ensemble average.
//
// By: Samuel Carter Boyd
// Date: 04/13/2020

#include <cmath>  
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <cstdlib>
#include <math.h> 
#include <vector>
#include <sstream>
#include <algorithm>
#include <ctime>
#include <stdio.h>
#include <string.h>

#define pi 3.14159265

using namespace std;

//**********************************************************************************
//Here the manning's n values which correspond to optimal model performance from
//previous tests are used. WEMO does not have any additional physical formulations

double mannings_n_SPM=0.01;
double mannings_n_TMA=0.01;
double mannings_n_CEM=0.02;

//**********************************************************************************

char content[10],temp[10];
int x,i,n,di_low,di_high,node_tot,Angres,dir,dir_tot,c,node;

double g=9.7918;	//gravity at meanlat (28.1091 deg)

int NDSETSE,NP,DTDP,NSPOOLGE,IRTYPE,TIME,IT;

bool check;

string s;

double Hw,Hwabove,Hwbelow,u,di_temp,di_input,y3,PRN, HWEMO,HSPM,HCEM,HTMA,WSPM,TSPM,TCEM,TTMA,Tw;

double HwaboveWEMO,HwbelowWEMO,HwaboveSPM,HwbelowSPM,HwaboveTMA,HwbelowTMA,HwaboveCEM,HwbelowCEM;
double TaboveSPM,TbelowSPM,TaboveTMA,TbelowTMA,TaboveCEM,TbelowCEM;

double WEMO_shoal,WEMO_fric,WEMO_cum,H_break;
double kSPM,LSPM,ksSPM,kfSPM,SPM_shoal,SPM_fric,SPM_cum;
double kCEM,LCEM,ksCEM,kfCEM,CEM_shoal,CEM_fric,CEM_cum;
double kTMA,LTMA,ksTMA,kfTMA,TMA_shoal,TMA_fric,TMA_cum;
double ENS_cum,H_break_SPM;

double sigma;
double tol;
double eps;
double knew,kold,k0;
double k;
double L;
double ks,kf,Tm,fm;


int line;

int Ftop,Ftop2,Favg,v=0;
double Dtop,Davg,dconst=0;

//*****************TMA formulation (EQ 19 from TMA paper)

double WaveTMA(double g, double f, double d, double u){

	double alpha;
	alpha = 0.076*pow( ((g*f)/(u*u)),-0.22);
	
	double fm;
	fm = 3.5*(g/u)*pow( (g*f/(u*u)),-0.33);
	
	double ET;
	ET = alpha*g*d/(16*pi*pi*0.9*0.9)*1/(fm*fm);

	double WTMA;
	WTMA = 4*pow(ET,0.5);

return WTMA;
}

double PeriodTMA(double g, double f, double d, double u){

	double fm = 3.5*(g/u)*pow(((g*f)/(u*u)),-0.33);

	TTMA = 1/fm;

return TTMA;
}

//***************CEM formulation (Eq II-2-36 in CEM)

double WaveCEM(double g, double f, double u){

	double Cd;
	Cd = 0.001*(0.009042*u + (-4.44*pow(10,-8)*f*f + 3.56*pow(10,-4)*f + 1.10949));	//Colvin
	//Cd = 0.001*(1.1+0.035*u);

	double Ustsq;
	Ustsq = Cd*u*u; 

	double WCEM;
	WCEM = (Ustsq/g)*0.0413*pow((g*f/Ustsq),0.5);

return WCEM;
}

double PeriodCEM(double g, double f, double u){

	double Cd;
	Cd = 0.001*(0.009042*u + (-4.44*pow(10,-8)*f*f + 3.56*pow(10,-4)*f + 1.10949));	//Colvin
	//Cd = 0.001*(1.1+0.035*u);

	double Ustsq;
	Ustsq = u*u*Cd; 

	double TCEM;
	TCEM = (pow(Ustsq,0.5)/g)*0.651*pow((g*f/Ustsq),0.3333333);

	//if(node==34283){
	//cout<<" Cd: "<<Cd<<" Ustsq: "<<Ustsq<<" TCEM: "<<TCEM<<endl;
	//}

return TCEM;
}


//**************SMB formulation in WEMO paper

double WaveWEMO(double g, double f, double d, double u){		//from WEMO

	Hw = 0.283*tanh(0.530*pow((g*d/(u*u)),0.75))*tanh((0.0125*pow((g*f/(u*u)),0.42))/
	           tanh(0.530*pow((g*d/(u*u)),0.75)))*(u*u)/g;

return Hw;
}

//*************SMB from SPM H and T formulations Eqs (3-39) and (3-40)

double PeriodSPM(double g, double f, double d, double u){

	double Ua = 0.713*pow(u,1.23);

	Tw =  7.54*tanh(0.833*pow((g*d/(Ua*Ua)),0.375))*tanh((0.0379*pow((g*f/(Ua*Ua)),(0.333333)))/
	           tanh(0.833*pow((g*d/(Ua*Ua)),0.375)))*(Ua/g);

return Tw;
}

double WaveSPM(double g, double f, double d, double u){			//from SPM

	double Ua = 0.713*pow(u,1.23);

	WSPM = 0.283*tanh(0.530*pow((g*d/(Ua*Ua)),0.75))*tanh((0.00565*pow((g*f/(Ua*Ua)),0.50))/
	             tanh(0.530*pow((g*d/(Ua*Ua)),0.75)))*(Ua*Ua)/g;

return WSPM;
}

//**************Linear Interpolation Function

double Interp(double x1, double x2, double x3, double y1, double y2){

	y3 = ((x2-x3)*y1 + (x3-x1)*y2) / (x2-x1);

return y3;
}

int main(){

	float clock_t,time_req;
	time_req=clock();			//initializing the clock for timing tasks


	ifstream infile1;
	infile1.open("Eff_fetch_2_deg.txt");		//reading in fetch data

	if(infile1.is_open()){

		cout<<endl;
		cout<<"Reading in fetch data"<<endl;

	}

	else{

		cout<<endl;
		cout<<"**************ERROR Fetch input file not read"<<endl;
		cout<<endl;

	}



//******These numbers need to be read in from EFF_fetch, make the format as such

	x=0;

	while(infile1 >> content){		//reading headers	

		x=x+1;					//content counter;

		if(x==1){

			node_tot=atoi(content);

		}

		if(x==2){

			Angres=atoi(content);
			dir_tot=360/Angres;
			break;

		}

	}

	cout<<"Node total: "<<node_tot<<endl;
	cout<<"Angular resolution: "<<Angres<<endl;
	cout<<"Number of angles: "<<dir_tot<<endl;

	double F[node_tot+1][dir_tot];
	double d[node_tot+1][dir_tot];
	double WVX[node_tot+1];
	double WVY[node_tot+1];
	double WDIR[node_tot+1];
	double WVEL[node_tot+1];
	double slope[node_tot][dir_tot];

	infile1.close();

	infile1.open("Eff_fetch_2_deg.txt");

	x=-3;
	n=0;

	while(infile1 >> content){

		x=x+1;

		if(x%(dir_tot+1)==0){		

			n=n+1;				//node counter increment
			i=-1;				//direction counter reset at new node 

		}

			i=i+1;				//direction counter increment
			F[n][i]=atof(content);		//fetch values

		//if(n<15){

		//	cout<<"F["<<n<<"]["<<i<<"]: "<<F[n][i]<<endl;

		//}
	
	}

	infile1.close();

	ifstream infile2;
	infile2.open("IDW_depths_2_deg.txt");		//reading in depth data

	if(infile2.is_open()){

		cout<<endl;
		cout<<"Reading in depth data"<<endl;

	}

	else{

		cout<<endl;
		cout<<"**************ERROR IDW depth input file not read"<<endl;
		cout<<endl;

	}

	n=0;					//reset counters
	x=-1;
	i=0;

	while(infile2 >> content){			

		x=x+1;					//content counter;

		if(x%(dir_tot+1)==0){		

			n=n+1;				//node counter increment
			i=-1;				//direction counter reset at new node 
					//dir index (i) = 0 corresponds to node number

		}
			i=i+1;				//direction counter increment
			d[n][i]=atof(content);		//depth values

		if(n==126765 || n==126769 || n==126770 || n==126771 || n==126772){	//for nan values

			d[n][i]=0;
			//cout<<n<<endl;

		}
	
	}


	infile2.close();

	ifstream infile3;				//reading wind field data
	infile3.open("PARAM_const_U_30.22");

	if(infile3.is_open()){

		cout<<endl;
		cout<<"Reading in wind data (fort.22)"<<endl;

	}

	else{

		cout<<endl;
		cout<<"**************ERROR wind input file not read"<<endl;
		cout<<endl;

	}
						
	NDSETSE=1;
	x=0;

	while(infile3 >> content){		//counting NDSETSE number of time steps	

		x=x+1;	
	
		if(x<=node_tot*4){

		}

		else{	
			
			x=1;
			NDSETSE=NDSETSE+1;	

		}

	}

	infile3.close();

	cout<<"NDSETSE: "<<NDSETSE<<endl;

//avg bottom slope

	ifstream infile4;
	infile4.open("Slope_upwind_2_deg.txt");		//reading in bottom slope data

	x=-1;
	n=0;

	while(infile4 >> content){

		x=x+1;

		if(x%(dir_tot+1)==0){		

			n=n+1;				//node counter increment
			i=-1;				//direction counter reset at new node 

		}

			i=i+1;				//direction counter increment
			slope[n][i]=atof(content);		//fetch values

		//if(n<15){

			//cout<<"slope["<<n<<"]["<<i<<"]: "<<slope[n][i]<<endl;

		//}
	
	}

	infile4.close();


	NP=node_tot;	//# of nodes
	DTDP=1;		//ADCIRC time step in seconds
	NSPOOLGE=3600;	//output is written to fort.63 every NSPOOLGE time steps
	IRTYPE=1;	//the record type (1 for elevation files, 2 for velocity files) 
	TIME=1;		//model time (in seconds) (TIME = STATIM*86400 + IT*DTDP)
	IT=1;		//model time step number since the beginning of the model run



	ofstream outfile14("WEMO_Hcum_const_wind.63");

	outfile14<<"Wave Height SMB TEST:3  U="<<u<<"(m/s)  Dir="<<di_temp<<"(deg N)"<<" from Mesh.grd"<<endl;
	outfile14<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
	outfile14<<NSPOOLGE<<"  "<<NSPOOLGE<<endl; //should be TIME<<"  "<<IT<<endl; if they are not equal

	ofstream outfile15("TMA_Hcum_const_wind.63");

	outfile15<<"Wave Height SMB TEST:3  U="<<u<<"(m/s)  Dir="<<di_temp<<"(deg N)"<<" from Mesh.grd"<<endl;
	outfile15<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
	outfile15<<NSPOOLGE<<"  "<<NSPOOLGE<<endl; //should be TIME<<"  "<<IT<<endl; if they are not equal

	ofstream outfile16("CEM_Hcum_const_wind.63");

	outfile16<<"Wave Height SMB TEST:3  U="<<u<<"(m/s)  Dir="<<di_temp<<"(deg N)"<<" from Mesh.grd"<<endl;
	outfile16<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
	outfile16<<NSPOOLGE<<"  "<<NSPOOLGE<<endl; //should be TIME<<"  "<<IT<<endl; if they are not equal

	ofstream outfile17("SPM_Hcum_const_wind.63");

	outfile17<<"Wave Height SMB TEST:3  U="<<u<<"(m/s)  Dir="<<di_temp<<"(deg N)"<<" from Mesh.grd"<<endl;
	outfile17<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
	outfile17<<NSPOOLGE<<"  "<<NSPOOLGE<<endl; //should be TIME<<"  "<<IT<<endl; if they are not equal

	ofstream outfile18("ENS_Hcum_const_wind.63");

	outfile18<<"Wave Height SMB TEST:3  U="<<u<<"(m/s)  Dir="<<di_temp<<"(deg N)"<<" from Mesh.grd"<<endl;
	outfile18<<NDSETSE<<"     "<<NP<<"     "<<DTDP*(NSPOOLGE)<<"     "<<NSPOOLGE<<"     "<<IRTYPE<<endl;
	outfile18<<NSPOOLGE<<"  "<<NSPOOLGE<<endl; //should be TIME<<"  "<<IT<<endl; if they are not equal
	

	cout<<"Time Step: "<<IT<<" at "<<NSPOOLGE<<" seconds"<<endl;

	infile3.open("PARAM_const_U_30.22");		//opening wind field data again

	//ofstream outfile3("fort.22_out");
	
	x=0;
	c=0;
	node=0;
	//int Q=0;

	//while(infile3 >> content){
	//if(infile3.is_open()){	

	x=0;
	c=0;
	node=0;
	line=1;
	IT=1;
	int Q=0;

	int y=0;

	HWEMO=-99999;
	HSPM= -99999;
	HCEM= -99999;
	HTMA= -99999;

	while(getline(infile3,s)){	

		//outfile3<<s<<endl;

		node=node+1;

		if(node<=126772){

			stringstream ss1;
			ss1<<s;			//read the line as string (s)
			ss1>>content;
			y=0;

			//if(node<15 && IT==1){

				//cout<<s<<endl;

			//}

			//cout<<s<<endl;

			while(!ss1.eof()){	//while values in the line (string ss)

				y=y+1;	

				ss1>>temp;	//read each value

				if(y==1){

					WVX[node]=atof(temp);

				}

				if(y==2){

					WVY[node]=atof(temp);

				}

				if(y==3){

					PRN=atof(temp);
					WVEL[node] = sqrt(WVX[node]*WVX[node] + WVY[node]*WVY[node]);


					if(WVX[node] <= 0 && WVY[node] < 0){		//Q1 [0,90)

						WDIR[node] = atan(WVX[node]/WVY[node]) * 180/pi;
						//WDIR[node] = 180+atan(WVX[node]/WVY[node]) * 180/pi;
					}

					if(WVX[node] < 0 && WVY[node] >= 0){		//Q4 [90,180)

						WDIR[node] = 90+atan(abs(WVY[node]/WVX[node]))*180/pi;
						//WDIR[node] = 270+atan(abs(WVY[node]/WVX[node]))*180/pi;
					}

					if(WVX[node] >= 0 && WVY[node] > 0){		//Q3 [180,270)

						WDIR[node] = 180+atan(WVX[node]/WVY[node]) * 180/pi;
						//WDIR[node] = atan(WVX[node]/WVY[node]) * 180/pi;
					}

					if(WVX[node] > 0 && WVY[node] <= 0){		//Q2 [270,360)

						WDIR[node] = 270+atan(abs(WVY[node]/WVX[node]))*180/pi;
						//WDIR[node] = 90+atan(abs(WVY[node]/WVX[node]))*180/pi;
					}

					u = WVEL[node];

					di_input=WDIR[node]/Angres;//di_input is interpolation direction index
					di_low=floor(di_input);	//rounds down to the direction index below
					di_high=ceil(di_input);	//rounds up to the direction index above

					//casting the double to int for the evenly divisible case

					dir = int(di_input);

		//if converged (evenly divisible case) no interpolation


if(node==1){
//cout<<"BEFORE, WDIR: "<<WDIR[node]<<"  di_input: "<<di_input<<"  di_low: "<<di_low<<"  di_high: "<<di_high<<"  dir: "<<dir<<" IT: "<<IT<<endl;
//cout<<"dir: "<<dir<<endl;
}


				//if( abs(((WDIR[node]+Angres)/Angres)-(2*IT)) < 0.01 || di_low==di_high){ 

				if(abs(di_input-di_low) < 0.01 || abs(di_input-di_high) < 0.01){

						if(node==1){

					cout<<"NO interpolation at time "<<IT; 

						}

						if(abs(di_input-di_low) < 0.01){	//round down

							WDIR[node]=floor(WDIR[node]);

						}

						if(abs(di_input-di_high) < 0.01){	//round up

							WDIR[node]=ceil(WDIR[node]);

						}

						di_input=WDIR[node]/Angres;

						//+1 is because 0 deg N corresponds to index 1

						dir = int(di_input)+1;



if(node==1){
cout<<"	   WVEL: "<<WVEL[node]<<"  WDIR: "<<WDIR[node]<<endl;
}

						HWEMO = WaveWEMO (g,F[node][dir],d[node][dir],u);
						HSPM  = WaveSPM  (g,F[node][dir],d[node][dir],u);
						HCEM  = WaveCEM  (g,F[node][dir],u);
						HTMA  = WaveTMA  (g,F[node][dir],d[node][dir],u);
						TSPM  = PeriodSPM(g,F[node][dir],d[node][dir],u);
						TCEM  = PeriodCEM(g,F[node][dir],u);
						TTMA  = PeriodTMA(g,F[node][dir],d[node][dir],u);

						if(F[node][dir]==0 || d[node][dir]==0 || HWEMO < 0.0001){

							HWEMO=-99999;
							HSPM =-99999;
							HTMA =-99999;
							HCEM =-99999;

						}
	
					}	//close if converged (evenly divisible no interpolation)

//*****case for an angle in between (not evenly divisible by Angres) using linear interpolation

					else{  

//case for setting maxindex(360)back to 0 deg

						if(WDIR[node]>360-Angres){ 

							di_high=1;

						}

						if(node==1){

		cout<<"     Interpolation between "<<di_low*Angres<<" and "<<di_high*Angres<<" degrees"<<endl;
//cout<<"F[node][di_high+1]: "<<F[node][di_high+1]<<" F[node][di_low+1]: "<<F[node][di_low+1]<<endl;
//cout<<"d[node][di_high+1]: "<<d[node][di_high+1]<<" d[node][di_low+1]: "<<d[node][di_low+1]<<endl;
//cout<<"WDIR: "<<WDIR[node]<<"  di_input: "<<di_input<<"  di_low: "<<di_low<<"  di_high: "<<di_high<<"  dir: "<<dir<<endl;

						}
//wave height and period values for interpolation (1 index above and 1 index below)

					HwaboveWEMO = WaveWEMO(g,F[node][di_high+1],d[node][di_high+1],u);
					HwbelowWEMO = WaveWEMO(g,F[node][di_low+1],d[node][di_low+1],u);
					HwaboveSPM  = WaveSPM(g,F[node][di_high+1],d[node][di_high+1],u);
					HwbelowSPM  = WaveSPM(g,F[node][di_low+1],d[node][di_low+1],u);
					HwaboveTMA  = WaveTMA(g,F[node][di_high+1],d[node][di_high+1],u);
					HwbelowTMA  = WaveTMA(g,F[node][di_low+1],d[node][di_low+1],u);
					HwaboveCEM  = WaveCEM(g,F[node][di_high+1],u);
					HwbelowCEM  = WaveCEM(g,F[node][di_low+1],u);

					TaboveSPM = PeriodSPM (g,F[node][di_high+1],d[node][di_high+1],u);
					TbelowSPM = PeriodSPM (g,F[node][di_low+1],d[node][di_low+1],u);
					TaboveTMA = PeriodTMA (g,F[node][di_high+1],d[node][di_high+1],u);
					TbelowTMA = PeriodTMA (g,F[node][di_low+1],d[node][di_low+1],u);
					TaboveCEM = PeriodCEM (g,F[node][di_high+1],u);
					TbelowCEM = PeriodCEM (g,F[node][di_low+1],u);

//resetting index to 36 to account for di_high

						if(di_temp>360-Angres){

							di_high=360/Angres;

						}

						HWEMO = Interp(di_low,di_high,di_input,HwbelowWEMO,HwaboveWEMO);
						HSPM  = Interp(di_low,di_high,di_input,HwbelowSPM,HwaboveSPM);
						HTMA  = Interp(di_low,di_high,di_input,HwbelowTMA,HwaboveTMA);
						HCEM  = Interp(di_low,di_high,di_input,HwbelowCEM,HwaboveCEM);
						TSPM  = Interp(di_low,di_high,di_input,TbelowSPM,TaboveSPM);
						TTMA  = Interp(di_low,di_high,di_input,TbelowTMA,TaboveTMA);
						TCEM  = Interp(di_low,di_high,di_input,TbelowSPM,TaboveCEM);

						if(F[node][dir]==0 || d[node][dir]==0 || HWEMO < 0.0001){

							HWEMO=-99999;
							HSPM =-99999;
							HTMA =-99999;
							HCEM =-99999;

						}

						//outfile3<<node<<"  "<<HWEMO<<endl;

					}	//close interpolation else statement


					check=false;

					if(F[node][dir]==0 || d[node][dir]==0 || HWEMO < 0.0001){

						HWEMO=-99999;
						HSPM= -99999;
						HCEM= -99999;
						HTMA= -99999;
						WEMO_cum=-99999;
						TMA_cum =-99999;
						SPM_cum =-99999;
						CEM_cum =-99999;
						ENS_cum =-99999;
						TSPM=0;
						TCEM=0;
						TTMA=0;
						check=true;

					}


//*******Dispersin Equation for k and L (WEMO not included bc it does not have a period formulation)

					if(check==false){//wet node with non-zero fetch and depth values


						//SPM
						Tm=TSPM;
						fm=1/Tm;
						sigma=2*pi*fm;		
						k0=sigma*sigma/g;
						kold=k0;
						eps=10;

						int b=0;

						while(eps>0.0001){

					 		knew=sigma*sigma/(g*tanh(kold*d[node][dir]));
							eps=abs(knew-kold);
							kold=knew;
							b=b+1;

							if(b>1000000){	//statement for no convergence
						//cout<<"SPM NOT CONVERGED at "<<node<<" eps: "<<eps<<endl;
								break;
							}
						}
						kSPM=knew;
						LSPM=2*pi/kSPM;

						//CEM
						Tm=TCEM;
						fm=1/Tm;
						sigma=2*pi*fm;		
						k0=sigma*sigma/g;
						kold=k0;
						eps=10;

						b=0;

						while(eps>0.0001){

					 		knew=sigma*sigma/(g*tanh(kold*d[node][dir]));
							eps=abs(knew-kold);
							kold=knew;
							b=b+1;

							if(b>1000000){
						//cout<<"CEM NOT CONVERGED at "<<node<<" eps: "<<eps<<endl;
								break;
							}
						}
						kCEM=knew;
						LCEM=2*pi/kCEM;

						//TMA
						Tm=TTMA;
						fm=1/Tm;
						sigma=2*pi*fm;		
						k0=sigma*sigma/g;
						kold=k0;
						eps=10;

						b=0;

						while(eps>0.0001){

					 		knew=sigma*sigma/(g*tanh(kold*d[node][dir]));
							eps=abs(knew-kold);
							kold=knew;
							b=b+1;

							if(b>1000000){
						//cout<<"TMA NOT CONVERGED at "<<node<<" eps: "<<eps<<endl;
								break;
							}
						}
						kTMA=knew;
						LTMA=2*pi/kTMA;


double Ch_SPM=(1.486/mannings_n_SPM)*pow(d[node][dir],1/6);
double Ch_TMA=(1.486/mannings_n_TMA)*pow(d[node][dir],1/6);
double Ch_CEM=(1.486/mannings_n_CEM)*pow(d[node][dir],1/6);
double ffactor_SPM=g/(Ch_SPM*Ch_SPM);
double ffactor_TMA=g/(Ch_TMA*Ch_TMA);
double ffactor_CEM=g/(Ch_CEM*Ch_CEM);


//shoaling coefficient ks (varies from 1 in deep water, dips to 0.9 in intermediate and grows larger than 1 as it approaches shallow water. Multiply H by this)

ksCEM = pow(((2*pow(cosh(kCEM*d[node][dir]),2)) / (sinh(2*kCEM*d[node][dir]) + 2*kCEM*d[node][dir]) ),0.5);
ksSPM = pow(((2*pow(cosh(kSPM*d[node][dir]),2)) / (sinh(2*kSPM*d[node][dir]) + 2*kSPM*d[node][dir]) ),0.5);
ksTMA = pow(((2*pow(cosh(kTMA*d[node][dir]),2)) / (sinh(2*kTMA*d[node][dir]) + 2*kTMA*d[node][dir]) ),0.5);

//friction factor kf (varies from 1 in deep water to larger than 1 in shallow. Divide H by this)

kfCEM = 1 + (64*pi*pi*pi)/(3*g*g) * (ffactor_CEM*HCEM*F[node][dir]/pow(TCEM,4)) * 
	(ksCEM*ksCEM/pow(sinh(2*pi*d[node][dir]/LCEM),3));
kfSPM = 1 + (64*pi*pi*pi)/(3*g*g) * (ffactor_SPM*HSPM*F[node][dir]/pow(TSPM,4)) * 
	(ksSPM*ksSPM/pow(sinh(2*pi*d[node][dir]/LSPM),3));
kfTMA = 1 + (64*pi*pi*pi)/(3*g*g) * (ffactor_TMA*HTMA*F[node][dir]/pow(TTMA,4)) * 
	(ksTMA*ksTMA/pow(sinh(2*pi*d[node][dir]/LTMA),3));


//SPM / WEMO

//WEMO_shoal = HWEMO*ksSPM;
//WEMO_fric  = HWEMO/kfSPM;
//WEMO_cum   = HWEMO - (HWEMO-WEMO_shoal) - (HWEMO-WEMO_fric);
WEMO_cum   = HWEMO;		//wind only

SPM_shoal  = HSPM*ksSPM;
SPM_fric   = HSPM/kfSPM;
SPM_cum    = HSPM - (HSPM-SPM_shoal) - (HSPM-SPM_fric);

//if(ksSPM>1 && IT==1){
//cout<<"SPM positive shoal effect at node: "<<node<<" ksSPM: "<<ksSPM<<endl;
//}

//checking for breaking criteria. If d/L < 0.05 then set cumulative Hs equal to Goda breaking formulation

double A=0.17;
double H_break=0;

if(d[node][dir]/LSPM < 0.05){

	H_break_SPM  = LSPM*A*(1-exp (-1.5* (pi*d[node][dir]/LSPM)*(1+15*pow(slope[node][dir],4/3)) ) );
	//if(IT==1){
	//cout<<"SPM BREAK!!! at node:"<<node<<" SPM_cum: "<<SPM_cum<<" SPM_break: "<<H_break_SPM<<endl;
	//}
	WEMO_cum = H_break_SPM;
	SPM_cum  = H_break_SPM;
}
else{
	H_break_SPM = -99999;
}

//CEM

CEM_shoal = HCEM*ksCEM;
CEM_fric  = HCEM/kfCEM;
CEM_cum   = HCEM - (HCEM-CEM_shoal) - (HCEM-CEM_fric);

if(d[node][dir]/LCEM < 0.05){

	H_break = LCEM*A*(1-exp (-1.5* (pi*d[node][dir]/LCEM)*(1+15*pow(slope[node][dir],4/3)) ) );
	CEM_cum = H_break;
	//cout<<"BREAK!!! node:"<<node<<endl;
}

//TMA

TMA_shoal = HTMA*ksTMA;
TMA_fric  = HTMA/kfTMA;
TMA_cum   = HTMA - (HTMA-TMA_shoal) - (HTMA-TMA_fric);

if(d[node][dir]/LTMA < 0.05){

	H_break = LTMA*A*(1-exp (-1.5* (pi*d[node][dir]/LTMA)*(1+15*pow(slope[node][dir],4/3)) ) );
	TMA_cum = H_break;
	//cout<<"BREAK!!! node:"<<node<<endl;
}


//statement to fix infinity values on some ks denominator values, in this case reset cumulative wave heights
//back to the wind wave height only value

if(isinf(sinh(2*kCEM*d[node][dir]) + 2*kCEM*d[node][dir])  ==1){	//denominator in ks, causing infinity
CEM_cum=HCEM;								//reset to wind wave height
}
if(isinf(sinh(2*kSPM*d[node][dir]) + 2*kSPM*d[node][dir])  ==1){	//denominator in ks, causing infinity
SPM_cum=HSPM;
WEMO_cum=HWEMO;								//reset to wind wave height
}
if(isinf(sinh(2*kTMA*d[node][dir]) + 2*kTMA*d[node][dir])  ==1){	//denominator in ks, causing infinity
TMA_cum=HTMA;								//reset to wind wave height
}


ENS_cum = (WEMO_cum + SPM_cum + TMA_cum + CEM_cum)/4;		//ensemble avg.

						}	//close if check=false


						else{	//fetch and depth are zero. Hshoal=-99999

						ks=1;	
						kf=1;					

						}


						outfile14<<node<<"  "<<WEMO_cum   <<endl;
						outfile15<<node<<"  "<<TMA_cum   <<endl;
						outfile16<<node<<"  "<<CEM_cum   <<endl;
						outfile17<<node<<"  "<<SPM_cum   <<endl;
						outfile18<<node<<"  "<<ENS_cum   <<endl;


if(node==35884){
//cout<<fixed<<setprecision(3)<<endl;
//cout<<"node: "<<node<<" HWEMO: "<<WEMO_cum<<" HSPM: "<<SPM_cum<<" HTMA: "<<TMA_cum<<" HCEM: "<<CEM_cum<< " ENS_cum: "<<ENS_cum<<endl;
//cout<<"dir: "<<dir<<" F: "<<F[node][dir]<<" d: "<<d[node][dir]<<" HSPM: "<<HSPM<<" SPM_cum: "<<SPM_cum<<endl;//" SPM_shoal: "<<SPM_shoal<<" SPM_fric: "<<SPM_fric<<endl;
//cout<<"node: "<<node<<" TWEMO: "<<TSPM<<" TSPM: "<<TSPM<<" TTMA: "<<TTMA<<" TCEM: "<<TCEM<<endl;
}

//statement to write wave height solutions of 0 to -99999 for ADCIRC and SMS

						HWEMO=0;
						HCEM =0;
						HTMA =0;
						HSPM =0;
						TSPM =0;
						TCEM =0;
						TTMA =0;
						WEMO_cum=-99999;
						TMA_cum =-99999;
						SPM_cum =-99999;
						CEM_cum =-99999;
						ENS_cum =-99999;



				}	//close if y==3

			}	//close while!eof

			//cout<<node<<"  "<<WVX[node]<<"  "<<WVY[node]<<"  "<<WDIR[node]<<"  "<<IT<<endl;

		}	//close if node <=126772


		if(node==126772){	//new timestep, reset node counter

			if(IT==NDSETSE){	//to break out of loop before writing next time step

				break;

			}
			
			IT=IT+1;			//model timestep number
			TIME=NSPOOLGE*IT; 		//model time (in seconds)
			node=0;


			outfile14<<TIME<<"  "<<TIME<<endl;
			outfile15<<TIME<<"  "<<TIME<<endl;
			outfile16<<TIME<<"  "<<TIME<<endl;
			outfile17<<TIME<<"  "<<TIME<<endl;
			outfile18<<TIME<<"  "<<TIME<<endl;

			cout<<"Time Step: "<<IT<<" at "<<TIME<<" seconds"<<endl;
			//cout<<"Number of SPM breaking nodes: "<<v<<endl;
			v=0;

		}


	}	//close while getline infile3

	//outfile3.close();

	outfile14.close();
	outfile15.close();
	outfile16.close();
	outfile17.close();
	outfile18.close();

	
	time_req=clock()- time_req;

	cout<<"***Time taken: "<<(float)time_req/CLOCKS_PER_SEC<<" seconds***"<<endl;

return 0;
}