diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index fe0e593fa26ac3cc245f147f4f57ce743aa8f863..3eb90d50935404fa29137db76f00eae65305193a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ build:
   stage: build
   script:
     - module load R/3.6.1-foss-2015a-bare
-#    - module load CDO
+    - module load CDO/1.9.8-foss-2015a
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual --run-donttest s2dv_*.tar.gz
     - R -e 'covr::package_coverage()'
diff --git a/R/CDORemap.R b/R/CDORemap.R
index fc25b527de4036257c462b28a0d5eb27eac89154..19e0197d3b9403c57116239445371b5087c63ba1 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -251,22 +251,27 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
       array_dims <- c(length(lats), length(lons))
       new_lon_dim_name <- 'lon'
       new_lat_dim_name <- 'lat'
-    } else {
-      array_dims <- dim(lons)
-      new_lon_dim_name <- 'i'
-      new_lat_dim_name <- 'j'
-    }
-    if (!is.null(names(dim(lons)))) {
-      if (any(known_lon_names %in% names(dim(lons)))) {
-        new_lon_dim_name <- known_lon_names[which(known_lon_names %in% names(dim(lons)))[1]]
+
+      if (!is.null(names(dim(lons)))) {
+        if (any(known_lon_names %in% names(dim(lons)))) {
+          new_lon_dim_name <- known_lon_names[which(known_lon_names %in% names(dim(lons)))[1]]
+        }
       }
-    }
-    if (!is.null(names(dim(lats)))) {
-      if (any(known_lat_names %in% names(dim(lats)))) {
-        new_lat_dim_name <- known_lat_names[which(known_lat_names %in% names(dim(lats)))[1]]
+      if (!is.null(names(dim(lats)))) {
+        if (any(known_lat_names %in% names(dim(lats)))) {
+          new_lat_dim_name <- known_lat_names[which(known_lat_names %in% names(dim(lats)))[1]]
+        }
+      }
+      names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name)
+
+    } else {  # irregular
+      array_dims <- dim(lons)
+      if (is.null(names(array_dims))) {
+        new_lon_dim_name <- 'i'
+        new_lat_dim_name <- 'j'
       }
     }
-    names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name)
+
     data_array <- array(as.numeric(NA), array_dims)
   }
   if (!(is.logical(data_array) || is.numeric(data_array)) || !is.array(data_array)) {
@@ -890,7 +895,16 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
       } else {
         new_dims <- dim(data_array)
         new_dims[c(lon_dim, lat_dim)] <- c(found_lon_dim_size, found_lat_dim_size)
-
+        if (is_irregular) {
+          lon_pos <- which(names(new_dims) == lon_dim)
+          lat_pos <- which(names(new_dims) == lat_dim)
+          if (lon_pos > lat_pos) {
+            new_pos <- 1:length(new_dims)
+            new_pos[lon_pos] <- lat_pos
+            new_pos[lat_pos] <- lon_pos
+            new_dims <- new_dims[new_pos]
+          }
+        }
         result_array <- ncvar_get(ncdf_remapped, 'var', collapse_degen = FALSE)
         dim(result_array) <- new_dims
       }
