diff --git a/R/DiffCorr.R b/R/DiffCorr.R
index c250814ce859739f4ad5a9e79e92aa14d8cd1fb4..95996b68e011a1be8834297889797cdf71b9bd95 100644
--- a/R/DiffCorr.R
+++ b/R/DiffCorr.R
@@ -4,10 +4,10 @@
 #'Positive values of the correlation difference indicate that the forecast is 
 #'more skillful than the reference forecast, while negative values mean that 
 #'the reference forecast is more skillful. The statistical significance of the 
-#'correlation differences is computed with a one-sided test for equality of 
-#'dependent correlation coefficients (Steiger, 1980; Siegert et al., 2017) using
-#'effective degrees of freedom to account for the autocorrelation of the time 
-#'series (von Storch and Zwiers, 1999).
+#'correlation differences is computed with a one-sided or two-sided test for 
+#'equality of dependent correlation coefficients (Steiger, 1980; Siegert et al.,
+#'2017) using effective degrees of freedom to account for the autocorrelation of
+#'the time series (von Storch and Zwiers, 1999).
 #'
 #'@param exp A named numerical array of the forecast data with at least time 
 #'  dimension. 
@@ -38,6 +38,12 @@
 #'  steps with no missing values in all "exp", "ref", and "obs" will be used. If
 #'  "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any
 #'  NA. The default value is "return.na".
+#'@param test.type A character string indicating the type of significance test. 
+#'  It can be "two-sided" (to assess whether the skill of "exp" and "ref" are 
+#'  significantly different with a z-test on Fisher z-transformed correlation 
+#'  coefficients, Hinkel et al, 1988) or "one-sided" (to assess whether the 
+#'  skill of "exp" is significantly higher than that of "ref" with the test 
+#'  described in Steiger, 1980). The default value is "two-sided".
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
@@ -66,13 +72,14 @@
 #' exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
 #' obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
 #' ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
-#' res <- DiffCorr(exp, obs, ref, memb_dim = 'member')
+#' res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', alpha = 0.05)
+#' res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided')
 #'
 #'@import multiApply
 #'@export
 DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate', 
                      memb_dim = NULL, method = 'pearson', alpha = NULL, 
-                     handle.na = 'return.na', ncores = NULL) {
+                     handle.na = 'return.na', test.type = "two-sided", ncores = NULL) {
   
   # Check inputs
   ## exp, ref, and obs (1)
@@ -143,6 +150,14 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate',
   if (!handle.na %in% c('return.na', 'only.complete.triplets', 'na.fail')) {
     stop('Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".')
   }
+  ## test.type
+  if (!test.type %in% c('two-sided', 'one-sided')) {
+    stop("Parameter 'test.type' must be 'two-sided' or 'one-sided'.")
+  }
+  #NOTE: warning can be removed in the next release
+  warning("The default significance test has changed after s2dv_1.2.0. The default method is 'two-sided'.")
+
+  ## ncores
   if (!is.null(ncores)) {
     if (any(!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
         length(ncores) > 1)) {
@@ -184,47 +199,67 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate',
                                        ref = time_dim, N.eff = NULL),
                     output_dims = output_dims,
                     fun = .DiffCorr, method = method,
-                    alpha = alpha, handle.na = handle.na, ncores = ncores)
+                    alpha = alpha, handle.na = handle.na, 
+                    test.type = test.type, ncores = ncores)
   } else { 
     output <- Apply(data = list(exp = exp, obs = obs, ref = ref),
                     target_dims = list(exp = time_dim, obs = time_dim, 
                                        ref = time_dim), 
                     output_dims = output_dims, N.eff = N.eff,
                     fun = .DiffCorr, method = method, 
-                    alpha = alpha, handle.na = handle.na, ncores = ncores)
+                    alpha = alpha, handle.na = handle.na, 
+                    test.type = test.type, ncores = ncores)
   }
 
   return(output)
 }
 
