rm(list=ls())
library(lubridate)
library(dplyr) 
library(data.table) 
library(stringr)
library(foreach)
library(doSNOW)
library(powerjoin)
library(ggplot2)
library(sf)

options(scipen = 999)
nearZero <- 1E-28

sNOW = str_replace_all(Sys.time(), "[[:punct:]]", "")

#----PCLAKE CODE---
days_of_summer <- expand.grid(seq(91, 274, 1), seq(0, 30))
summer_vec <- days_of_summer$Var1 + (365 * days_of_summer$Var2)

dirHome <- file.path("C:", "Users", "SvenT", "Documents", "PCModel","PCModel-master", "Licence_agreement", "I_accept")
#dirHome <- file.path("C:", "Users", "LilithK", "Documents", "PCModel", "PCModel-master", "Licence_agreement", "I_accept")
#dirHome <- "C:/Users/SvenT/Documents/PCModel/PCModel-master/Licence_agreement/I_accept/"	# location of the PCModel1350 folder
dirShell <- file.path(dirHome, "PCModel1350", "PCModel", "3.00", "Models", "PCLake+", "6.13.16", "PCShell")
dirCpp_root <- file.path(dirHome, "PCModel1350", "PCModel", "3.00", "Frameworks", "Osiris", "3.01", "PCLake_plus")
nameWorkCase <- "PCLake_plus_DavidsonRun2"
fileDATM <- file.path(dirHome, "PCModel1350", "PCModel", "3.00", "Models", "PCLake+", "6.13.16", "PL613162PLUS_Davidson.xls")

dirSave = file.path("D:", "PCAEM","PCModel_TEO","PCModel-master","Licence_agreement","I_accept","PCModel1350", "PCModel", "3.00", "Models", "PCLake+", "6.13.16", "PCShell")

#library("LAGOSNE")

## load all the functions
source(file.path(dirShell, "scripts", "R_system", "functions.R"))  	#load base functions by Luuk van Gerven (2012-2016)
source(file.path(dirShell, "scripts", "R_system", "functions_PCLake.R")) 

## NB: you cannot use different combinations of forcings. e.g. 
##      - in model run 1: mQInEpi & mQInHyp
##      - in model run 2: mPLoad
##     If you want to do this you'll have to do for both: mQInEpi, mQInHyp AND mPLoad. 
##     This is due to the fact that the include file for the parameters (...rp.cpp file) has to be 
##     adjusted before compilation - the file will be built into the code.  
##     This adjustment entails changing parameters that are given as forcings should become "dummy". 
##     BTW: this has never been possible! 
##     How is this possible in the Excel version?

## Order of actions
##   1. Making folder structure for running the model
##   2. Load file 
##   < Make adjustments to the model > 
##   3. Make cpp files
##   4. Compile model
##   5. Initialize model
##   6. Run model


## 1. Making folder structure 
PCModelWorkCaseSetup(dirSHELL = dirShell, 
                     dirCPP_ROOT = dirCpp_root,
                     nameWORKCASE = nameWorkCase)

## 2. Load file
lDATM_SETTINGS <- PCModelReadDATMFile_PCLakePlus(fileXLS = fileDATM,
                                                 locDATM = "excel",
                                                 locFORCING = "excel",
                                                 readAllForcings = F)

## Optional: change sediment settings
# lDATM_SETTINGS$params <- adjustSedimentParamSettings(lDATM_SETTINGS$params, paramset = 2, sediment_type = "clay")


## 3. Make and adjust cpp files
##    - nRUN_SET determines which forcings are switched on
PCModelAdjustCPPfiles(dirSHELL = dirShell,
                      nameWORKCASE = nameWorkCase,
                      lDATM = lDATM_SETTINGS,
                      nRUN_SET = 0)

## 4. Compile model
PCModelCompileModelWorkCase(dirSHELL = dirShell,
                            nameWORKCASE = nameWorkCase)

# ## 5. Initialize model
# ##    - make all initial states according to the run settings
# InitStates <- PCModelInitializeModel(lDATM = lDATM_SETTINGS,
#                                      dirSHELL = dirShell,
#                                      nameWORKCASE = nameWorkCase)
# 
# ## 6. run one model
# ##    - Error catching on run_state & restart (if run_state = 0 & you use restart should you be able to do so?)
# PCModel_run01 <- PCModelSingleRun(lDATM = lDATM_SETTINGS,
#                                   nRUN_SET = 0,
#                                   dfSTATES = InitStates,
#                                   integrator_method = "rk45ode",
#                                   dirSHELL = dirShell,
#                                   nameWORKCASE = nameWorkCase)

#---BIFURCATIONS ALONG FETCH AND DEPTH----
vDEPTH = seq(0.25, 3.0, by=0.25)
vFETCH = seq(100, 6000, by=300)

dfGRID = expand.grid(vDEPTH, vFETCH)
colnames(dfGRID)=c('depth','fetch')

