### Read weather types (Schwander et al. 2017), note that this file has the entire year, below we select summer 
y <- read.table("CAP7_1763-2009.txt",header=T)

### read the weights for Basel and Rhone
fb <- read.table("fpi_weights_bas.txt")
fr <- read.table("fpi_weights_rho.txt")

#### weights for time lags
wd <- c(1/16, 2/16, 3/16, 4/16, 4/16, 2/16)

yr <- as.numeric(substr(y[,1],1,4))
mon <- as.numeric(substr(y[,1],6,7))
day <- as.numeric(substr(y[,1],9,10))

### process Rhine/Basel
selseas <- mon>4&mon<11
wt <- y[selseas,2]
yrm <- yr[selseas]
fpidbas <- rep(0,length(wt))
for (i in 6:length(wt)){
  fpidbas[i] <- (wd[1]*fb[1,wt[i-5]])+(wd[2]*fb[2,wt[i-4]])+(wd[3]*fb[3,wt[i-3]])+(wd[4]*fb[4,wt[i-2]])+(wd[5]*fb[5,wt[i-1]])+(wd[6]*fb[6,wt[i]])
}
fpiabas <- matrix(NA,nrow=247,ncol=5)
fpiabas[,1] <- c(1763:2009)
for (i in 1763:2009){
 fpiabas[i-1762,2] <- mean(fpidbas[yrm==i])
 fpiabas[i-1762,3:5] <- quantile(fpidbas[yrm==i],c(0.5,0.75,0.9))
}
#### write daily index
write.table(cbind(yrm,mon[selseas],day[selseas],fpidbas),"fpi_daily_Basel_MJJASO.txt",sep=" ",row.names=F,col.names=c("yr","mon","day","fpi_daily_Basel"))
### write annual index
write.table(fpiabas,"fpi_annual_Basel_MJJASO.txt",sep=" ",row.names=F,col.names=c("yr","mean","50th_percentile","75th_percentile","90th_percentile"))

### process Rhone/Beaucaire
selseas <- mon<5|mon>10
yw <- yr
yw[mon>10] <- yw[mon>10]+1
wt <- y[selseas,2]
yrm <- yw[selseas]
fpidrho <- rep(0,length(wt))
for (i in 6:length(wt)){
  fpidrho[i] <- (wd[1]*fr[1,wt[i-5]])+(wd[2]*fr[2,wt[i-4]])+(wd[3]*fr[1,wt[i-3]])+(wd[4]*fr[4,wt[i-2]])+(wd[5]*fr[5,wt[i-1]])+(wd[6]*fr[6,wt[i]])
}

fpiarho <- matrix(NA,nrow=247,ncol=5)
fpiarho[,1] <- c(1763:2009)
for (i in 1763:2009){
 fpiarho[i-1762,2] <- mean(fpidrho[yrm==i])
 fpiarho[i-1762,3:5] <- quantile(fpidrho[yrm==i],c(0.5,0.75,0.9))
}

#### write daily index
write.table(cbind(yrm,mon[selseas],day[selseas],fpidrho),"fpi_daily_Rhone_NDJFMA.txt",sep=" ",row.names=F,col.names=c("yr","mon","day","fpi_daily_PonteTresa"))

### write annual index
write.table(fpiarho,"fpi_annual_Rhone_NDJFMA.txt",sep=" ",row.names=F,col.names=c("yr","mean","50th_percentile","75th_percentile","90th_percentile"))