-.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, handle.na = 'return.na') {
-
-  .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL) {
+.DiffCorr <- function(exp, obs, ref, N.eff = NA, method = 'pearson', alpha = NULL, 
+                      handle.na = 'return.na', test.type = 'two.sided') {
 
+  .diff.corr <- function(exp, obs, ref, method = 'pearson', N.eff = NA, alpha = NULL, test.type = 'two.sided') {
     # Correlation difference
     cor.exp <- cor(x = exp, y = obs, method = method)
     cor.ref <- cor(x = ref, y = obs, method = method)
     output <- NULL
     output$diff.corr <- cor.exp - cor.ref
-
-    # Significance with one-sided test for equality of dependent correlation coefficients (Steiger, 1980)
-    r12 <- cor.exp
-    r13 <- cor.ref
-    r23 <- cor(exp, ref)
     if (is.na(N.eff)) {
       N.eff <- .Eno(x = obs, na.action = na.pass) ## effective degrees of freedom
     }
-    R <- (1 - r12 * r12 - r13 * r13 - r23 * r23) + 2 * r12 * r13 * r23
-    t <- (r12 - r13) * sqrt((N.eff - 1) * (1 + r23) / (2 * ((N.eff - 1) / (N.eff - 3)) * R + 0.25 * (r12 + r13)^2 * (1 - r23)^3))
-    p.value <- 1 - pt(t, df = N.eff - 3)
-    if (is.null(alpha)) {
-      output$p.val <- p.value
+
+    if (test.type == 'one-sided') {
+      # Significance with one-sided test for equality of dependent correlation coefficients (Steiger, 1980)
+      r12 <- cor.exp
+      r13 <- cor.ref
+      r23 <- cor(exp, ref)
+      R <- (1 - r12 * r12 - r13 * r13 - r23 * r23) + 2 * r12 * r13 * r23
+      t <- (r12 - r13) * sqrt((N.eff - 1) * (1 + r23) / (2 * ((N.eff - 1) / (N.eff - 3)) * R + 0.25 * (r12 + r13)^2 * (1 - r23)^3))
+      p.value <- pt(t, df = N.eff - 3, lower.tail = FALSE)
+      if (is.null(alpha)) {
+        output$p.val <- p.value
+      } else {
+        output$sign <- ifelse(!is.na(p.value) & p.value <= alpha & output$diff.corr > 0, TRUE, FALSE)
+      }
+
+    } else if (test.type == 'two-sided') {
+      # Significance with two-sided z-test on Fisher z-transformed correlation coefficients (Hinkel et al, 1988)
+      Z <- abs(0.5 * log((1 + cor.exp) / (1 - cor.exp)) - 0.5 * log((1 + cor.ref) / (1 - cor.ref))) / sqrt(1/ (N.eff - 3) + 1/(N.eff - 3))
+      # sign <- Z >= qnorm(p = alpha/2, lower.tail = FALSE)
+      p.value <- pnorm(q = Z, lower.tail = FALSE)
+
+      if (is.null(alpha)) {
+        output$p.val <- p.value
+      } else {
+        output$sign <- ifelse(!is.na(p.value) & p.value <= alpha / 2, TRUE, FALSE)
+      }
+ 
     } else {
-      output$sign <- ifelse(!is.na(p.value) & p.value <= alpha, TRUE, FALSE)
+      stop("Parameter 'test.type' is not supported.")
     }
+
     return(output)
   }
-
+  
   #==================================================
   
   if (anyNA(exp) | anyNA(obs) | anyNA(ref)) { ## There are NAs 
@@ -237,7 +272,7 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate',
       ref <- ref[!nna]
 
       output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, 
-                           N.eff = N.eff, alpha = alpha)
+                           N.eff = N.eff, alpha = alpha, test.type = test.type)
       
     } else if (handle.na == 'return.na') {
       # Data contain NA, return NAs directly without passing to .diff.corr
@@ -250,7 +285,8 @@ DiffCorr <- function(exp, obs, ref, N.eff = NA, time_dim = 'sdate',
     
   } else { ## There is no NA  
     output <- .diff.corr(exp = exp, obs = obs, ref = ref, method = method, 
-                         N.eff = N.eff, alpha = alpha)
+                         N.eff = N.eff, alpha = alpha, test.type = test.type)
   }
+
   return(output)
 }