for(nROW in 1:nrow(dfGRID)){
  fFETCH = dfGRID[nROW,'fetch']
  fDEPTH = dfGRID[nROW,'depth']
  
  ##Cluster computing variant
  library(doSNOW)
  library(foreach)
  #make a cluster for calculations
  vLOAD_MULT = c(seq(0.0001,	0.01, length.out = 21))
  
  nTHREADS=7
  snowCLUSTER <- makeCluster(nTHREADS)
  clusterExport(snowCLUSTER, c())
  registerDoSNOW(snowCLUSTER)
  pb<-txtProgressBar(0,length(vLOAD_MULT),style=3)
  progress<-function(n){
    setTxtProgressBar(pb,n)
  }
  opts<-list(progress=progress)
  
  comb <- function(x, ...) {
    lapply(seq_along(x),
           function(i) c(x[[i]], lapply(list(...), function(y) y[[i]])))
  }
  #vLOAD_MULT
  lOUT <-foreach(i=vLOAD_MULT,.combine='rbind', .multicombine=TRUE,
                 .packages=c('data.table', 'plyr', "dplyr","stringr"),.export = c("PCModelSingleRun"),.options.snow=opts) %dopar% {
                   
                   source(file.path(dirShell, "scripts", "R_system", "functions.R"))  ## load base functions by Luuk van Gerven (2012-2016)
                   source(file.path(dirShell, "scripts", "R_system", "functions_PCLake.R")) 
                   # i <- 0.2
                   lDATM_SETTINGS2=lDATM_SETTINGS
                   #set depth
                   lDATM_SETTINGS2$params$sDefault0[which(rownames(lDATM_SETTINGS$params)=='cDepthWInit0')] = fDEPTH
                   #set fetch
                   lDATM_SETTINGS2$params$sDefault0[which(rownames(lDATM_SETTINGS$params)=='cFetch')] = fFETCH
                   lDATM_SETTINGS2$forcings$sDefault0$mPLoadEpi$value <- i#*lDATM_SETTINGS$forcings$sDefault0$mPLoadEpi$value
                   fPLOAD =lDATM_SETTINGS2$forcings$sDefault0$mPLoadEpi$value[1]
                   sPLOAD = str_replace(as.character(fPLOAD),"\\.","-")
                   
                   print(paste0("This is a base run with pload ", i, "."))
                   
                   ## 5. Initialize model 
                   ##    - make all initial states according to the run settings
                   InitStates <- PCModelInitializeModel(lDATM = lDATM_SETTINGS2,
                                                        dirSHELL = dirShell,
                                                        nameWORKCASE = nameWorkCase)
                   
                   ## 6. run one model - turbid
                   ##    - Error catching on run_state & restart (if run_state = 0 & you use restart should you be able to do so?)
                   PCModel_run01 <- PCModelSingleRun(lDATM = lDATM_SETTINGS2,
                                                     nRUN_SET = 0,
                                                     dfSTATES = InitStates,
                                                     integrator_method = "rk45ck",
                                                     dirSHELL = dirShell,
                                                     nameWORKCASE = nameWorkCase)
                   
                   
                   
                   
                   ## 6. run one model - clear
                   ##    - Error catching on run_state & restart (if run_state = 0 & you use restart should you be able to do so?)
                   PCModel_run02 <- PCModelSingleRun(lDATM = lDATM_SETTINGS2,
                                                     nRUN_SET = 1,
                                                     dfSTATES = InitStates,
                                                     integrator_method = "rk45ck",
                                                     dirSHELL = dirShell,
                                                     nameWORKCASE = nameWorkCase)
                   #results = results
                   
                   # Define output
                   temp_res <- PCModel_run01
                   
                   temp_res$depth <- fDEPTH
                   temp_res$fetch <- fFETCH
                   temp_res$period <- "winter"
                   temp_res[temp_res$time %in% summer_vec, "period"] <- "summer"
                   temp_res$year <- cut(temp_res$time, breaks = seq(91, 30*365, 365), labels = F)
                   temp_res$PLoad=fPLOAD
                   temp_res$initstate="turbid"
                   
                   data.table::fwrite(temp_res, 
                                      file = file.path(dirSave, "work_cases", nameWorkCase, "output","single_runs", paste0("Depth", "_", fDEPTH,"Fetch", "_", fFETCH,"PLoad", "_", sPLOAD,"_turbid.txt")))
                   # Define output
                   temp_res2 <- PCModel_run02
                   
                   temp_res2$depth <- fDEPTH
                   temp_res2$fetch <- fFETCH
                   temp_res2$period <- "winter"
                   temp_res2[temp_res2$time %in% summer_vec, "period"] <- "summer"
                   temp_res2$year <- cut(temp_res2$time, breaks = seq(91, 30*365, 365), labels = F)
                   temp_res2$PLoad=fPLOAD
                   temp_res2$initstate="clear"
                   
                   data.table::fwrite(temp_res2, 
                                      file = file.path(dirSave, "work_cases", nameWorkCase, "output","single_runs", paste0("Depth", "_", fDEPTH,"Fetch", "_", fFETCH,"PLoad", "_", sPLOAD,"_clear.txt")))
                   
                   dfOUT = rbind.data.frame(temp_res,temp_res2)
                   
                   return(dfOUT)
                   
                 }
  stopCluster(snowCLUSTER)
  
  #write output per bifurc
  data.table::fwrite(lOUT, 
                     file = file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", paste0("Depth", "_", fDEPTH,"Fetch", "_", fFETCH,"bifurc_full.txt")))
}

#make last year only files
vFILES = list.files(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs"))#, pattern = "_grid_", ignore.case=TRUE)

#vCOLS = c("time", "sDepthW",     "aDVeg",     "oNH4WEpi",      "oNO3WEpi", "oPO4WEpi",  "oO2WEpi",  "oDDiatWEpi",  "oDGrenWEpi", "oDBlueWEpi",   "oDZooEpi", "oDSestWEpi",   "oPTotWEpi", "oNTotWEpi",   "oChlaEpi",    "int", "depth", "fetch", "period", "year",  "PLoad", "initstate")
vCOLS = c("PLoad", "initstate", "period", "oChlaEpi", "oPTotWEpi", "oNTotWEpi", "aDVeg", "depth", "fetch" )
dfOUT_FULL = data.frame(matrix(NA,0,length(vCOLS),))
colnames(dfOUT_FULL) = vCOLS

for(sFILE in vFILES){
  #sFILE=vFILES[1]
  dtFILE = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", sFILE))
  
  vTIME = c((max(dtFILE$time)-365):max(dtFILE$time))
  
  data.table::fwrite( dtFILE[which(dtFILE$time %in% vTIME),], 
                      file = file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs/last_year", sFILE))
  
  #make total summer average output
  dtFILE_SEL= dtFILE[which(dtFILE$time %in% vTIME),]
  
  dtFILE_SEL[dtFILE_SEL$time %in% summer_vec, "period"] <- "summer"
  output_calc=dtFILE_SEL
  output_calc <- dtFILE_SEL %>% 
    group_by(PLoad, initstate, period) %>%
    summarise(across(colnames(output_calc)[c(grep("oChlaEpi", colnames(output_calc)),grep("oPTotWEpi", colnames(output_calc)),grep("oNTotWEpi", colnames(output_calc)),grep("aDVeg", colnames(output_calc)))], mean))
  
  output_calc_sel = output_calc[which(output_calc$period=="summer"),]
  
  dfOUT_FULL= rbind.data.frame(dfOUT_FULL, cbind(output_calc_sel, depth=as.numeric(str_match(sFILE, "Depth_\\s*(.*?)\\s*Fetch")[,2]), fetch=as.numeric(str_match(sFILE, "Fetch_\\s*(.*?)\\s*bifurc")[,2])))
                     
}

data.table::fwrite(dfOUT_FULL, 
                    file = file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs","last_year", "bifurc_full.txt"))


##CRITICAL LIMITS -----
#extract maximal change for each (i.e. critical limit) for each bifurcation run (both clear and turbid)
#make combined last year grid files filtering
dfOUT_FULL =data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs","last_year", "bifurc_full.txt"))

#declare the outputs 
dfOUT_CRIT = data.frame(matrix(NA,0,5))
colnames(dfOUT_CRIT)=c('depth','fetch','size_ha', 'kP_clear_turbid','kP_turbid_clear')

vCOLS = c("PLoad","initstate","period","oChlaEpi","oPTotWEpi","oNTotWEpi","aDVeg","depth","fetch")
dfOUT_CRIT_ALL = data.frame(matrix(NA,0,length(vCOLS)))
colnames(dfOUT_CRIT_ALL)=vCOLS

#make grid of all combinations
dfGRID_ALL = expand.grid(unique(dfOUT_FULL$depth), unique(dfOUT_FULL$fetch))

for(nSEL in 1:nrow(dfGRID_ALL)){
  #nSEL=35
  #select data within a specific run of fetch and depth
  dfSEL_CLEAR = dfOUT_FULL[which(dfOUT_FULL$depth == dfGRID_ALL$Var1[nSEL] & dfOUT_FULL$fetch == dfGRID_ALL$Var2[nSEL] & dfOUT_FULL$initstate == "clear"),]
  
  dfSEL_TURBID = dfOUT_FULL[which(dfOUT_FULL$depth == dfGRID_ALL$Var1[nSEL] & dfOUT_FULL$fetch == dfGRID_ALL$Var2[nSEL] & dfOUT_FULL$initstate == "turbid"),] 
  #for each run derive the critical limit as the maximum change in Chla or plants 
  vCLEAR_CHLA_DIFF = diff(dfSEL_CLEAR$oChlaEpi)
  vCLEAR_VEG_DIFF = diff(dfSEL_CLEAR$aDVeg)
  fMAX_DIFF_VEG_CLEAR = min(diff(dfSEL_CLEAR$aDVeg))
  fMAX_DIFF_CHLA_CLEAR = max(diff(dfSEL_CLEAR$oChlaEpi))
  
  vTURBID_CHLA_DIFF = diff(dfSEL_TURBID$oChlaEpi)
  vTURBID_VEG_DIFF = diff(dfSEL_TURBID$aDVeg)
  fMAX_DIFF_VEG_TURBID = min(diff(dfSEL_TURBID$aDVeg))
  fMAX_DIFF_CHLA_TURBID = max(diff(dfSEL_TURBID$oChlaEpi))
                              
   #we assume that it lake turbid when chl-a high (right-side loading just after the maximal change in oChla) [clear->turbid critical limit]
  if(fMAX_DIFF_VEG_CLEAR <0){
    dfSEL_CLEAR_CRIT = dfSEL_CLEAR[min(which(vCLEAR_CHLA_DIFF == fMAX_DIFF_CHLA_CLEAR)+1, nrow(dfSEL_CLEAR)),]
    dfOUT_CRIT_ALL = rbind.data.frame(dfOUT_CRIT_ALL, dfSEL_CLEAR_CRIT)
    vCRIT_CLEAR = dfSEL_CLEAR_CRIT$PLoad
  }else{
    vCRIT_CLEAR = NA
  }
  #we assume that the lake is clear when plants have re-emerged after coming from a turbid situation (left-hand side of the maximal change in vegetation biomass)
  if(fMAX_DIFF_VEG_TURBID <0){
    dfSEL_TURBID_CRIT = dfSEL_TURBID[max(which(vTURBID_VEG_DIFF == fMAX_DIFF_VEG_TURBID)-1,1),]
    dfOUT_CRIT_ALL = rbind.data.frame(dfOUT_CRIT_ALL, dfSEL_TURBID_CRIT)
    vCRIT_TURBID = dfSEL_TURBID_CRIT$PLoad
  }else{
    vCRIT_TURBID = NA
  }
  
  dfOUT_CRIT = rbind.data.frame(dfOUT_CRIT, cbind.data.frame(depth=dfGRID_ALL$Var1[nSEL], fetch=dfGRID_ALL$Var2[nSEL], size_ha=(dfGRID_ALL$Var2[nSEL]^2)/10000 , kP_clear_turbid=vCRIT_CLEAR,kP_turbid_clear=vCRIT_TURBID))
  
}


##grid plot of width of hysteresis field
dfOUT_CRIT$hyst_width = dfOUT_CRIT$kP_clear_turbid - dfOUT_CRIT$kP_turbid_clear
dfOUT_CRIT$hyst_width[dfOUT_CRIT$hyst_width<0]=0
dtPLOT= dfOUT_CRIT

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","Hyst_width.pdf"),width=6*1.4, height=5*1.4, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(format(round(size_ha*0.01,2), nsmall = 2)),as.factor(depth))) +
    geom_tile(aes(fill = hyst_width))+
    #geom_raster(aes(fill = hyst_width), interpolate = TRUE)+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("lightyellow", "tan", "darkblue"), values = scales::rescale(c(0,0.35, 1)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("size (ha)")+
    ylab("depth (m)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_classic()+
    theme(axis.text = element_text(size=18),
          axis.title = element_text(size=18),
          legend.text=element_text(size=18),
          legend.title=element_text(size=18),
      legend.position = 'top', legend.key.width = unit(2, 'cm'))+ #change legend key width)+
    labs(x=expression(paste('Size (km'^2,')')), 
         y=expression(paste('Depth (m)')),fill=expression(paste('Width of the hysteretic zone (gP m'^-2,' d'^-1,')')))+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()  

##PTot difference ----

#declare the outputs 

#declare the outputs 
vCOLS = c("PLoad","initstate","period","oChlaEpi","oPTotWEpi","oNTotWEpi","aDVeg","depth","fetch")
dfOUT_CRIT_ALL_PTOT = data.frame(matrix(NA,0,length(vCOLS)))
colnames(dfOUT_CRIT_ALL_PTOT)=vCOLS

dfOUT_CRIT_PTOT = data.frame(matrix(NA,0,5))
colnames(dfOUT_CRIT_PTOT)=c('depth','fetch','size_ha', 'kP_clear_turbid','kP_turbid_clear')

#make grid of all combinations
dfGRID_ALL_PTOT = expand.grid(unique(dfOUT_FULL$depth), unique(dfOUT_FULL$fetch))

for(nSEL in 1:nrow(dfGRID_ALL_PTOT)){
  #nSEL=35
  #select data within a specific run of fetch and depth
  dfSEL_CLEAR = dfOUT_FULL[which(dfOUT_FULL$depth == dfGRID_ALL$Var1[nSEL] & dfOUT_FULL$fetch == dfGRID_ALL_PTOT$Var2[nSEL] & dfOUT_FULL$initstate == "clear"),]
  
  dfSEL_TURBID = dfOUT_FULL[which(dfOUT_FULL$depth == dfGRID_ALL_PTOT$Var1[nSEL] & dfOUT_FULL$fetch == dfGRID_ALL_PTOT$Var2[nSEL] & dfOUT_FULL$initstate == "turbid"),] 
  #for each run derive the critical limit as the maximum change in Chla or plants 
  vCLEAR_CHLA_DIFF = diff(dfSEL_CLEAR$oChlaEpi)
  vCLEAR_VEG_DIFF = diff(dfSEL_CLEAR$aDVeg)
  fMAX_DIFF_VEG_CLEAR = min(diff(dfSEL_CLEAR$aDVeg))
  fMAX_DIFF_CHLA_CLEAR = max(diff(dfSEL_CLEAR$oChlaEpi))
  
  vTURBID_CHLA_DIFF = diff(dfSEL_TURBID$oChlaEpi)
  vTURBID_VEG_DIFF = diff(dfSEL_TURBID$aDVeg)
  fMAX_DIFF_VEG_TURBID = min(diff(dfSEL_TURBID$aDVeg))
  fMAX_DIFF_CHLA_TURBID = max(diff(dfSEL_TURBID$oChlaEpi))
  
  #we assume that it lake turbid when chl-a high (right-side loading just after the maximal change in oChla) [clear->turbid critical limit]
  if(fMAX_DIFF_VEG_CLEAR <0){
    dfSEL_CLEAR_CRIT = dfSEL_CLEAR[min(which(vCLEAR_CHLA_DIFF == fMAX_DIFF_CHLA_CLEAR)+1, nrow(dfSEL_CLEAR)),]
    dfOUT_CRIT_ALL = rbind.data.frame(dfOUT_CRIT_ALL_PTOT, dfSEL_CLEAR_CRIT)
    vCRIT_CLEAR = dfSEL_CLEAR_CRIT$oPTotWEpi
  }else{
    vCRIT_CLEAR = NA
  }
  #we assume that the lake is clear when plants have re-emerged after coming from a turbid situation (left-hand side of the maximal change in vegetation biomass)
  if(fMAX_DIFF_VEG_TURBID <0){
    dfSEL_TURBID_CRIT = dfSEL_TURBID[max(which(vTURBID_VEG_DIFF == fMAX_DIFF_VEG_TURBID)-1,1),]
    dfOUT_CRIT_ALL = rbind.data.frame(dfOUT_CRIT_ALL_PTOT, dfSEL_TURBID_CRIT)
    vCRIT_TURBID = dfSEL_TURBID_CRIT$oPTotWEpi
  }else{
    vCRIT_TURBID = NA
  }
  
  dfOUT_CRIT_PTOT = rbind.data.frame(dfOUT_CRIT_PTOT, cbind.data.frame(depth=dfGRID_ALL$Var1[nSEL], fetch=dfGRID_ALL$Var2[nSEL], size_ha=(dfGRID_ALL$Var2[nSEL]^2)/10000 , kP_clear_turbid=vCRIT_CLEAR,kP_turbid_clear=vCRIT_TURBID))
  
}


##grid plot of width of hysteresis field
dfOUT_CRIT_PTOT$hyst_width = dfOUT_CRIT_PTOT$kP_clear_turbid - dfOUT_CRIT_PTOT$kP_turbid_clear
dfOUT_CRIT_PTOT$hyst_width[dfOUT_CRIT$hyst_width<0]=0
dtPLOT= dfOUT_CRIT_PTOT

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","Hyst_width_PTot.pdf"),width=6*1.4, height=5*1.4, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(format(round(size_ha*0.01,2), nsmall = 2)),as.factor(depth))) +
    geom_tile(aes(fill = hyst_width*1000))+
    #geom_raster(aes(fill = hyst_width), interpolate = TRUE)+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("lightyellow", "tan", "darkblue"), values = scales::rescale(c(0,0.35, 1)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("size (ha)")+
    ylab("depth (m)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_classic()+
    theme(axis.text = element_text(size=18),
          axis.title = element_text(size=18),
          legend.text=element_text(size=18),
          legend.title=element_text(size=18),
          legend.position = 'top', legend.key.width = unit(2, 'cm'))+ #change legend key width)+
    labs(x=expression(paste('Size (km'^2,')')), 
         y=expression(paste('Depth (m)')),fill=expression(paste('Width of the hysteretic zone (', mu,'gP L'^-1,')')))+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()  


#density plots
dfPLOT=dfOUT_FULL

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","density_plot2.pdf"),width=5, height=4, useDingbats=FALSE)
print(
ggplot(dfPLOT, aes(x=oChlaEpi))+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
  theme_classic()+ 
  theme(axis.text = element_text(size=13),
        axis.title = element_text(size=13),
        legend.text=element_text(size=13),
        legend.title=element_text(size=13))+
  xlab(expression(paste('Chl a (', mu,'g L'^-1,')')))+
  ylab("Density")
)
dev.off()

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","density_plot_veg.pdf"),width=5, height=4, useDingbats=FALSE)
print(
  ggplot(dfPLOT, aes(x=ifelse(aDVeg/0.45>100,100,aDVeg/0.45)))+
    geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
    theme_classic()+ 
    theme(axis.text = element_text(size=13),
          axis.title = element_text(size=13),
          legend.text=element_text(size=13),
          legend.title=element_text(size=13))+
    xlab(expression(paste('Vegetation cover (%)')))+
    ylab("Density")
)
dev.off()

ggplot(dfPLOT, aes(x=oPTotWEpi))+
  geom_density(aes(x=PLoad*20), fill="grey20", color="grey40", alpha=0.5)+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)+
  theme_bw()+    
  xlab(expression(paste("TP concentrations", " mg/L")))+
  ylab("density")


#----HYDROLAKES-----
shapeHYDROLAKES <- read_sf(dsn = file.path("C:","Users","SvenT","Documents","HYDROLakes","HydroLAKES_polys_v10_shp"), layer = "HydroLAKES_polys_v10")
#readOGR(dsn = file.path("C:","Users","SvenT","Documents","HYDROLakes","HydroLAKES_polys_v10_shp"), layer = "HydroLAKES_polys_v10")

#make list of depth and size combinations in which we expect ample hysteresis
dfCRIT_OCC = dfOUT_CRIT[which(dfOUT_CRIT$hyst_width >0),]
dfCRIT_OCC$depth_round = round(dfCRIT_OCC$depth,1)
dfCRIT_OCC$area_round = round((dfCRIT_OCC$fetch^2)/1000000,2)

#go through all unique areas and compute a step by step increase of fetch
dfOUT_CRIT_SEL = data.frame(matrix(NA,0,2))
colnames(dfOUT_CRIT_SEL) = c("depth_m", "area_km2")
for(fAREA in unique(dfCRIT_OCC$area_round)){
  #fAREA = unique(dfCRIT_OCC$area_round)[4]
  dfCRIT_OCC_SEL = dfCRIT_OCC[dfCRIT_OCC$area_round==fAREA,]
  vDEPTHS = seq(min(dfCRIT_OCC$depth_round), max(dfCRIT_OCC$depth_round), by=0.1)
  dfOUT_CRIT_SEL = rbind.data.frame(dfOUT_CRIT_SEL, cbind.data.frame(depth_m=vDEPTHS, area_km2 = fAREA))
  
}

#go through all depths and extend the matrix with missing values
dfOUT_CRIT_SEL2 = data.frame(matrix(NA,0,2))
colnames(dfOUT_CRIT_SEL2) = c("depth_m", "area_km2")
for(fDEPTH in unique(dfOUT_CRIT_SEL$depth_m)){
  #fDEPTH = unique(dfOUT_CRIT_SEL$depth_m)[4]
  dfCRIT_OCC_SEL = dfOUT_CRIT_SEL[dfOUT_CRIT_SEL$depth_m==fDEPTH,]
  vAREAS = seq(min(dfCRIT_OCC_SEL$area_km2), max(dfCRIT_OCC_SEL$area_km2), by=0.01)
  dfOUT_CRIT_SEL2 = rbind.data.frame(dfOUT_CRIT_SEL2, cbind.data.frame(depth_m=fDEPTH, area_km2 = vAREAS))
  
}
dfOUT_CRIT_SEL2$depth_area_comb = paste0(as.character(dfOUT_CRIT_SEL2$depth_m), "_", as.character(dfOUT_CRIT_SEL2$area_km2))

#select entries in HYDROLAKES that are within the range where there is hysteresis according to PCLake analysis
vHYDROLAKES_DEPTH_AREA_COMB = paste0(as.character(shapeHYDROLAKES$Depth_avg), "_", as.character(shapeHYDROLAKES$Lake_area))
shapeHYDROLAKES_SELECT = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL2$depth_area_comb),]
nLAKES = nrow(shapeHYDROLAKES)
#select different size classes
#<1km2
dfOUT_CRIT_SEL_LE1 = dfOUT_CRIT_SEL2[which(dfOUT_CRIT_SEL2$area_km2<1),]
shapeHYDROLAKES_SELECT_LE1 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_LE1$depth_area_comb),]
fHL_DEPTH_HYST_LE1 =nrow(shapeHYDROLAKES_SELECT_LE1)/length(which(shapeHYDROLAKES$Lake_area<1))*100 
nLAKES_LE1 = length(which(shapeHYDROLAKES$Lake_area<1))
nLAKES_HYST_LE1 = nrow(shapeHYDROLAKES_SELECT_LE1)

#1-10km2
dfOUT_CRIT_SEL_LE10 = dfOUT_CRIT_SEL2[which(dfOUT_CRIT_SEL2$area_km2>=1 & dfOUT_CRIT_SEL2$area_km2<10),]
shapeHYDROLAKES_SELECT_LE10 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_LE10$depth_area_comb),]
fHL_DEPTH_HYST_LE10 =nrow(shapeHYDROLAKES_SELECT_LE10)/length(which(shapeHYDROLAKES$Lake_area>=1 & shapeHYDROLAKES$Lake_area<10))*100 
nLAKES_LE10 = length(which(shapeHYDROLAKES$Lake_area>=1 & shapeHYDROLAKES$Lake_area<10))
nLAKES_HYST_LE10 = nrow(shapeHYDROLAKES_SELECT_LE10)

