
# Code developed by Kabi Raj Khatiwada (Email: Kabiraj.Khatiwada@icimod.org)
# This R code checks the water balance and produce the visual graph of long term ..
# water balance which needs to be zero at the daily scale
# This R code can be applied to both glaciers and non-glaciers catchments

rm(list=ls())
library(data.table)
library(ggplot2)
library(tidyverse)
library(knitr)
library(scales)

TimeloopFolder  <- "C:\\jams\\data\\J2K_Dudh Koshi_MtEverest_region\\output\\current\\"

#read header names
variable_head <- unlist(strsplit(scan(paste(TimeloopFolder,"TimeLoop.dat",sep=""),"",nlines = 1,skip = 5,sep = "\n"),split = "\t"))
#read data
timeloop <- fread(paste(TimeloopFolder,"TimeLoop.dat",sep=""),skip = 10)
timeloop <- timeloop[,1:length(variable_head)]
#assign column names
colnames(timeloop) <- c("Date",variable_head[-1])
timeloop$Date <- as.Date(timeloop$Date, format= "%Y-%m-%d")
timeloop$year <- format(as.Date(timeloop$Date), format= "%Y")
#creating an empty column if there is not any variable related to the glacier melt process 
if(!"iceRunoff"  %in% colnames(timeloop)) { timeloop$iceRunoff =0;}
if(!"snowStorage_G"  %in% colnames(timeloop)) {timeloop$snowStorage_G =0;}
if(!"glacStorage"  %in% colnames(timeloop)){ timeloop$glacStorage =0;}
if(!"glacierRunoff"  %in% colnames(timeloop)){timeloop$glacierRunoff =0;}
if(!"snowMelt_G"  %in% colnames(timeloop)){timeloop$snowMelt_G =0;}
if(!"rainRunoff"  %in% colnames(timeloop)){timeloop$rainRunoff =0;}

#Calculation for the Daily WaterBalance for 
WaterBalance <- timeloop %>%
  arrange(Date) %>%
  mutate(	Year = year,
			Date = Date,
			Input = precip + iceRunoff,#
			Output = actET + catchmentRD1_w + catchmentRD2_w + catchmentRG1_w + catchmentRG2_w,
			diffIntStor = intercStorage - lag(intercStorage, default = first(intercStorage)),
			diffSWE = snowTotSWE - lag(snowTotSWE, default = first(snowTotSWE)), 
			diffSnowS_G = snowStorage_G - lag(snowStorage_G, default = first(snowStorage_G)),
			diffactMPS = actMPS - lag(actMPS, default = first(actMPS)),
			diffactLPS = actLPS - lag(actLPS, default = first(actLPS)),
			diffactRG1 = actRG1 - lag(actRG1, default = first(actRG1)),
			diffactRG2 = actRG2 - lag(actRG2, default = first(actRG2)),
			diffchannelStorage = channelStorage_w - lag(channelStorage_w, default = first(channelStorage_w)),
			diffactDPS = actDPS - lag(actDPS, default = first(actDPS)),
			glacestor = glacStorage, 
			Storage = diffIntStor + diffSWE + diffSnowS_G + diffactMPS+ diffactLPS + diffactRG1 + diffactRG2 + diffactDPS + glacestor + diffchannelStorage, 
			WatBal = Input - Output - Storage 
			)
WaterBalance <- WaterBalance[-1,]

#Visualization of the timeLoop water balance 
ggplot(data=WaterBalance, mapping=aes(y =  WatBal, x = Date))+
	geom_line(color="blue") +
	theme_bw() + ylab("mm") + 
	ggtitle("Water Balance of the Basin")+  
	theme(plot.title = element_text(hjust = 0.5)) 
	# scale_y_continuous(labels = comma)
	
# Export the water balance figure
ggsave("WaterBalance.png", width=900, height=600,units = "mm", dpi=600,scale = .6, limitsize = F)
	
	
	
	
	
	
