diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..8c72037204e6310847da6daed5bc7fb7b336f59d
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,22 @@
+stages:  # List of stages for jobs, and their order of execution
+  - test
+
+unit-test-seasonal:  # This job runs in the test stage.
+  stage: test 
+  script:
+    - echo "Loading modules..."
+    - module load R/4.1.2-foss-2015a-bare
+    - module load CDO/1.9.8-foss-2015a
+    - module list
+    - echo "Running seasonal unit tests..."
+    - Rscript ./tests/test_seasonal.R
+
+unit-test-decadal:  # This job runs in the test stage.
+  stage: test
+  script:
+    - echo "Loading modules..."
+    - module load R/4.1.2-foss-2015a-bare
+    - module load CDO/1.9.8-foss-2015a
+    - module list
+    - echo "Running decadal unit tests..."
+    - Rscript ./tests/test_decadal.R
diff --git a/modules/Loading/Loading_decadal.R b/modules/Loading/Loading_decadal.R
index fe615594d41f61c80595524cc88c5775ce1e3302..f94d134315bfc8f144a738d2d8a645a6ebc14f12 100644
--- a/modules/Loading/Loading_decadal.R
+++ b/modules/Loading/Loading_decadal.R
@@ -2,11 +2,11 @@
 # 1. archive.yml
 # 2. recipe.yml
 # 3. Load_decadal.R (V)
-setwd('/esarchive/scratch/aho/git/auto-s2s/')
+#setwd('/esarchive/scratch/aho/git/auto-s2s/')
 
 ## TODO: remove paths to personal scratchs 
 source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R")
-source("/esarchive/scratch/vagudets/repos/cstools/R/s2dv_cube.R")
+#source("/esarchive/scratch/vagudets/repos/cstools/R/s2dv_cube.R")
 # Load required libraries/funs
 source("modules/Loading/get_daily_time_ind.R")
 source("modules/Loading/check_latlon.R")
