# -----------------------------------------------------
# Loading tas and tos for EC-Earth decadal simulations:
# Authors: Carlos Delgado and Núria Pérez-Zanón
# ------------------------------------------------------

# Three ways to load the same data are provided:
# 1) single Start call providing two paths and two variable names
# 2) single Start call providing one path and two variabe names
# 3) two Start call one for each path and variable


# Case 1) returns dimensions 'dataset' and 'var' with length 2 , but only the positions of the diagonal are filled: 
# tas is stored in {dataset = 1, var = 1}
# tos is stored in {dataset = 2, var = 2}
# NOTE!!! check {datastet = 1, var = 2} because an issue in ESMValTool:https://earth.bsc.es/gitlab/es/auto-ecearth3/issues/1258

library(startR)

paths = list(list(path = '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Amon/$var$/gr/v20190713/$var$_Amon_EC-Earth3_dcppA-hindcast_s$sdate$-$member$_gr_$fyear$.nc'),
             list(path = '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Omon/$var$/gr/v20190713/$var$_Omon_EC-Earth3_dcppA-hindcast_s$sdate$-$member$_gr_$fyear$.nc'))
data1 <- Start(dataset = paths,
               var = c('tas','tos'),
               sdate = paste0(1960:1962),
               fmonth = 1,
               lat = values(list(0, 10)),
               lon = values(list(0, 10)),
               fyear = 'all',
               member = indices(1),
               fyear_depends = 'sdate',
               fmonth_across = 'fyear',
               merge_across_dims = TRUE,
               synonims = list(fmonth = c('fmonth','time'),
                               lon = c('lon', 'longitude'),
                               lat = c('lat', 'latitude')),
               return_vars = list(lat = 'dataset', lon = 'dataset'),
               num_procs = 1, retrieve = TRUE)

dim(data1)
#dataset     var   sdate  fmonth     lat     lon  member 
#      2       2       3      12      14      15       1 
# Check empty and filled dimensions:
sum(is.na(data1[1,1,,,,,])) == (3*12*14*15)
#[1] FALSE
sum(is.na(data1[1,2,,,,,])) == (3*12*14*15)
#[1] TRUE 
sum(is.na(data1[2,2,,,,,])) == (3*12*14*15)
#[1] FALSE
sum(is.na(data1[2,1,,,,,])) == (3*12*14*15)
#[1] TRUE

lat1 <- as.vector(attributes(data1)$Variables$dat1$lat)
lon1 <- as.vector(attributes(data1)$Variables$dat1$lon)

# ---------------------------------------------------------------

# Case 2) using a single path, {dataset = 1, var = 2, type = 1}
# 'type' dimension is necessary to distinguish between 'Amon' and 'Omon'. 

library(startR)

path = '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/$type$/$var$/gr/v20190713/$var$_$type$_EC-Earth3_dcppA-hindcast_s$sdate$-$member$_gr_$fyear$.nc'
data2 <- Start(dataset = path,
               var = c('tas', 'tos'),
               type = 'all',
               type_depends = 'var',
               sdate = paste0(1960:1962),
               fmonth = 1,
               lat = values(list(0, 10)),
               lon = values(list(0, 10)),
               fyear = indices(1),
               member = indices(1),
               fyear_depends = 'sdate',
               fmonth_across = 'fyear',
               merge_across_dims = TRUE,
               synonims = list(fmonth = c('fmonth','time'),
                               lon = c('lon', 'longitude'),
                               lat = c('lat', 'latitude')),
               return_vars = list(lat = 'dataset', lon = 'dataset'),
               num_procs = 1, retrieve = TRUE)

dim(data2)
# dataset     var    type   sdate  fmonth     lat     lon  member 
#      1       2       1       3       1      14      15       1 


# ---------------------------------------------------------------

# Case 3) Two different Start calls can save data_tas and data_tos both with {dataset = 1 and var = 1} dimensions and avoiding extra dimensions like 'type'.

path = '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Amon/$var$/gr/v20190713/$var$_Amon_EC-Earth3_dcppA-hindcast_s$sdate$-$member$_gr_$fyear$.nc'
data_tas <- Start(dataset = path,
               var = 'tas',
               sdate = paste0(1960:1962),
               fmonth = 1,
               lat = values(list(0, 10)),
               lon = values(list(0, 10)),
               fyear = indices(1),
               member = indices(1),
               fyear_depends = 'sdate',
               fmonth_across = 'fyear',
               merge_across_dims = TRUE,
               synonims = list(fmonth = c('fmonth','time'),
                               lon = c('lon', 'longitude'),
                               lat = c('lat', 'latitude')),
               return_vars = list(lat = 'dataset', lon = 'dataset'),
               num_procs = 1, retrieve = TRUE)

dim(data_tas)
#dataset     var   sdate  fmonth     lat     lon  member 
#      1       1       3       1      14      15       1 

path = '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/$member$/Omon/$var$/gr/v20190713/$var$_Omon_EC-Earth3_dcppA-hindcast_s$sdate$-$member$_gr_$fyear$.nc'
data_tos <- Start(dataset = path,
               var = 'tos',
               sdate = paste0(1960:1962),
               fmonth = 1,
               lat = values(list(0, 10)),
               lon = values(list(0, 10)),
               fyear = indices(1),
               member = indices(1),
               fyear_depends = 'sdate',
               fmonth_across = 'fyear',
               merge_across_dims = TRUE,
               synonims = list(fmonth = c('fmonth','time'),
                               lon = c('lon', 'longitude'),
                               lat = c('lat', 'latitude')),
               return_vars = list(lat = 'dataset', lon = 'dataset'),
               num_procs = 1, retrieve = TRUE)

dim(data_tos)
#dataset     var   sdate  fmonth     lat     lon  member 
#      1       1       3       1      14      15       1 

# ---------------------------------------------------------------------

# Comparison cases 1) to 3):
#----------------------------------------------------------------------

all(data1[1, 1, , , , , ] == data_tas[1, 1, , , , , ])
all((data1[2, 2, , , , , ]) == data_tos[1, 1, , , , , ], na.rm = TRUE)
all(data2[1, 1, 1, , , , , ] == data_tas[1, 1, , , , ,])
all((data2[1, 2, 1, , , , , ]) == data_tos[1, 1, , , , , ], na.rm = TRUE)