#10-100
dfOUT_CRIT_SEL_LE100 = dfOUT_CRIT_SEL2[which(dfOUT_CRIT_SEL2$area_km2>=10 & dfOUT_CRIT_SEL2$area_km2<100),]
shapeHYDROLAKES_SELECT_LE100 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_LE100$depth_area_comb),]
fHL_DEPTH_HYST_LE100 =nrow(shapeHYDROLAKES_SELECT_LE100)/length(which(shapeHYDROLAKES$Lake_area>=10 & shapeHYDROLAKES$Lake_area<100))*100 
nLAKES_LE100 = length(which(shapeHYDROLAKES$Lake_area>=10 & shapeHYDROLAKES$Lake_area<100))
nLAKES_HYST_LE100 = nrow(shapeHYDROLAKES_SELECT_LE100)

#100-1000
dfOUT_CRIT_SEL_LE1000 = dfOUT_CRIT_SEL2[which(dfOUT_CRIT_SEL2$area_km2>=100 & dfOUT_CRIT_SEL2$area_km2<1000),]
shapeHYDROLAKES_SELECT_LE1000 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_LE1000$depth_area_comb),]
fHL_DEPTH_HYST_LE1000 =nrow(shapeHYDROLAKES_SELECT_LE1000)/length(which(shapeHYDROLAKES$Lake_area>=100 & shapeHYDROLAKES$Lake_area<1000))*100 
nLAKES_LE1000 = length(which(shapeHYDROLAKES$Lake_area>=100 & shapeHYDROLAKES$Lake_area<1000))
nLAKES_HYST_LE1000 = nrow(shapeHYDROLAKES_SELECT_LE1000)

