diff --git a/NEWS.md b/NEWS.md
index 9979a52e4c114c1e8a2e1534f5b4f22b8bd9e27d..155781620c9b322f77e20506e4f246494dc00fbf 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,7 +8,8 @@
 - Add p-value by ANOVA in Trend().
 - Change MeanDims() na.rm default to FALSE to be in line with mean()
 - Remove unecessary parameter checks in Clim()
-
+- Bugfix for Clim() per-paired NA exclusion.
+ 
 # s2dv 0.0.1 (Release date: 2020-02-07)
 - The package is the advanced version of package 's2dverification', adopting the regime of package 'multiApply' for all the analytic functions. Most of the other functions for plotting and data retrieval in 's2dverification' are also preserved in this package.
 - Because of the adoption of 'multiApply' regime, the functions work well with package 'startR'. 
diff --git a/R/Clim.R b/R/Clim.R
index d879fc4f4fdba572efeb4308ded28ea1ce41b799..7d4ab88d68062e088fdd0bd055ae24dee62cf737 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -160,29 +160,105 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
   order_obs <- match(name_exp, name_obs)
   obs <- Reorder(obs, order_obs)
   
-
+#print('dim(exp) 1')
+#print(dim(exp))
+#print('dim(obs) 1')
+#print(dim(obs))
+#print('exp[1,,,1,1,1] 1')
+#print(exp[1,,,1,1,1])
+#print('obs[1,1,,,1,1] 1')
+#print(obs[1,1,,,1,1])
   ###############################
   # Calculate Clim
 
   #----------------------------------
   #  Remove all sdate if not complete along dat_dim
 