diff --git a/man/DiffCorr.Rd b/man/DiffCorr.Rd
index d8ff65cbeda9efc4a6a345610971527774851ad0..e0e5ec37ada8d7dc45d90ac83df33a340cc15c16 100644
--- a/man/DiffCorr.Rd
+++ b/man/DiffCorr.Rd
@@ -14,6 +14,7 @@ DiffCorr(
   method = "pearson",
   alpha = NULL,
   handle.na = "return.na",
+  test.type = "two-sided",
   ncores = NULL
 )
 }
@@ -56,6 +57,13 @@ steps with no missing values in all "exp", "ref", and "obs" will be used. If
 "na.fail", an error will arise if any of "exp", "ref", or "obs" contains any
 NA. The default value is "return.na".}
 
+\item{test.type}{A character string indicating the type of significance test. 
+It can be "two-sided" (to assess whether the skill of "exp" and "ref" are 
+significantly different with a z-test on Fisher z-transformed correlation 
+coefficients, Hinkel et al, 1988) or "one-sided" (to assess whether the 
+skill of "exp" is significantly higher than that of "ref" with the test 
+described in Steiger, 1980). The default value is "two-sided".}
+
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
@@ -81,16 +89,17 @@ Compute the correlation difference between two deterministic forecasts.
 Positive values of the correlation difference indicate that the forecast is 
 more skillful than the reference forecast, while negative values mean that 
 the reference forecast is more skillful. The statistical significance of the 