#1000-10000
dfOUT_CRIT_SEL_LE10000 = dfOUT_CRIT_SEL2[which(dfOUT_CRIT_SEL2$area_km2>=1000 & dfOUT_CRIT_SEL2$area_km2<10000),]
shapeHYDROLAKES_SELECT_LE10000 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_LE10000$depth_area_comb),]
fHL_DEPTH_HYST_LE10000 =nrow(shapeHYDROLAKES_SELECT_LE10000)/length(which(shapeHYDROLAKES$Lake_area>=1000 & shapeHYDROLAKES$Lake_area<10000))*100 
nLAKES_LE10000 = length(which(shapeHYDROLAKES$Lake_area>=1000 & shapeHYDROLAKES$Lake_area<10000))
nLAKES_HYST_LE10000 = nrow(shapeHYDROLAKES_SELECT_LE10000)

#>=10000
dfOUT_CRIT_SEL_GE10000 = dfOUT_CRIT_SEL2[which( dfOUT_CRIT_SEL2$area_km2>=10000),]
shapeHYDROLAKES_SELECT_GE10000 = shapeHYDROLAKES[which(vHYDROLAKES_DEPTH_AREA_COMB %in% dfOUT_CRIT_SEL_GE10000$depth_area_comb),]
fHL_DEPTH_HYST_GE10000 =nrow(shapeHYDROLAKES_SELECT_GE10000)/length(which(shapeHYDROLAKES$Lake_area>10000))*100 
nLAKES_GE10000 = length(which(shapeHYDROLAKES$Lake_area>10000))
nLAKES_HYST_GE10000 = nrow(shapeHYDROLAKES_SELECT_GE10000)