diff --git a/tests/testthat/test-CDORemap.R b/tests/testthat/test-CDORemap.R
new file mode 100644
index 0000000000000000000000000000000000000000..4d70ca7f39eed917fe0f3ade387116a4e759d803
--- /dev/null
+++ b/tests/testthat/test-CDORemap.R
@@ -0,0 +1,312 @@
+context("s2dv::CDORemap tests")
+
+# data1: regular grid
+data1 <- array(1:360*181*2, dim = c(lon = 360, lat = 181, time = 2))
+lons1 <- seq(0, 359)
+lats1 <- seq(-90, 90)
+data1_1 <- s2dv:::.aperm2(data1, c(3, 1, 2))
+data1_2 <- s2dv::InsertDim(data1, 1, 1, name = 'var')
+data1_3 <- s2dv::InsertDim(data1_2, 5, 1, name = 'dat')
+
+# data2: irregular grid
+path2 <- '/esarchive/exp/ipsl-cm6a-lr/cmip6-dcppA-hindcast_i1p1/original_files/cmorfiles/DCPP/IPSL/IPSL-CM6A-LR/dcppA-hindcast/r1i1p1f1/SImon/siconc/gn/v20200101/$var$_SImon_IPSL-CM6A-LR_dcppA-hindcast_s2016-r1i1p1f1_gn_201701-202612.nc'
+library(startR)
+suppressWarnings(
+data2 <- Start(dat = path2,
+              var = 'siconc',
+              x = 'all', #indices(2:361),
+              y = 'all', #indices(2:331),
+              time = indices(1), #'all',
+              return_vars = list(#x = NULL, y = NULL, 
+                                 nav_lat = NULL, nav_lon = NULL),
+              retrieve = TRUE, silent = T)
+)
+lons2 <- attributes(data2)$Variables$common$nav_lon#[2:361, 2:331]
+lats2 <- attributes(data2)$Variables$common$nav_lat#[2:361, 2:331]
+
+data2_1 <- ClimProjDiags::Subset(data2, 1, 1, drop = 'selected')
+data2_2 <- ClimProjDiags::Subset(data2, 5, 1, drop = 'selected')
+data2_3 <- ClimProjDiags::Subset(data2, c(1, 2), list(1, 1), drop = 'selected')
+data2_4 <- ClimProjDiags::Subset(data2, c(1, 2, 5), list(1, 1, 1), drop = 'selected')
+data2_5 <- s2dv:::.aperm2(data2, c(3,4,1,2,5))
+data2_6 <- s2dv:::.aperm2(data2, c(3,4,5,2,1))
+data2_7 <- s2dv:::.aperm2(data2, c(4,3,1,2,5))
+data2_8 <- s2dv:::.aperm2(data2_4, c(2,1))
+
+# data3: irregular grid
+path3 <- '/esarchive/exp/ecearth/a1ua/cmorfiles/DCPP/EC-Earth-Consortium/EC-Earth3/dcppA-hindcast/r1i1p1f1/Omon/$var$/gn/v20190713/$var$_*_s$sdate$-$member$_gn_$aux$.nc'
+suppressWarnings(
+data3 <- Start(dataset = path3,
+               var = 'tos',
+               sdate = paste0(1960),
+               aux = indices(1),
+               aux_depends = 'sdate',
+               j = 'all',
+               i = 'all',
+               time = indices(1),
+               member = 'r1i1p1f1',
+               return_vars = list(j = NULL, i = NULL,
+                                  latitude = NULL, longitude = NULL),
+               retrieve = T, silent = T)
+)
+lons3 <- attributes(data3)$Variables$common$longitude
+lats3 <- attributes(data3)$Variables$common$latitude
+
+data3_1 <- drop(data3)
+data3_2 <- s2dv:::.aperm2(data3_1, c(2, 1))
+
+# data4: irregular grid
+path4 <- paste0('/esarchive/exp/CMIP6/dcppA-hindcast/cmcc-cm2-sr5/cmip6-dcppA-hindcast_i1p1/',
+               'DCPP/CMCC/CMCC-CM2-SR5/dcppA-hindcast/$member$/Omon/$var$/gn/v20210312/',
+               '$var$_*_s$sdate$-$member$_gn_$aux$.nc')
+suppressWarnings(
+data4 <- Start(dataset = path4,
+              var = 'tos',
+              sdate = '1960',
+              aux = 'all',
+              aux_depends = 'sdate',
+              j = 'all',
+              i = 'all',
+              time = indices(1),
+              member = 'r1i1p1f1',
+              return_vars = list(j = NULL, i = NULL,
+                                 latitude = NULL, longitude = NULL),
+              retrieve = T, silent = T)
+)
+lons4 <- attributes(data4)$Variables$common$longitude
+lats4 <- attributes(data4)$Variables$common$latitude
+
+data4_1 <- drop(data4)
+data4_2 <- ClimProjDiags::Subset(data4, c(1,2,3,4,8), list(1,1,1,1,1), drop = 'selected')
+data4_3 <- s2dv:::.aperm2(data4_2, c(1, 3, 2))
+
+##############################################
+
+test_that("1. Input checks", {
+
+expect_error(
+CDORemap(data_array = NULL, lons = 'lons', lats = 'lats'),
+"Expected numeric 'lons' and 'lats'."
+)
+expect_error(
+CDORemap(data_array = NULL, lons = c(NA, lons1), lats = lats1),
+"Found invalid values in 'lons'."
+)
+expect_error(
+CDORemap(data_array = NULL, lons = lons1, lats = c(NA, lats1)),
+"Found invalid values in 'lats'."
+)
+
+})
+
+test_that("2. regular regrid", {
+suppressWarnings(
+res <- CDORemap(data1, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data1_1, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data1_2, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data1_3, lons1, lats1, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+
+expect_equal(
+as.vector(res),
+as.vector(s2dv:::.aperm2(res1, c(2,3,1)))
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+
+})
+
+##################################################################
+
+test_that("3. dat2: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data2, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data2_1, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data2_2, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data2_3, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res4 <- CDORemap(data2_4, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res5 <- CDORemap(data2_5, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res6 <- CDORemap(data2_6, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res7 <- CDORemap(data2_7, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res8 <- CDORemap(data2_8, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res9 <- CDORemap(NULL, lons2, lats2, grid = 'r100x50', method = 'bil', crop = F)
+)
+
+expect_equal(
+drop(res)[95:100, 50],
+c(99.54009, 99.53733, 99.53453, 99.53169, 99.52792, 99.52304),
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+expect_equal(
+as.vector(res),
+as.vector(res4)
+)
+expect_equal(
+as.vector(res),
+as.vector(res5)
+)
+expect_equal(
+as.vector(res),
+as.vector(res6)
+)
+expect_equal(
+as.vector(res),
+as.vector(res7)
+)
+expect_equal(
+as.vector(res),
+as.vector(res8)
+)
+expect_equal(
+res9$data_array,
+NULL
+)
+expect_equal(
+as.vector(res9$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res9$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})
+
+#########################################################
+
+test_that("4. dat3: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data3, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data3_1, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data3_2, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res3 <- CDORemap(NULL, lons3, lats3, grid = 'r100x50', method = 'bil', crop = F)
+)
+
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+res1[1:5, 50],
+c(-1.588068, -1.586225, -1.584405, -1.582579, -1.580713),
+tolerance = 0.0001
+)
+expect_equal(
+res3$data_array,
+NULL
+)
+expect_equal(
+as.vector(res3$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res3$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})
+
+############################################################
+
+test_that("5. dat4: irregular regrid", {
+suppressWarnings(
+res <- CDORemap(data4, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res1 <- CDORemap(data4_1, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res2 <- CDORemap(data4_2, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+suppressWarnings(
+res3 <- CDORemap(data4_3, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)$data_array
+)
+# no data
+suppressWarnings(
+res4 <- CDORemap(NULL, lons4, lats4, grid = 'r100x50', method = 'bil', crop = F)
+)
+expect_equal(
+as.vector(res),
+as.vector(res1)
+)
+expect_equal(
+as.vector(res),
+as.vector(res2)
+)
+expect_equal(
+as.vector(res),
+as.vector(res3)
+)
+expect_equal(
+res1[1:5, 50],
+c(-1.795211, -1.795057, -1.794898, -1.794707, -1.794533),
+tolerance = 0.0001
+)
+expect_equal(
+res4$data_array,
+NULL
+)
+expect_equal(
+as.vector(res4$lons)[1:5],
+c(0.0, 3.6, 7.2, 10.8, 14.4)
+)
+expect_equal(
+as.vector(res4$lats)[1:5],
+c(-88.2, -84.6, -81.0, -77.4, -73.8)
+)
+
+})