# This is a script to calculate changes in SPEI02 for 6 LOCA models x 2 RCPs
# 2006-2099 seasonal mean timeseries (relative to 1986-2005)
#===================================================================================================================
library(reshape2)
source('./my.functions.R')
#===================================================================================================================

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

  for (i in 1:length(models)){
  print(models[i])
    load(paste('./data/',scenarios[j],'/spei/spei02_',models[i],'_195001-209912.Rdata',sep='')) # res
    lon=res$lon; lat=res$lat; time=res$time
    
    for (n in 1:length(seas)){
      if(n==1){xlim=c(-117,-105); ylim=c(30,38)}
      if(n==2){xlim=c(-115,-103); ylim=c(28,38)}
      if(n==3){xlim=c(-115,-103); ylim=c(29,40)}
      if(n==4){xlim=c(-113,-103); ylim=c(30,39)}   
     
      ilon=which(lon >= xlim[1] & lon <= xlim[2]); nlon=lon[ilon]
      ilat=which(lat >= ylim[1] & lat <= ylim[2]); nlat=lat[ilat]
      itime=which(time >= "1985-12-01" & time <= "2099-11-01"); ntime=time[itime]
      
      temp=res$data[ilon,ilat,itime] # lon x lat x months
      df = array(temp, dim=c(length(ilon),length(ilat),12,dim(temp)[3]/12))  
      df2 = array(df, dim=c(length(ilon),length(ilat),3,4,dim(temp)[3]/12))
      df3 = apply(df2[,,,n,], c(1,2,4), FUN=mean, na.rm=TRUE) # seasmean
      df4 = apply(df3, 3, FUN=function(xx) calc_awa(xx, lon=nlon, lat=nlat)) # regmean      

      #---Calculate diff
      tnew=c(1986:2099); 
      t1 = which(tnew %in% c(1986:2005))
      t2 = which(tnew %in% c(2006:2099))

      vp=df4[t1]; vf=df4[t2]
      mp = mean(vp, na.rm=TRUE)
      diff[,n,i] = vf-mp
  }}  
  ss = list('diff'=diff, 'years'=years, 'seas'=seas, 'models'=models)
  save(ss, file=paste('./data/fpdiff/diff_spei02_',scenarios[j],'.Rdata',sep=''))
}

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