diff --git a/tests/recipes/recipe-decadal_monthly_1.yml b/tests/recipes/recipe-decadal_monthly_1.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5acd642d1e6b1b6ad8512bf4d394671de2b5f8bb
--- /dev/null
+++ b/tests/recipes/recipe-decadal_monthly_1.yml
@@ -0,0 +1,46 @@
+Description:
+  Author: An-Chi Ho
+  '': split version
+Analysis:
+  Horizon: Decadal
+  Variables:
+    name: tas
+    freq: monthly_mean
+  Datasets:
+    System: 
+      name: EC-Earth3-i4 #CanESM5 
+      member: r1i4p1f1,r2i4p1f1
+    Multimodel: no
+    Reference:
+      name: ERA5 #JRA-55
+  Time:
+    sdate:
+    fcst:  2021
+#      fcst_sday: '1101'
+    hcst_start: 1991 #'1993'
+    hcst_end: 1994 #'2016'
+#    season: 'Annual'
+    leadtimemin: 0
+    leadtimemax: 2
+  Region:
+    latmin: 17 #-90
+    latmax: 20 #90
+    lonmin: 12
+    lonmax: 15 #359.9
+  Regrid:
+    method: bilinear
+    type: to_system #to_reference
+  Workflow:
+    Calibration:
+      method: bias
+    Skill:
+      metric: RPSS
+    Indicators:
+      index: FALSE
+  Output_format: S2S4E
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/aho/git/auto-s2s/out-logs/
+  code_dir: /esarchive/scratch/aho/git/auto-s2s/
+
diff --git a/tests/recipes/recipe-seasonal_monthly_1.yml b/tests/recipes/recipe-seasonal_monthly_1.yml
new file mode 100644
index 0000000000000000000000000000000000000000..3631c783ab08f5b895e62e0f2e70e82dee34aa8d
--- /dev/null
+++ b/tests/recipes/recipe-seasonal_monthly_1.yml
@@ -0,0 +1,43 @@
+Description:
+  Author: V. Agudetse
+
+Analysis:
+  Horizon: Seasonal
+  Variables:
+    name: tas
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: system7c3s
+    Multimodel: False
+    Reference:
+      name: era5
+  Time:
+    sdate:
+      fcst_syear: '2020'
+      fcst_sday: '1101'
+    hcst_start: '1993'
+    hcst_end: '1996'
+    leadtimemin: 0
+    leadtimemax: 2
+  Region:
+    latmin: 17
+    latmax: 20
+    lonmin: 12
+    lonmax: 15
+  Regrid:
+    method: bilinear
+    type: to_system
+  Workflow:
+    Calibration:
+      method: mse_min
+    Skill:
+      metric: RPSS
+    Indicators:
+      index: no
+  Output_format: S2S4E
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/
+  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
diff --git a/tests/test_decadal.R b/tests/test_decadal.R
new file mode 100644
index 0000000000000000000000000000000000000000..84a23b6275c68be35bb1d53820f906530d3f95ed
--- /dev/null
+++ b/tests/test_decadal.R
@@ -0,0 +1,8 @@
+library(testthat)
+
+path_testthat <- file.path('./tests/testthat/') 
+files_testthat <- list.files('./tests/testthat/', pattern = 'decadal')
+
+for (i_file in 1:length(files_testthat)) {
+  source(paste0('./tests/testthat/', files_testthat[i_file]))
+}
diff --git a/tests/test_seasonal.R b/tests/test_seasonal.R
new file mode 100644
index 0000000000000000000000000000000000000000..4718e3d4d1b224ba390c3a7767519ca42492c4d2
--- /dev/null
+++ b/tests/test_seasonal.R
@@ -0,0 +1,9 @@
+library(testthat)
+
+path_testthat <- file.path('./tests/testthat/')
+files_testthat <- list.files('./tests/testthat/', pattern = 'seasonal')
+
+for (i_file in 1:length(files_testthat)) {
+  source(paste0('./tests/testthat/', files_testthat[i_file]))
+}
+
diff --git a/tests/testthat/test-decadal_monthly.R b/tests/testthat/test-decadal_monthly.R
new file mode 100644
index 0000000000000000000000000000000000000000..19560921e95b6265e522a7f5351a63dfbb52bb77
--- /dev/null
+++ b/tests/testthat/test-decadal_monthly.R
@@ -0,0 +1,197 @@
+context("Decadal monthly data - 1")
+
+###########################################
+
+source("modules/Loading/Loading_decadal.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Save_data.R")
+
+recipe_file <- "tests/recipes/recipe-decadal_monthly_1.yml"
+
+# Load datasets
+suppressWarnings({invisible(capture.output(
+data <- load_datasets(recipe_file)
+))})
+
+recipe <- read_yaml(recipe_file)
+# Calibrate datasets
+suppressWarnings({invisible(capture.output(
+  calibrated_data <- calibrate_datasets(data, recipe)
+))})
+
+# Compute skill metrics
+suppressWarnings({invisible(capture.output(
+skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs,
+                                       recipe, na.rm = T, ncores = 4)
+))})
+
+#======================================
+
+test_that("1. Loading", {
+
+expect_equal(
+is.list(data),
+TRUE
+)
+expect_equal(
+names(data),
+c("hcst", "fcst", "obs")
+)
+expect_equal(
+class(data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(data$fcst),
+"s2dv_cube"
+)
+expect_equal(
+class(data$obs),
+"s2dv_cube"
+)
+expect_equal(
+names(data$hcst),
+c("data", "Variable", "Datasets", "Dates", "when", "source_files", "load_parameters")
+)
+expect_equal(
+names(data$hcst),
+names(data$fcst)
+)
+expect_equal(
+names(data$hcst),
+names(data$obs)
+)
+expect_equal(
+dim(data$hcst$data),
+c(dat = 1, var = 1, ensemble = 2, sday = 1, sweek = 1, syear = 4, time = 3, latitude = 5, longitude = 4)
+)
+expect_equal(
+dim(data$fcst$data),
+c(dat = 1, var = 1, ensemble = 2, sday = 1, sweek = 1, syear = 1, time = 3, latitude = 5, longitude = 4)
+)
+expect_equal(
+dim(data$hcst$Dates$start),
+c(sday = 1, sweek = 1, syear = 4, time = 3)
+)
+expect_equal(
+as.vector(drop(data$hcst$data)[,1:2,1,2,3]),
+c(291.3831, 291.6227, 292.3012, 290.9779),
+tolerance = 0.0001
+)
+expect_equal(
+mean(data$hcst$data),
+287.3804,
+tolerance = 0.0001
+)
+expect_equal(
+range(data$hcst$data),
+c(281.7395, 294.2467),
+tolerance = 0.0001
+)
+expect_equal(
+(data$hcst$Dates$start)[1],
+as.POSIXct("1991-11-16", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$Dates$start)[2],
+as.POSIXct("1992-11-16", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$Dates$start)[5],
+as.POSIXct("1991-12-16 12:00:00", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$Dates$start)[10],
+as.POSIXct("1993-01-16 12:00:00", tz = 'UTC')
+)
+
+})
+
+#======================================
+test_that("2. Calibration", {
+
+expect_equal(
+is.list(calibrated_data),
+TRUE
+)
+expect_equal(
+names(calibrated_data),
+c("hcst", "fcst")
+)
+expect_equal(
+class(calibrated_data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(calibrated_data$fcst),
+"s2dv_cube"
+)
+expect_equal(
+dim(calibrated_data$hcst$data),
+c(dat = 1, var = 1, ensemble = 2, sday = 1, sweek = 1, syear = 4, time = 3, latitude = 5, longitude = 4)
+)
+expect_equal(
+dim(calibrated_data$fcst$data),
+c(dat = 1, var = 1, ensemble = 2, sday = 1, sweek = 1, syear = 1, time = 3, latitude = 5, longitude = 4)
+)
+expect_equal(
+mean(calibrated_data$fcst$data),
+291.8375,
+tolerance = 0.0001
+)
+expect_equal(
+mean(calibrated_data$hcst$data),
+289.6679,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(calibrated_data$hcst$data)[1, , 2, 3, 4]),
+c(286.3895, 286.6408, 290.6652, 288.3759),
+tolerance = 0.0001
+)
+expect_equal(
+range(calibrated_data$fcst$data),
+c(287.2173, 297.4578),
+tolerance = 0.0001
+)
+
+})
+
+
+#======================================
+test_that("3. Metrics", {
+
+expect_equal(
+is.list(skill_metrics),
+TRUE
+)
+expect_equal(
+names(skill_metrics),
+c("rpss", "rpss_significance")
+)
+expect_equal(
+class(skill_metrics$rpss[[1]]),
+"array"
+)
+expect_equal(
+dim(skill_metrics$rpss[[1]]),
+c(dat = 1, var = 1, sday = 1, sweek = 1, time = 3, latitude = 5, longitude = 4)
+)
+expect_equal(
+dim(skill_metrics$rpss_significance[[1]]),
+dim(skill_metrics$rpss[[1]])
+)
+expect_equal(
+as.vector(drop(skill_metrics$rpss[[1]])[, 2, 3]),
+c(-0.2857143, -1.2500000, -1.8928571),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(skill_metrics$rpss_significance[[1]])[, 2, 3]),
+rep(FALSE, 3)
+)
+
+})
+
+
diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R
new file mode 100644
index 0000000000000000000000000000000000000000..4bdf833d44ff33ba185508e23b2e49dad5447ea5
--- /dev/null
+++ b/tests/testthat/test-seasonal_monthly.R
@@ -0,0 +1,190 @@
+context("Seasonal monthly data")
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Save_data.R")
+
+recipe_file <- "tests/recipes/recipe-seasonal_monthly_1.yml"
+
+# Load datasets
+suppressWarnings({invisible(capture.output(
+data <- load_datasets(recipe_file)
+))})
+
+recipe <- read_yaml(recipe_file)
+
+suppressWarnings({invisible(capture.output(
+calibrated_data <- calibrate_datasets(data, recipe)
+))})
+
+# Compute skill metrics
+suppressWarnings({invisible(capture.output(
+skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs,
+                                       recipe, na.rm = T, ncores = 4)
+))})
+
+test_that("1. Loading", {
+
+expect_equal(
+is.list(data),
+TRUE
+)
+expect_equal(
+names(data),
+c("hcst", "fcst", "obs")
+)
+expect_equal(
+class(data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(data$fcst),
+"s2dv_cube"
+)
+expect_equal(
+class(data$obs),
+"s2dv_cube"
+)
+expect_equal(
+names(data$hcst),
+c("data", "lon", "lat", "Variable", "Datasets", "Dates", "when", "source_files", "load_parameters")
+)
+expect_equal(
+names(data$hcst),
+names(data$fcst)
+)
+expect_equal(
+names(data$hcst),
+names(data$obs)
+)
+expect_equal(
+dim(data$hcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 4, time = 3, latitude = 3, longitude = 3, ensemble = 25)
+)
+expect_equal(
+dim(data$fcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 1, time = 3, latitude = 3, longitude = 3, ensemble = 51)
+)
+expect_equal(
+dim(data$hcst$Dates$start),
+c(sday = 1, sweek = 1, syear = 4, time = 3)
+)
+expect_equal(
+as.vector(drop(data$hcst$data)[1:2,1:2,1,2,3]),
+c(293.9651, 295.9690, 290.6771, 290.7957),
+tolerance = 0.0001
+)
+expect_equal(
+mean(data$hcst$data),
+290.8758,
+tolerance = 0.0001
+)
+expect_equal(
+range(data$hcst$data),
+c(284.7413, 299.6219),
+tolerance = 0.0001
+)
+expect_equal(
+(data$hcst$Dates$start)[1],
+as.POSIXct("1993-12-01", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$Dates$start)[2],
+as.POSIXct("1994-12-01", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$Dates$start)[5],
+as.POSIXct("1994-01-01", tz = 'UTC')
+)
+expect_equal(
+(data$obs$Dates$start)[10],
+as.POSIXct("1995-02-14", tz = 'UTC')
+)
+
+})
+
+test_that("2. Calibration", {
+
+expect_equal(
+is.list(calibrated_data),
+TRUE
+)
+expect_equal(
+names(calibrated_data),
+c("hcst", "fcst")
+)
+expect_equal(
+class(calibrated_data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(calibrated_data$fcst),
+"s2dv_cube"
+)
+expect_equal(
+dim(calibrated_data$hcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 4, time = 3, latitude = 3, longitude = 3, ensemble = 25)
+)
+expect_equal(
+dim(calibrated_data$fcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 1, time = 3, latitude = 3, longitude = 3, ensemble = 51)
+)
+expect_equal(
+mean(calibrated_data$fcst$data),
+291.1218,
+tolerance = 0.0001
+)
+expect_equal(
+mean(calibrated_data$hcst$data),
+289.8596,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(calibrated_data$hcst$data)[1, , 2, 3, 4]),
+c(287.7982, 287.0422, 290.4297),
+tolerance = 0.0001
+)
+expect_equal(
+range(calibrated_data$fcst$data),
+c(283.5374, 306.2353),
+tolerance = 0.0001
+)
+
+})
+
+
+#======================================
+test_that("3. Metrics", {
+
+expect_equal(
+is.list(skill_metrics),
+TRUE
+)
+expect_equal(
+names(skill_metrics),
+c("rpss", "rpss_significance")
+)
+expect_equal(
+class(skill_metrics$rpss[[1]]),
+"array"
+)
+expect_equal(
+dim(skill_metrics$rpss[[1]]),
+c(dat = 1, var = 1, sday = 1, sweek = 1, time = 3, latitude = 3, longitude = 3)
+)
+expect_equal(
+dim(skill_metrics$rpss_significance[[1]]),
+dim(skill_metrics$rpss[[1]])
+)
+expect_equal(
+as.vector(drop(skill_metrics$rpss[[1]])[, 2, 3]),
+c(-1.153829, -1.114743, -1.392457),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(skill_metrics$rpss_significance[[1]])[, 2, 3]),
+rep(FALSE, 3)
+)
+
+})