#### Aggregating data to daily frequency ####

# Here, we provide code for the CDO command line tool to aggregate higher frequency data to daily means
# Both ERA5 temperature and precipitation are converted from hourly to daily
# For CMIP5, FGOALS-g3 tas is converted from 3-hourly to daily
# For CMIP6, AWI-CM-1-1-MR pr is converted from 3-hourly to daily
# After running this script, hourly and 3-hourly data can be deleted

library(tidyverse)
library(ncdf4)

years = paste(1979:2005)
months = c('01','02','03','04','05','06','07','08','09','10','11','12')

data_dir = '/mnt/r/historical/'

#### ERA5 Temperature ####

for(y in years){
  for(m in months){
    era5 <- nc_open(paste0(data_dir,'era5/tas/hourly/era5_tas_hourly_',y,'_',m,'.nc'))
    hourly <- ncvar_get(era5, "t2m") - 273.15
    
    #Calculate hours and days
    n_hours = dim(hourly)[3]
    n_days = n_hours/24
    
    #iterate and record each daily mean
    daily = array(NA,dim = c(dim(hourly)[1:2],n_days))
    for(d in 1:n_days){
      daily[,,d] = apply(hourly[,,(d-1)*24 + 1:24],1:2,mean)
    }
    
    rm(hourly)
    gc()
    
    #Reformat to day, lat, long
    daily = aperm(daily,3:1)
    
    #Save daily data
    saveRDS(daily,paste0(data_dir,'era5/tas/daily/era5_tas_daily_',y,'_',m,'.rds'))
    rm(daily)
    gc()
  }
}

#### ERA5 Precipitation ####

for(y in years){
  for(m in months){
    era5 <- nc_open(paste0(data_dir,'era5/pr/hourly/era5_pr_hourly_',y,'_',m,'.nc'))
    hourly <- ncvar_get(era5, "tp")
    
    #Calculate hours and days
    n_hours = dim(hourly)[3]
    n_days = n_hours/24
    
    #iterate and record each daily mean
    daily = array(NA,dim = c(dim(hourly)[1:2],n_days))
    for(d in 1:n_days){
      daily[,,d] = apply(hourly[,,(d-1)*24 + 1:24],1:2,sum)
    }
    
    rm(hourly)
    gc()
    
    #Reformat to day, lat, long
    daily = aperm(daily,3:1)
    
    #Save daily data
    saveRDS(daily,paste0(data_dir,'era5/pr/daily/era5_pr_daily_',y,'_',m,'.rds'))
    rm(daily)
    gc()
  }
}

#### Convert FGOALS-g2 3-hourly to daily ####

setwd(paste0(data_dir,'3hr/FGOALS-g2/'))
fgoals_files = sort(list.files())[1:27] #sort by year and subset files we need

#use cdo to create our upsampling grid and save the template for later
for(i in 1:length(1979:2005)){
  in_file = fgoals_files[i]
  out_file = gsub('3hr','day',in_file)
  out_file = paste0(substr(out_file,1,44),'-',as.numeric(substr(out_file,50,53))-1,'1231','.nc')
  system(paste0('cdo shifttime,1sec -daymean -shifttime,-1sec ',data_dir,'3hr/FGOALS-g2/',
                in_file, ' ',
                data_dir,'cmip5tas/',
                out_file))
}

#### Convert AWI-CM-1-1-MR 3-hourly to daily ####

setwd(paste0(data_dir,'3hr/AWI-CM-1-1-MR/'))
awi_files = sort(list.files())[1:27] #sort by year and subset files we need

#use cdo to create our upsampling grid and save the template for later
for(i in 1:length(1979:2005)){
  in_file = awi_files[i]
  out_file = gsub('3hr','day',in_file)
  out_file = paste0(substr(out_file,1,52),substr(out_file,57,65),'.nc')
  system(paste0('cdo daymean ',data_dir,'3hr/AWI-CM-1-1-MR/',
         in_file, ' ',
         data_dir,'cmip6pr/',
         out_file))
}