dfHYDROLAKES_SIZE_TABLE = cbind.data.frame(size_class = c("<1","1-10","10-100","100-1000","1000-10000",">=10000"), no_lakes=c(nLAKES_LE1,nLAKES_LE10,nLAKES_LE100,nLAKES_LE1000,nLAKES_LE10000,nLAKES_GE10000), no_lakes_with_hysteresis=c(nLAKES_HYST_LE1,nLAKES_HYST_LE10,nLAKES_HYST_LE100,nLAKES_HYST_LE1000,nLAKES_HYST_LE10000,nLAKES_HYST_GE10000))

write.csv(dfHYDROLAKES_SIZE_TABLE,file.path(dirSave, "work_cases", nameWorkCase, "output","HYDROLAKES_SIZE_HYST.csv"))
#percentage
nrow(shapeHYDROLAKES_SELECT)/nrow(shapeHYDROLAKES)*100 #=43.86834%


#load land cover data from HYDROAtlas
shapeHYDROATLAS_WEST <- read_sf(dsn = file.path("C:","Users","SvenT","Documents","HYDROLakes","LakeATLAS_v10_shp"), layer = "LakeATLAS_v10_pol_west")
shapeHYDROATLAS_EAST <- read_sf(dsn = file.path("C:","Users","SvenT","Documents","HYDROLakes","LakeATLAS_v10_shp"), layer = "LakeATLAS_v10_pol_east")
shapeHYDROLAKES_SELECT$Hylak_id

#select relevant shapes in HYDROATLAS
shapeHYDROATLAS_WEST_SEL = shapeHYDROATLAS_WEST[which(shapeHYDROATLAS_WEST$Hylak_id %in% shapeHYDROLAKES_SELECT$Hylak_id),]

shapeHYDROATLAS_EAST_SEL = shapeHYDROATLAS_EAST[which(shapeHYDROATLAS_EAST$Hylak_id %in% shapeHYDROLAKES_SELECT$Hylak_id),]

library(readxl)
#read in excel data base with legends for land use
dfHYDROATLAS_GLC_LEGEND <- read_excel(file.path("C:","Users","SvenT","Documents","HYDROLakes", "HydroATLAS_v10_Legends.xlsx"), sheet = "glc_cl")

