# TODO: Add comment
# 
# Author: jaredwo
###############################################################################
PARAMS_SNOW_SEELEY <- list()
PARAMS_SNOW_SEELEY$SFTMP<-1.576050
PARAMS_SNOW_SEELEY$SMTMP<-3.149728
PARAMS_SNOW_SEELEY$TIMP<-0.126985
PARAMS_SNOW_SEELEY$SMFMN<-0.000000*0.1 #convert from mm to cm
PARAMS_SNOW_SEELEY$SMFMX<-3.608737*0.1 #convert from mm to cm

calib.swat_snow <-function(x,stn.data)
{
	params = to.params(x)
	params$SFTMP<-x[1]
	params$SMTMP<-x[2]
	params$TIMP<-x[3]
	params$SMFMN<-x[4]*0.1 #convert from mm to cm
	params$SMFMX<-x[5]*0.1 #convert from mm to cm
	
	swe <- swat_snow(stn.data,params)
	return(sum((swe-stn.data$SWE)^2))
}

calib.swat_snow.3 <-function(x,stn.data)
{
	params <- list()
	params$SFTMP<-1
	params$SMTMP<-0
	params$TIMP<-x[1]
	params$SMFMN<-x[2]*0.1 #convert from mm to cm
	params$SMFMX<-x[3]*0.1 #convert from mm to cm
	
	swe <- swat_snow(stn.data,params)
	return(sum((swe-stn.data$SWE)^2))
}

to.params <- function(x)
{
	params = list()
	params$SFTMP<-x[1]
	params$SMTMP<-x[2]
	params$TIMP<-x[3]
	params$SMFMN<-x[4]*0.1 #convert from mm to cm
	params$SMFMX<-x[5]*0.1 #convert from mm to cm
	return(params)
}

swat_snow <- function(metdata,params=PARAMS_SNOW_SEELEY)
{
	swe = vector(mode="numeric",nrow(metdata))
	snw = vector(mode="numeric",nrow(metdata))
	yday = metdata$yday
	
	sftmp = params$SFTMP
	smtmp = params$SMTMP
	timp = params$TIMP
	smfmn = params$SMFMN
	smfmx = params$SMFMX
	
	tsnow_prev = 0.0
	swe_prev = 0.0
	
	for (x in seq(length(swe)))
	{
		snw_fall = calc_snw_fall(metdata$PRCP[x], metdata$TAVG[x],sftmp)
		sm_rate = calc_snw_melt_rate(smfmn, smfmx,yday[x])
		tsnow = calc_tsnow(tsnow_prev,timp, metdata$TAVG[x])
		snw_melt = calc_snw_melt(sm_rate, tsnow, metdata$TAVG[x],metdata$TMAX[x],smtmp)
		delta_snw = snw_fall - snw_melt

		swe_cur = swe_prev + delta_snw
		
		if (!is.finite(swe_cur))
			print(params)
	
		if ((swe_cur) < 0) 
		{
			swe_cur = 0.0
		}
		swe[x] = swe_cur
		snw[x] = snw_fall

		tsnow_prev = tsnow
		swe_prev = swe_cur
	}
	return(swe)	
}

calc_snw_fall<-function(prcp,tavg,SFTMP)
{
	if (tavg < SFTMP)
	{
		return(prcp)
	}				
	else
	{
		return(0.0)
	}	
}

calc_snw_melt<-function(sm_rate,tsnow,tavg,tmax,SMTMP)
{
	if (tavg <= SMTMP)
	{
		return(0.0)	
	}
	else
	{
		return(sm_rate*(((tsnow+tavg)/2.0) - SMTMP))
	}					
}


calc_snw_melt_rate<-function(SMFMN,SMFMX,yday)
{
	return (((SMFMX+SMFMN)/2.0)+(sin((yday-81)/58.09)*((SMFMX-SMFMN)/2.0)))
	#365/2pi = 58.09
}
		
calc_tsnow<-function(tsnow_prev,TIMP,tavg)
{
	return ((tsnow_prev*(1.0-TIMP)) + (tavg*TIMP))	
}