-    pos <- rep(0, length(dat_dim))
-    for (i in 1:length(dat_dim)) {  #[dat, sdate]
-      ## dat_dim: [dataset, member]
-      pos[i] <- which(names(dim(obs)) == dat_dim[i])
+  # First, for loop through dat_dim to create per-paired NA array
+  dat_dim_ind <- match(dat_dim, name_exp)
+  na_array <- array(0, dim = dim(exp)[-dat_dim_ind])
+
+  na_array <- na_array + MeanDims(exp, dat_dim_ind, na.rm = FALSE)
+  na_array <- na_array + MeanDims(obs, dat_dim_ind, na.rm = FALSE)
+#print('na_array')
+#print(dim(na_array))
+#print(summary(na_array))
+ 
+#  na_array <- as.logical(na_array)  # TRUE is non-NA, NA is NA
+
+  exp_dat_dim_length <- dim(exp)[dat_dim_ind]  #[dat = 2, memb = 3, dat = 2] 
+  exp_dat_dim_length <- lapply(lapply(exp_dat_dim_length, seq, 1), sort)   # $dat: 1 2; $memb: 1 2 3
+  exp_ind_list_for_Subset <- expand.grid(exp_dat_dim_length)
+#print('exp_ind_list_for_Subset')
+#print(exp_ind_list_for_Subset)
+  obs_dat_dim_length <- dim(obs)[dat_dim_ind]  
+  obs_dat_dim_length <- lapply(lapply(obs_dat_dim_length, seq, 1), sort)
+  obs_ind_list_for_Subset <- expand.grid(obs_dat_dim_length)
+#print('obs_ind_list_for_Subset')
+#print(obs_ind_list_for_Subset)
+
+  list_exp <- vector('list', length = length(exp_ind_list_for_Subset[, 1]))
+  list_obs <- vector('list', length = length(obs_ind_list_for_Subset[, 1]))
+
+  # exp
+  for (kk in 1:length(list_exp)) {
+    tmp_exp <- ClimProjDiags::Subset(exp, dat_dim, as.list(exp_ind_list_for_Subset[kk, ]), drop = 'selected')
+    tmp_exp[which(is.na(na_array))] <- NA
+    list_exp[[kk]] <- tmp_exp
+  }
+  # obs
+  for (kk in 1:length(list_obs)) {
+    tmp_obs <- ClimProjDiags::Subset(obs, dat_dim, as.list(obs_ind_list_for_Subset[kk, ]), drop = 'selected')
+    tmp_obs[which(is.na(na_array))] <- NA
+    list_obs[[kk]] <- tmp_obs
+  }
+
+#print('str(list_exp)')
+#print(str(list_exp))
+#print('str(list_obs)')
+#print(str(list_obs))
+
+  exp <- array(unlist(list_exp), dim = c(dim(exp)[-dat_dim_ind], dim(exp)[dat_dim_ind]))
+  obs <- array(unlist(list_obs), dim = c(dim(obs)[-dat_dim_ind], dim(obs)[dat_dim_ind]))
+  
+  # Reorder exp and obs back to the original one
+  ori_order <- match(name_exp, names(dim(exp)))
+  exp <- Reorder(exp, ori_order)
+  obs <- Reorder(obs, ori_order)
+#print('dim(exp) 2')
+#print(dim(exp))
+#print('dim(obs) 2')
+#print(dim(obs))
+#print('exp[1,,,1,1,1] 2')
+#print(exp[1,,,1,1,1])
+#print('obs[1,1,,,1,1] 2')
+#print(obs[1,1,,,1,1])
+#--------------------------------------------
+
+    if (ftime_dim %in% names(dim(exp))) {
+      pos <- rep(0, length(time_dim) + length(ftime_dim))
+      ## time_dim: [sdate]; ftime_dim: [ftime]
+      pos[1] <- which(names(dim(obs)) == time_dim)
+      pos[2] <- which(names(dim(obs)) == ftime_dim)
+    } else {
+      pos <- rep(0, length(time_dim))
+      ## time_dim: [sdate]
+      pos[1] <- which(names(dim(obs)) == time_dim)
     }
-    outrows_exp <- MeanDims(exp, pos, na.rm = FALSE) + 
-                   MeanDims(obs, pos, na.rm = FALSE)
+
+    avg_dims <- c(1:length(dim(obs)))[-pos]
+    outrows_exp <- MeanDims(exp, avg_dims, na.rm = FALSE) + 
+                   MeanDims(obs, avg_dims, na.rm = FALSE)
     outrows_obs <- outrows_exp
+    # outrows: [sdate]
 
-    for (i in 1:length(pos)) {
-      outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]])
-      outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]])
+    for (i in 1:length(avg_dims)) {
+      outrows_exp <- InsertDim(outrows_exp, avg_dims[i], dim(exp)[avg_dims[i]])
+      outrows_obs <- InsertDim(outrows_obs, avg_dims[i], dim(obs)[avg_dims[i]])
     }
+
     exp[which(is.na(outrows_exp))] <- NA
     obs[which(is.na(outrows_obs))] <- NA