#glc_pc_uXX is in percentage of watershed area
#We assume a watershed is heavily modified when it has more than 20% non-natural landscape area (16 Cultivated and managed areas, 17 Mosaic: cropland/tree cover/other natural vegetation 18 Mosaic: cropland/shrub and/or herbaceous cover 22 Artificial surfaces and associated areas)
## while correcting for the area of water in the watershed (20 Water bodies (natural and artificial))
shapeHYDROATLAS_EAST_SEL$glc_pc_u_art = ((shapeHYDROATLAS_EAST_SEL$glc_pc_u14+shapeHYDROATLAS_EAST_SEL$glc_pc_u16+shapeHYDROATLAS_EAST_SEL$glc_pc_u17+shapeHYDROATLAS_EAST_SEL$glc_pc_u18+shapeHYDROATLAS_EAST_SEL$glc_pc_u19+shapeHYDROATLAS_EAST_SEL$glc_pc_u22)/(100-shapeHYDROATLAS_EAST_SEL$glc_pc_u20))*100

shapeHYDROATLAS_WEST_SEL$glc_pc_u_art = ((shapeHYDROATLAS_WEST_SEL$glc_pc_u14+shapeHYDROATLAS_WEST_SEL$glc_pc_u16+shapeHYDROATLAS_WEST_SEL$glc_pc_u17+shapeHYDROATLAS_WEST_SEL$glc_pc_u18+shapeHYDROATLAS_WEST_SEL$glc_pc_u19+ shapeHYDROATLAS_WEST_SEL$glc_pc_u22)/(100-shapeHYDROATLAS_WEST_SEL$glc_pc_u20))*100

shapeHYDROATLAS_SEL = rbind(shapeHYDROATLAS_WEST_SEL, shapeHYDROATLAS_EAST_SEL)

#percentage in heavily modified or disturbed landscapes:
(sum((shapeHYDROATLAS_SEL$glc_pc_u_art > 25)+0, na.rm=TRUE)/sum(is.na(shapeHYDROATLAS_SEL$glc_pc_u_art)==FALSE, na.rm=TRUE))*100#=49.91339% 



#shapeHYDROATLAS_EAST_SEL$glc_pc_u_all = shapeHYDROATLAS_EAST_SEL$glc_pc_u01+shapeHYDROATLAS_EAST_SEL$glc_pc_u02+shapeHYDROATLAS_EAST_SEL$glc_pc_u03+shapeHYDROATLAS_EAST_SEL$glc_pc_u04+shapeHYDROATLAS_EAST_SEL$glc_pc_u05+shapeHYDROATLAS_EAST_SEL$glc_pc_u06+shapeHYDROATLAS_EAST_SEL$glc_pc_u07+shapeHYDROATLAS_EAST_SEL$glc_pc_u08+shapeHYDROATLAS_EAST_SEL$glc_pc_u09+shapeHYDROATLAS_EAST_SEL$glc_pc_u10+shapeHYDROATLAS_EAST_SEL$glc_pc_u11+shapeHYDROATLAS_EAST_SEL$glc_pc_u12+shapeHYDROATLAS_EAST_SEL$glc_pc_u13+shapeHYDROATLAS_EAST_SEL$glc_pc_u14+shapeHYDROATLAS_EAST_SEL$glc_pc_u15+shapeHYDROATLAS_EAST_SEL$glc_pc_u16+shapeHYDROATLAS_EAST_SEL$glc_pc_u17+shapeHYDROATLAS_EAST_SEL$glc_pc_u18+shapeHYDROATLAS_EAST_SEL$glc_pc_u19+shapeHYDROATLAS_EAST_SEL$glc_pc_u20+shapeHYDROATLAS_EAST_SEL$glc_pc_u21+shapeHYDROATLAS_EAST_SEL$glc_pc_u22


table(shapeHYDROATLAS_WEST_SEL$glc_cl_vmj)+table(shapeHYDROATLAS_EAST_SEL$glc_cl_vmj) 

dfLAND_USE_FREQ = merge(dfHYDROATLAS_GLC_LEGEND, as.data.frame(table(shapeHYDROATLAS_WEST_SEL$glc_cl_vmj)), by.x="GLC_ID", by.y="Var1", keep.x=TRUE)

dfLAND_USE_FREQ = merge(dfLAND_USE_FREQ, as.data.frame(table(shapeHYDROATLAS_EAST_SEL$glc_cl_vmj)), by.x="GLC_ID", by.y="Var1", keep.x=TRUE)
dfLAND_USE_FREQ$count = dfLAND_USE_FREQ$Freq.x+dfLAND_USE_FREQ$Freq.y
dfLAND_USE_FREQ$perc = dfLAND_USE_FREQ$count/sum(dfLAND_USE_FREQ$count)*100

dfHYDROATLAS_GLC_LEGEND$GLC_Name

seq(min(dfCRIT_OCC$area_round), max(dfCRIT_OCC$area_round), by=0.01)

shapeHYDROLAKES$Depth_avg
shapeHYDROLAKES$Lake_area





dfPLOT$macrofyte_dom = (dfPLOT$aDVeg > 30)+0


vSEL = sample(1:nrow(dfPLOT), 1000)

ggplot(dfPLOT, aes(x=log10(oPTotWEpi), y=log10(oChlaEpi)))+
  geom_point()  + facet_wrap(~macrofyte_dom, scales="free")
ggplot(dfPLOT[vSEL,], aes(x=oPTotWEpi, y=oChlaEpi))+
  geom_point()+

ggplot(dfPLOT, aes(x=oChlaEpi))+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)

ggplot(dfPLOT, aes(x=log10(aDVeg)))+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8)

ggplot(dfPLOT, aes(x=oPTotWEpi))+
  geom_histogram()


ggplot(dfPLOT, aes(x=oPTotWEpi, y=oChlaEpi))+
  geom_point()

