# This is a script to calculate changes in FD using slr results for 6 LOCA models x 2 RCPs
# 2006-2099 seasonal mean timeseries (relative to 1986-2005)
#===================================================================================================================
library(reshape2)
#===================================================================================================================

scenarios = c('rcp45','rcp85')
models = c('canesm2','ccsm4','gfdl-cm3','giss-e2-r','hadgem2-es','miroc5')
seas = c('DJF','MAM','JJA','SON')
years=c(2006:2099)
  
for (j in 1:2){ # RCPs 
  print(scenarios[j])
  diff = array(NA, dim=c(length(years), length(seas), length(models))) # year x seas x models

  #--- Load spei diff
  load(paste('./data/fpdiff/FD/spei02/diff_spei02_',scenarios[j],'.Rdata',sep='')) # ss
  spei = ss$diff # years x seas x models
  models = ss$models

  for (n in 1:length(seas)){

  #--- Load SLR results
    load(paste('../historical/data/slr/FD/spei02/slr_2000-2013-',seas[n],'.RData',sep=''))
    slr = ss$SLR
    diff[,n,] = slr$coefficients[1] + slr$coefficients[2]*spei[,n,]
  }

  ss = list('diff'=diff, 'years'=years, 'seas'=seas, 'models'=models)
  save(ss, file=paste('./data/fpdiff/FD/spei02/diff_FD_',scenarios[j],'.Rdata',sep=''))
}

#===================================================================================================================

#--- Calculate annual mean and prepare .csv output

years=c(2006:2099)
scenarios = c('rcp45','rcp85')
models = c('canesm2','ccsm4','gfdl-cm3','giss-e2-r','hadgem2-es','miroc5')

data = array(NA, dim=c(length(years),12))

for (j in 1:2){ # RCPs
  load(paste('./data/fpdiff/FD/spei02/diff_FD_',scenarios[j],'.Rdata',sep=''))
  df = ss$diff
  df2 = apply(df, c(1,3), FUN=mean, na.rm=TRUE) # annual mean

  if(j==1){data[,1:6]=df2}
  if(j==2){data[,7:12]=df2}
}

cnames=vector()
for (i in 1:6){cnames[i] = paste(scenarios[1],'.',models[i],sep='')}
for (i in 7:12){cnames[i] = paste(scenarios[2],'.',models[i-6],sep='')}

colnames(data)=cnames
rownames(data)=years

save(data, file='./data/fpdiff/FD/spei02/diff_FD_annual.Rdata')
write.csv(data, file='../outputs/diff-FD_2006-2099-annual.csv')

#===================================================================================================================

#--- Calculate 20-year averages from annual means
load('./data/fpdiff/FD/spei02/diff_FD_annual.Rdata')

years=c(2006:2099)
i1 = which(years %in% c(2020:2039))
i2 = which(years %in% c(2040:2059))
i3 = which(years %in% c(2060:2079))
i4 = which(years %in% c(2080:2099))

ss = array(NA, dim=c(4,12))
ss[1,] = apply(data[i1,], 2, FUN=mean, na.rm=TRUE)
ss[2,] = apply(data[i2,], 2, FUN=mean, na.rm=TRUE)
ss[3,] = apply(data[i3,], 2, FUN=mean, na.rm=TRUE)
ss[4,] = apply(data[i4,], 2, FUN=mean, na.rm=TRUE)
rownames(ss)=c('2030','2050','2070','2090')
colnames(ss)=colnames(data)

save(ss, file='./data/fpdiff/FD/spei02/diff_FD_20y.Rdata')
write.csv(ss, file='../outputs/diff-FD_20y.csv')

#===================================================================================================================