-
   #-----------------------------------
 
   if (method == 'clim') {
diff --git a/tests/testthat/test-Clim.R b/tests/testthat/test-Clim.R
index 1d2ff3e6b013cdefd1cd3375a34193d2bbc0545d..d96a42a337ff0b9696b2c9f47c963021385146a3 100644
--- a/tests/testthat/test-Clim.R
+++ b/tests/testthat/test-Clim.R
@@ -3,20 +3,16 @@ context("s2dv::Clim tests")
 ##############################################
   # dat1
   set.seed(1)
-  exp1 <- array(rnorm(360), dim = c(dataset = 1, member = 3, sdate = 5, 
-                                    ftime = 3, lon = 2, lat = 4))
+  exp1 <- array(rnorm(360), dim = c(dataset = 2, member = 3, sdate = 5, 
+                                    ftime = 3, lon = 2, lat = 2))
   set.seed(2)
-  obs1 <- array(rnorm(120), dim = c(dataset = 1, member = 1,
-                                   ftime = 3, lon = 2, lat = 4, sdate = 5))
+  obs1 <- array(rnorm(60), dim = c(dataset = 1, member = 1,
+                                   ftime = 3, lon = 2, lat = 2, sdate = 5))
   # dat2
   exp2 <- exp1
-  set.seed(1)
-  na <- floor(runif(5, min = 1, max = 360))
-  exp2[na] <- NA
+  exp2[1, 2, , 1, 1, 1] <- NA
   obs2 <- obs1
-  set.seed(2)
-  na <- floor(runif(30, min = 1, max = 120))
-  obs2[na] <- NA
+  obs2[1, 1, 1, 1, 1, 2] <- NA
 
   # dat3
   set.seed(1)
@@ -112,59 +108,59 @@ test_that("2. Output checks: dat1", {
 
   expect_equal(
     dim(Clim(exp1, obs1)$clim_exp),
-    c(dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4)
+    c(dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, memb = FALSE)$clim_exp),
-    c(dataset = 1, ftime = 3, lon = 2, lat = 4)
+    c(dataset = 2, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, time_dim = 'lon')$clim_exp),
-    c(dataset = 1, member = 3, sdate = 5, ftime = 3, lat = 4)
+    c(dataset = 2, member = 3, sdate = 5, ftime = 3, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, method = 'kharin')$clim_exp),
-    c(sdate = 5, dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4)
+    c(sdate = 5, dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, method = 'NDV')$clim_exp),
-    c(sdate = 5, dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4)
+    c(sdate = 5, dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1)$clim_obs),
-    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4)
+    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, method = 'kharin')$clim_obs),
-    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4)
+    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     dim(Clim(exp1, obs1, method = 'NDV')$clim_obs),
-    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4)
+    c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2)
   )
   expect_equal(
     (Clim(exp1, obs1)$clim_obs)[1:5],
-     c(0.14831161, -0.60462627, 0.06609153, -0.23378059, 0.50553522),
+     c(-0.6389627, -1.0815757, 0.4899925, -0.4436428, 0.1289446),
      tolerance = 0.001
   )
   expect_equal(
     (Clim(exp1, obs1, memb = FALSE)$clim_exp)[1:5],
-     c(0.10084284, 0.06407350, 0.09028584, 0.17526332, 0.18004387),
+     c(0.24870405, -0.08378771, 0.23637962, 0.02916955, 0.18643424),
      tolerance = 0.001
   )
   expect_equal(
     max(Clim(exp1, obs1)$clim_exp, na.rm = T),
-    1.026186,
+    1.065416,
     tolerance = 0.001
   )
   expect_equal(
     max(Clim(exp1, obs1, method = 'kharin')$clim_exp, na.rm = T),
-    2.282634,
+    2.634209,
     tolerance = 0.001
   )
   expect_equal(
     min(Clim(exp1, obs1, method = 'NDV')$clim_exp, na.rm = T),
-    -4.025745,
+    -3.156395,
     tolerance = 0.001
   )
 
@@ -175,21 +171,20 @@ test_that("3. Output checks: dat2", {
 
   expect_equal(
     (Clim(exp2, obs2)$clim_obs)[1:5],
-     c(0.23142987, -0.60462627, -0.03669491, -0.14193572, 0.61163024),
+     c(NaN, -1.0815757, 0.4899925, NaN, 0.1289446),
      tolerance = 0.001
   )
   expect_equal(
-    (Clim(exp2, obs2)$clim_exp)[1:5],
-     c(0.01054951, -0.04744191, -0.03533071, 0.14149945, 0.02359945),
-     tolerance = 0.001
+    all(is.na(Clim(exp2, obs2)$clim_exp[, , 1, 2, 2])),
+    TRUE
   )
   expect_equal(
     length(which(is.na(Clim(exp2, obs2, na.rm = FALSE)$clim_exp))),
-    57
+    24
   )
   expect_equal(
     length(which(is.na(Clim(exp2, obs2, na.rm = FALSE)$clim_obs))),
-    19
+    8
   )
   expect_equal(
     length(which(is.na(Clim(exp2, obs2)$clim_obs))),
@@ -212,9 +207,5 @@ test_that("4. Output checks: dat3", {
      tolerance = 0.00001
   )
 
-
-
-
-
 })