#extract maximal change for each (i.e. critical limit) for each bifurcation run (both clear and turbid)
#make combined last year grid files filtering
vFILES = list.files(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year"), ignore.case=TRUE)
sFILE=vFILES[1]
dtFILE = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year", sFILE))
dfOUT = data.frame(matrix(NA,0,6))
colnames(dfOUT)=c('depth','fetch','diff_oChla', 'diff_aDVeg','diff_oPTot','diff_oNTot')

for(sFILE in vFILES){
    dtFILE = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year", sFILE))
  
    #select summer period conform the definition of davidson (May-September)
    dtFILE$day = dtFILE$time-(365*29)
    nDAY_START =121
    nDAY_END = 273
    dtFILE_SEL = dtFILE[which(dtFILE$day<=nDAY_END & dtFILE$day>=nDAY_START),]
    
    output_calc=dtFILE_SEL
    output_calc <- dtFILE_SEL %>% 
      group_by(PLoad, initstate) %>%
      summarise(across(colnames(output_calc)[c(grep("oChlaEpi", colnames(output_calc)),grep("oPTotWEpi", colnames(output_calc)),grep("oNTotWEpi", colnames(output_calc)),grep("aDVeg", colnames(output_calc)))], mean))

    dfDIFF=cbind.data.frame(diff_oChla=diff(output_calc$oChlaEpi)[seq(1,nrow(output_calc), by=2)],diff_aDVeg=diff(output_calc$aDVeg)[seq(1,nrow(output_calc), by=2)],diff_oNTot=diff(output_calc$oNTotWEpi)[seq(1,nrow(output_calc), by=2)],diff_oPTot=diff(output_calc$oPTotWEpi)[seq(1,nrow(output_calc), by=2)],PLoad=output_calc$PLoad[seq(1,nrow(output_calc), by=2)])
    nMAX_CHLA = which(dfDIFF$diff_oChla == max(dfDIFF$diff_oChla))

    dfOUT = rbind.data.frame(dfOUT, cbind.data.frame(depth = dtFILE$depth[1], fetch=dtFILE$fetch[1], dfDIFF[nMAX_CHLA,c('diff_oChla', 'diff_aDVeg','diff_oPTot','diff_oNTot')]))
}

dtPLOT = dfOUT

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","FD_TPdiff.pdf"),width=6, height=5, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(fetch),as.factor(depth))) +
    geom_tile(aes(fill = diff_oPTot))+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("lightblue4", "maroon1", "maroon4"), values = scales::rescale(c(0,0.1, 0.5)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("fetch (m)")+
    ylab("depth (m)")+
    labs(fill = "Difference in TP in hysteretic zone (mg/L)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_bw()+
    theme(legend.position="top")+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()    


pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","FD_TNdiff.pdf"),width=6, height=5, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(fetch),as.factor(depth))) +
    geom_tile(aes(fill = diff_oNTot))+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("lightblue4", "royalblue", "navy"), values = scales::rescale(c(0,0.1, 0.5)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("fetch (m)")+
    ylab("depth (m)")+
    labs(fill = "Difference in TN in hysteretic zone (mg/L)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_bw()+
    theme(legend.position="top")+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()  

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","FD_Chladiff.pdf"),width=6, height=5, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(fetch),as.factor(depth))) +
    geom_tile(aes(fill = diff_oChla))+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("lightyellow", "seagreen1", "seagreen4"), values = scales::rescale(c(0,100, 300)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("fetch (m)")+
    ylab("depth (m)")+
    labs(fill = "Difference in chlorophyll-a in hysteretic zone (µg/L)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_bw()+
    theme(legend.position="top")+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()  

pdf(file.path(dirSave, "work_cases", nameWorkCase, "output/grid_plots","FD_Vegdiff.pdf"),width=6, height=5, useDingbats=FALSE)
print(
  ggplot(dtPLOT, aes(as.factor(fetch),as.factor(depth))) +
    geom_tile(aes(fill = diff_aDVeg))+
    #geom_rect(aes(xmin = 0.5, xmax = 15.5 , ymin = -0.5, ymax = 0.5),fill=NA, color="darkred", size=1.0)+      
    #coord_flip()+
    scale_fill_gradientn(colours = c("darkred", "lightyellow", "darkgreen"), values = scales::rescale(c(-400,0,100)), limits=c(NA,NA))+
    #facet_wrap(~ initstate, nrow=2)+
    xlab("fetch (m)")+
    ylab("depth (m)")+
    labs(fill = "Difference in vegetation in hysteretic zone (gDW)")+
    #scale_y_reverse(breaks=seq(0,-10,-2))+
    scale_x_discrete(guide = guide_axis(angle = 90)) +
    theme_bw()+
    theme(legend.position="top")+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
)
dev.off()  


###PLOTS CNP LOAD----
dtPLOT_FULL = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs","last_year","bifurc_full.txt"))

output_calc=dtPLOT_FULL
output_calc <- dtPLOT_FULL %>% 
  group_by(PLoad, initstate) %>%
  summarise(across(colnames(output_calc)[2:75], mean))

dtPLOT = output_calc
dfPLOT_BIF_SEL = as.data.frame(dtPLOT)

for(sVAR in colnames(dtPLOT)[3:76]){
 #sVAR = colnames(dtPLOT)[4] 
  plot_bifurc <- ggplot(data = dfPLOT_BIF_SEL, aes(x=PLoad*1000, y=get(sVAR), color=initstate)) +
    geom_point(aes(shape=initstate))+
    geom_line(size = 0.5, aes(linetype=initstate))+
    scale_y_continuous(limits = c(min(dfPLOT_BIF_SEL[sVAR],0), NA)) +
    scale_x_continuous(limits = c(0, NA)) +  
    #scale_x_date(breaks = date_breaks("3 months"), labels = date_format("%b"))+
    scale_color_manual(values=c("green4", "tan4"))+
    ylab(sVAR) +
    xlab("Nutrient loading (mgP/m²/dag)") + 
    theme_bw()+
    theme(legend.position="top", strip.background = element_rect(fill = "white", colour = "white"))+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
  
  pdf(file.path(dirShell, "work_cases", nameWorkCase, "output","plots",paste0(sVAR,"_bifur.pdf")),width=6, height=4, useDingbats=FALSE)
  print(plot_bifurc)
  dev.off()
}

##PLoad/Pconc vs oChla -----


##CRITICAL LIMITS -----
#extract maximal change for each (i.e. critical limit) for each bifurcation run (both clear and turbid)
#make combined last year grid files filtering
vFILES = list.files(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year"), ignore.case=TRUE)
sFILE=vFILES[1]
dtFILE = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year", sFILE))
dfOUT = data.frame(matrix(NA,0,6))
colnames(dfOUT)=c('depth','fetch','size_ha', 'kPCT_Plant','kPTC_Plant','kPTC_Phyt','kPCT_Phyt')

for(sFILE in vFILES){
  dtFILE = data.table::fread(file.path(dirSave, "work_cases", nameWorkCase, "output","bifurc_runs", "last_year", sFILE))
  
  #select summer period conform the definition of davidson (May-September)
  dtFILE$day = dtFILE$time-(365*29)
  nDAY_START =121
  nDAY_END = 273
  dtFILE_SEL = dtFILE[which(dtFILE$day<=nDAY_END & dtFILE$day>=nDAY_START),]
  
  output_calc=dtFILE_SEL
  output_calc <- dtFILE_SEL %>% 
    group_by(PLoad, initstate) %>%
    summarise(across(colnames(output_calc)[c(grep("oChlaEpi", colnames(output_calc)),grep("oPTotWEpi", colnames(output_calc)),grep("oNTotWEpi", colnames(output_calc)),grep("aDVeg", colnames(output_calc)))], mean))
  
  dfDIFF=cbind.data.frame(diff_oChla=diff(output_calc$oChlaEpi)[seq(1,nrow(output_calc), by=2)],diff_aDVeg=diff(output_calc$aDVeg)[seq(1,nrow(output_calc), by=2)],diff_oNTot=diff(output_calc$oNTotWEpi)[seq(1,nrow(output_calc), by=2)],diff_oPTot=diff(output_calc$oPTotWEpi)[seq(1,nrow(output_calc), by=2)],PLoad=output_calc$PLoad[seq(1,nrow(output_calc), by=2)])
  nMAX_CHLA = which(dfDIFF$diff_oChla == max(dfDIFF$diff_oChla))
  
  dfOUT = rbind.data.frame(dfOUT, cbind.data.frame(depth = dtFILE$depth[1], fetch=dtFILE$fetch[1], dfDIFF[nMAX_CHLA,c('diff_oChla', 'diff_aDVeg','diff_oPTot','diff_oNTot')]))
}

dtPLOT = dfOUT








#---EUTROPHICATION: NP LOAD ONLY----
##Cluster computing variant
library(doSNOW)
library(foreach)
#make a cluster for calculations
vLOAD_MULT = c(seq(0.0001,	0.01, length.out = 21))
lDATM_SETTINGS$params[which(rownames(lDATM_SETTINGS$params)=='cCPLoadMeas'),'sDefault0']=0 #normally, turn on ReadCLoad, then set mCLoad to zero, but that means making a whole new forcing which is a bit off a hassle

nTHREADS=7
snowCLUSTER <- makeCluster(nTHREADS)
clusterExport(snowCLUSTER, c())
registerDoSNOW(snowCLUSTER)
pb<-txtProgressBar(0,length(vLOAD_MULT),style=3)
progress<-function(n){
  setTxtProgressBar(pb,n)
}
opts<-list(progress=progress)

comb <- function(x, ...) {
  lapply(seq_along(x),
         function(i) c(x[[i]], lapply(list(...), function(y) y[[i]])))
}
#vLOAD_MULT
lOUT <-foreach(i=vLOAD_MULT,.combine='rbind', .multicombine=TRUE,
               .packages=c('data.table', 'plyr', "dplyr","stringr"),.export = c("PCModelSingleRun"),.options.snow=opts) %dopar% {
                 
                 source(file.path(dirShell, "scripts", "R_system", "functions.R"))  ## load base functions by Luuk van Gerven (2012-2016)
                 source(file.path(dirShell, "scripts", "R_system", "functions_PCLake.R")) 
                 # i <- 0.2
                 lDATM_SETTINGS2=lDATM_SETTINGS
                 lDATM_SETTINGS2$forcings$sDefault0$mPLoadEpi$value <- i#*lDATM_SETTINGS$forcings$sDefault0$mPLoadEpi$value
                 fPLOAD =lDATM_SETTINGS2$forcings$sDefault0$mPLoadEpi$value[1]
                 sPLOAD = str_replace(as.character(fPLOAD),"\\.","-")
                 
                 print(paste0("This is a base run with pload ", i, "."))
                 
                 ## 5. Initialize model 
                 ##    - make all initial states according to the run settings
                 InitStates <- PCModelInitializeModel(lDATM = lDATM_SETTINGS2,
                                                      dirSHELL = dirShell,
                                                      nameWORKCASE = nameWorkCase)
                 
                 ## 6. run one model - turbid
                 ##    - Error catching on run_state & restart (if run_state = 0 & you use restart should you be able to do so?)
                 PCModel_run01 <- PCModelSingleRun(lDATM = lDATM_SETTINGS2,
                                                   nRUN_SET = 0,
                                                   dfSTATES = InitStates,
                                                   integrator_method = "rk45ck",
                                                   dirSHELL = dirShell,
                                                   nameWORKCASE = nameWorkCase)
                 
                 
                 
                 
                 ## 6. run one model - clear
                 ##    - Error catching on run_state & restart (if run_state = 0 & you use restart should you be able to do so?)
                 PCModel_run02 <- PCModelSingleRun(lDATM = lDATM_SETTINGS2,
                                                   nRUN_SET = 1,
                                                   dfSTATES = InitStates,
                                                   integrator_method = "rk45ck",
                                                   dirSHELL = dirShell,
                                                   nameWORKCASE = nameWorkCase)
                 #results = results
                 
                 # Define output
                 temp_res <- PCModel_run01
                 
                 temp_res$period <- "winter"
                 temp_res[temp_res$time %in% summer_vec, "period"] <- "summer"
                 temp_res$year <- cut(temp_res$time, breaks = seq(91, 100*365, 365), labels = F)
                 temp_res$PLoad=fPLOAD
                 temp_res$initstate="turbid"
                 
                 data.table::fwrite(temp_res, 
                                    file = file.path(dirShell, "work_cases", nameWorkCase, "output","single_runs", paste0("PLoad", "_", sPLOAD,"_turbidNP.txt")))
                 # Define output
                 temp_res2 <- PCModel_run02
                 
                 temp_res2$period <- "winter"
                 temp_res2[temp_res2$time %in% summer_vec, "period"] <- "summer"
                 temp_res2$year <- cut(temp_res2$time, breaks = seq(91, 100*365, 365), labels = F)
                 temp_res2$PLoad=fPLOAD
                 temp_res2$initstate="clear"
                 
                 data.table::fwrite(temp_res2, 
                                    file = file.path(dirShell, "work_cases", nameWorkCase, "output","single_runs", paste0("PLoad", "_", sPLOAD,"_clearNP.txt")))
                 
                 dfOUT = rbind.data.frame(temp_res,temp_res2)
                 
                 return(dfOUT)
                 
               }
stopCluster(snowCLUSTER)

#write output per bifurc
data.table::fwrite(lOUT, 
                   file = file.path(dirShell, "work_cases", nameWorkCase, "output","bifurc_runs", paste0("bifurc_fullNP.txt")))


#make last year only files
vFILES = list.files(file.path(dirShell, "work_cases", nameWorkCase, "output","bifurc_runs"))#, pattern = "_grid_", ignore.case=TRUE)

for(sFILE in vFILES){
  #sFILE=vFILES[1]
  dtFILE = data.table::fread(file.path(dirShell, "work_cases", nameWorkCase, "output","bifurc_runs", sFILE))
  
  vTIME = c((max(dtFILE$time)-365):max(dtFILE$time))
  
  data.table::fwrite( dtFILE[which(dtFILE$time %in% vTIME),], 
                      file = file.path(dirShell, "work_cases", nameWorkCase, "output","bifurc_runs/last_year", sFILE))
  
}

###PLOTS NP LOAD----

dtPLOT_FULL = data.table::fread(file.path(dirShell, "work_cases", nameWorkCase, "output","bifurc_runs","last_year","bifurc_fullNP.txt"))

output_calc=dtPLOT_FULL
output_calc <- dtPLOT_FULL %>% 
  group_by(PLoad, initstate) %>%
  summarise(across(colnames(output_calc)[2:75], mean))

dtPLOT = output_calc
dfPLOT_BIF_SEL = as.data.frame(dtPLOT)

for(sVAR in colnames(dtPLOT)[3:76]){
  #sVAR = colnames(dtPLOT)[4] 
  plot_bifurc <- ggplot(data = dfPLOT_BIF_SEL, aes(x=PLoad*1000, y=get(sVAR), color=initstate)) +
    geom_point(aes(shape=initstate))+
    geom_line(size = 0.5, aes(linetype=initstate))+
    scale_y_continuous(limits = c(min(dfPLOT_BIF_SEL[sVAR],0), NA)) +
    scale_x_continuous(limits = c(0, NA)) +  
    #scale_x_date(breaks = date_breaks("3 months"), labels = date_format("%b"))+
    scale_color_manual(values=c("green4", "tan4"))+
    ylab(sVAR) +
    xlab("Nutrient loading (mgP/m²/dag)") + 
    theme_bw()+
    theme(legend.position="top", strip.background = element_rect(fill = "white", colour = "white"))+
    guides(fill = guide_colourbar(title.position="top", title.hjust = 0.5))
  
  pdf(file.path(dirShell, "work_cases", nameWorkCase, "output","plots",paste0(sVAR,"_bifurNP.pdf")),width=6, height=4, useDingbats=FALSE)
  print(plot_bifurc)
  dev.off()
}