-correlation differences is computed with a one-sided test for equality of 
-dependent correlation coefficients (Steiger, 1980; Siegert et al., 2017) using
-effective degrees of freedom to account for the autocorrelation of the time 
-series (von Storch and Zwiers, 1999).
+correlation differences is computed with a one-sided or two-sided test for 
+equality of dependent correlation coefficients (Steiger, 1980; Siegert et al.,
+2017) using effective degrees of freedom to account for the autocorrelation of
+the time series (von Storch and Zwiers, 1999).
 }
 \examples{
 exp <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
 obs <- array(rnorm(1000), dim = c(lat = 3, lon = 2, sdate = 50))
 ref <- array(rnorm(1000), dim = c(lat = 3, lon = 2, member = 10, sdate = 50))
-res <- DiffCorr(exp, obs, ref, memb_dim = 'member')
+res_two.sided_sign <- DiffCorr(exp, obs, ref, memb_dim = 'member', alpha = 0.05)
+res_one.sided_pval <- DiffCorr(exp, obs, ref, memb_dim = 'member', test.type = 'one-sided')
 
 }
 \references{
diff --git a/tests/testthat/test-DiffCorr.R b/tests/testthat/test-DiffCorr.R
index e0834e1d05196c34d6acdbe84b2359fbfee7b5ba..b8c8029ebbc41e09156d485a69b6c407418654ed 100644
--- a/tests/testthat/test-DiffCorr.R
+++ b/tests/testthat/test-DiffCorr.R
@@ -20,7 +20,6 @@ set.seed(3)
 obs2 <- array(rnorm(10), dim = c(sdate = 10))
 
 
-
 ##############################################
 test_that("1. Input checks", {
   # exp and obs (1)
@@ -89,9 +88,14 @@ test_that("1. Input checks", {
   DiffCorr(exp2, obs2, ref2, handle.na = TRUE),
   'Parameter "handle.na" must be "return.na", "only.complete.triplets" or "na.fail".'
   )
+  # test.type
+  expect_error(
+  DiffCorr(exp2, obs2, ref2, test.type = "two.sided"),
+  "Parameter 'test.type' must be 'two-sided' or 'one-sided'."
+  )
   # ncores
   expect_error(
-  DiffCorr(exp2, obs2, ref2, ncores = 1.5),
+  suppressWarnings(DiffCorr(exp2, obs2, ref2, ncores = 1.5)),
   'Parameter "ncores" must be either NULL or a positive integer.'
   )
 
@@ -99,6 +103,8 @@ test_that("1. Input checks", {
 
 ##############################################
 test_that("2. Output checks: dat1", {
+suppressWarnings({
+
 expect_equal(
 names(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb')),
 c("diff.corr", "p.val")
@@ -118,7 +124,7 @@ tolerance = 0.0001
 )
 expect_equal(
 as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb')$p),
-c(0.26166060, 0.15899774, 0.39264452, 0.27959883, 0.34736305, 0.07479832),
+c(0.28984782, 0.14395946, 0.41374677, 0.32253784, 0.32677872, 0.03659053),
 tolerance = 0.0001
 )
 expect_equal(
@@ -131,17 +137,25 @@ c(0.27347087, 0.50556882, 0.08855968, 0.24199701, 0.22935182, 0.88336336),
 tolerance = 0.0001
 )
 expect_equal(
+as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided")$diff.corr),
+as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "two-sided")$diff.corr)
+)
+expect_equal(
 as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05)$sign),
 rep(FALSE, 6)
 )
 expect_equal(
+as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', alpha = 0.05, test.type = "one-sided")$sign),
+rep(FALSE, 6)
+)
+expect_equal(
 suppressWarnings(as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$diff.corr)),
 c(0.07272727, 0.54545455, 0.10909091, -0.01212121, -0.03636364, 1.01818182),
 tolerance = 0.0001
 )
 expect_equal(
 suppressWarnings(as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', method = "spearman")$p)),
-c(0.4358970, 0.1341575, 0.3448977, 0.5114262, 0.5264872, 0.0437861),
+c(0.44479470, 0.11854232, 0.36667620, 0.49095264, 0.46962814, 0.01714977),
 tolerance = 0.0001
 )
 expect_equal(
@@ -151,7 +165,7 @@ tolerance = 0.0001
 )
 expect_equal(
 as.vector(DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', N.eff = Neff1)$p),
-c(0.27841537, 0.15899774, 0.40096749, 0.27959883, 0.35889690, 0.07479832),
+c(0.30406438, 0.14395946, 0.42005522, 0.32253784, 0.33887632, 0.03659053),
 tolerance = 0.0001
 )
 
@@ -181,12 +195,15 @@ DiffCorr(exp1, obs1, ref1, memb_dim = 'memb', handle.na = 'only.complete.triplet
 "There is no complete set of forecasts and observations."
 )
 
-
+}) # suppressWarnings
 })
 
 ##############################################
 
 test_that("3. Output checks: dat2", {
+
+suppressWarnings({
+
 expect_equal(
 names(DiffCorr(exp2, obs2, ref2)),
 c("diff.corr", "p.val")
@@ -201,13 +218,27 @@ NULL
 )
 expect_equal(
 DiffCorr(exp2, obs2, ref2)$p,
+0.3201317,
+tolerance = 0.0001
+)
+expect_equal(
+DiffCorr(exp2, obs2, ref2, test.type = 'one-sided')$p,
 0.6577392,
 tolerance = 0.0001
 )
 expect_equal(
+DiffCorr(exp2, obs2, ref2, test.type = 'one-sided', alpha = 0.5)$sign,
+FALSE
+)
+expect_equal(
+DiffCorr(exp2, obs2, ref2, test.type = 'two-sided', alpha = 0.65)$sign,
+TRUE
+)
+expect_equal(
 DiffCorr(exp2, obs2, ref2)$diff,
 -0.2434725,
 tolerance = 0.0001
 )
 
+}) # suppressWarnings
 })