# This is a script to read and process LOCA modeled data
# Resolution: 1/8 x 1/8 degree, 1950-2099
#===================================================================================================================
library(ncdf4)
library(caTools)
library(abind)
#===================================================================================================================

#---Extract 1990-2005 (same modeled baseline for both RCPs) and 2006-2099 separately, then bind them

scenarios=c('rcp45','rcp85')
vars = c('pr','pet') 
models = c('canesm2','ccsm4','gfdl-cm3','giss-e2-r','hadgem2-es','miroc5')

dir = '/groups/anenberggrp/PA/dust/LOCA/projections/'
xlim=c(-120,-100); ylim=c(25,45) # Extract SW+Mex
years=c(1950:2099)
ts = seq(as.Date("1950/01/01"), as.Date("2099/12/01"), by = "month")
tmin = format(min(ts), format='%Y%m'); tmax = format(max(ts), '%Y%m')

for (n in 1:length(scenarios)){
  for (j in 1:length(vars)){
    for (k in 1:length(models)){
      print(paste(scenarios[n],' +  ',vars[j],' + ',models[k],sep=''))
     
      #--- 1950-2005
      files = (Sys.glob(paste(dir,'hindcast/',vars[j],'/',models[k],'/*.nc', sep='')))
      raw.dat = nc_open(files[1])
      lon = ncvar_get(raw.dat, "centroidLons")
      lat = ncvar_get(raw.dat, "centroidLats")
      ind1=which(lon>=xlim[1] & lon<=xlim[2])
      ind2=which(lat>=ylim[1] & lat<=ylim[2])
      sub.lon = lon[ind1]; sub.lat=lat[ind2]

      df = ncvar_get(raw.dat, start=c(ind2[1],ind1[1],1), count=c(length(ind2),length(ind1),12))
      # lat x lon x months
      nc_close(raw.dat)

      for (i in 2:length(files)){
      #print(i)
      fname=files[i]
      raw.dat = nc_open(fname)
      x = ncvar_get(raw.dat, start=c(ind2[1],ind1[1],1), count=c(length(ind2),length(ind1),12))
      df = abind(df, x, along=3)
      nc_close(raw.dat)
      }
      print(dim(df))

      #--- 2006-2099
      files=(Sys.glob(paste(dir,scenarios[n],'/',vars[j],'/',models[k],'/*.nc', sep=''))); i=1
      for (i in 1:length(files)){
      #print(i)
      fname=files[i]
      raw.dat = nc_open(fname)
      x = ncvar_get(raw.dat, start=c(ind2[1],ind1[1],1), count=c(length(ind2),length(ind1),12))
      df = abind(df, x, along=3)
      nc_close(raw.dat) 
      }
      print(dim(df))
  
  new.lat = rev(sub.lat)
  df = aperm(df, c(2,1,3)) # change to lon x lat x months
  df = df[,dim(df)[2]:1,] # reverse lat to be from S to N

  if (j==2){ # need to convert daily to monthly
    m1=c(1,3,5,7,8,10,12); m2=c(4,6,9,11); m3=2
    temp = array(df, dim=c(length(sub.lon), length(new.lat), 12, dim(df)[3]/12))
    temp[,,m1,]=temp[,,m1,]*31; temp[,,m2,]=temp[,,m2,]*30; temp[,,m3,]=temp[,,m3,]*28
    df=array(temp, dim=dim(df))  
  }

  ss = list('data'=df, 'lon'=sub.lon, 'lat'=new.lat, 'time'=ts)
  print(dim(ss$data))

    #fillvalue = 1e+20 # missing value
    #x[x == fillvalue] = NA # set missing values to zero 

  fout = paste('./data/',scenarios[n],'/',vars[j],'/',vars[j],'_',models[k],'_',tmin,'-',tmax,'.Rdata', sep='') # filename
  save(ss, file=fout)
  print(fout)
  rm(df); rm(ss);
    }
  }
}

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