###################################################
## Prepare input data for Hydrus 1D simulations ###
###################################################

# load R libraries
library(zoo); Sys.setenv(TZ = "GMT")
library(xts)

#####################
## set directory for output of Hydrus 1D input files
dir_output = "testdir/"
#####################

#####################
## necessary function to properly convert from zoo-object to data.frame
#####################
# has to be run in order to make the function known for later usage (further down)
zootodf = function(value) {
    df = as.data.frame(value)
    df$time = index(value) #create a Date column
    rownames(df) = NULL #so row names not filled with dates
    df = df[,c(ncol(df), 1:(ncol(df)-1))] #reorder columns so Date first
    return(df)
}

#####################
## generate synthetic input data
#####################
# complete artificial setup
# serves as virtual experiment to verify theoretical umbrella effect and topographic effect on gravity signal
# following structure:
# 15 days: only rain (const)
# 15 days: only EVT (const)
## precipitation
inf.start = as.POSIXct("2015-01-17 00:00:00 GMT") 
inf.end = as.POSIXct("2015-02-01 00:00:00 GMT") 
precip_synt = zoo(1, order.by=seq(inf.start,inf.end, by="hour"))
## evt
evt.start = as.POSIXct("2015-02-01 00:00:00 GMT") 
evt.end = as.POSIXct("2015-02-16 00:00:00 GMT") 
# with EVT
evt_synt = zoo(0.1, order.by=seq(evt.start,evt.end, by="hour"))
# without EVT
evt_synt = zoo(0, order.by=seq(evt.start,evt.end, by="hour"))

# combine both time series
timeseries = merge(precip_synt,evt_synt, all=T, fill=0)
timeseries = zootodf(timeseries)
# set column names
colnames(timeseries) = c("Datetime","Precipitation","ET0")

# specify name of output file
# for Hydrus 1D input, this has to be ATMOSPH.IN
file.tsinput.save = paste0(dir_output, "hydrus1d_forcing_InfExperiment_hydrusFormat_withoutEVT.IN")
# file.tsinput.save = paste0(dir_output, "ATMOSPH.IN")

## write hydrus input
# output directly to hydrus ATMOSPH.IN-file
# atmospheric data rounded to 7 digits after coma
# without groundwater !!
numberofdata = length(index(timeseries))
data.out.hydrus = data.frame(tAtm = 1:length(index(timeseries)),Prec = coredata(timeseries$Precipitation)/1000,rSoil =coredata(timeseries$ET0)/1000,rRoot=0,hCritA = 150,rB=0,hB = 0,ht=0, RootDepth=0)

# formating output of numbers
data.out.hydrus_format = data.out.hydrus
data.out.hydrus_format$Prec = format(data.out.hydrus_format$Prec, digits=3)
data.out.hydrus_format$rSoil = format(data.out.hydrus_format$rSoil, digits=3)
data.out.hydrus_format$hB = format(data.out.hydrus_format$hB, digits=3)

# hydrus 1D & (old) 2D
text.start = rbind("Pcp_File_Version=4","*** BLOCK I: ATMOSPHERIC INFORMATION  **********************************","   MaxAL                    (MaxAL = number of atmospheric data-records)",paste("  ",numberofdata,sep="")," DailyVar  SinusVar  lDummy  lDummy  lDummy  lDummy  lDummy  lDummy  lDummy  lDummy","       f       f       f       f       f       f       f       f       f       f"," hCritS                 (max. allowed pressure head at the soil surface)","      0")
text.end = c("end*** END OF INPUT FILE 'ATMOSPH.IN' **********************************")
write.table(text.start, file=file.tsinput.save, row.names=F, col.names=F, quote=F)
write.table(data.out.hydrus_format, file=file.tsinput.save, row.names=F, col.names=T, sep="      ", dec=".", append=T, quote=F)
write.table(text.end, file=file.tsinput.save, row.names=F, col.names=F, append=T, quote=F)

#####################
## generate printout times file
#####################
# the output has to be *manually* inserted in file: Selektor.IN
ptimes_max = 721
pstart = 1
pend = 721
# rounded to whole numbers
ptimes = round(seq(pstart, pend, length.out=ptimes_max),0)
# not rounded
# this should be avoided as hydrus could have problems finding data for this exact timestep
# write data
write.table(ptimes, file=paste0(dir_output, "hydrus1d_completeTS_ptimes.csv"), col.names=F, row.names=F, sep="     ")

