diff --git a/DESCRIPTION b/DESCRIPTION
index ec571f4ce6bd0ba036ae274b7e085237d6c19863..c5a8139e8707202a51025c59db4cf38fd360a2f2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -50,5 +50,5 @@ BugReports: https://earth.bsc.es/gitlab/es/s2dv/-/issues
 LazyData: true
 SystemRequirements: cdo
 Encoding: UTF-8
-RoxygenNote: 7.2.0
+RoxygenNote: 7.2.3
 Config/testthat/edition: 3
diff --git a/R/NAO.R b/R/NAO.R
index 5586985e20add9ee9e1026d477001293b4af68d9..207af2b32c908d137ac5bf3ab34f7703e520af7d 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -7,11 +7,12 @@
 #'anomalies onto the EOF pattern of the other years of the forecast. 
 #'By default (ftime_avg = 2:4), NAO() computes the NAO index for 1-month 
 #'lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
-#'cross-validated PCs of the NAO index for forecast (exp) and observations 
-#'(obs) based on the leading EOF pattern.
+#'cross-validated PCs of the NAO index for hindcast (exp) and observations 
+#'(obs) based on the leading EOF pattern, or, if forecast (exp_cor) is provided,
+#'the NAO index for forecast and the corresponding data (exp and obs).
 #'
 #'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
-#'  forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+#'  hindcast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
 #'  dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
 #'  If only NAO of observational data needs to be computed, this parameter can
 #'  be left to NULL. The default value is NULL.
@@ -20,6 +21,12 @@
 #'  dimensions 'time_dim', 'ftime_dim', and 'space_dim' at least.
 #'  If only NAO of experimental data needs to be computed, this parameter can 
 #'  be left to NULL. The default value is NULL.
+#'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
+#'  forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+#'  dimension 'time_dim' of length 1 (as in the case of an operational 
+#'  forecast), 'memb_dim', 'ftime_dim', and 'space_dim' at least.
+#'  If only NAO of reference period needs to be computed, this parameter can
+#'  be left to NULL. The default value is NULL. 
 #'@param lat A vector of the latitudes of 'exp' and 'obs'.
 #'@param lon A vector of the longitudes of 'exp' and 'obs'.
 #'@param time_dim A character string indicating the name of the time dimension
@@ -37,25 +44,32 @@
 #'  value is 2:4, i.e., from 2nd to 4th forecast time steps.
 #'@param obsproj A logical value indicating whether to compute the NAO index by
 #'  projecting the forecast anomalies onto the leading EOF of observational 
-#'  reference (TRUE) or compute the NAO by first computing the leading 
-#'  EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
-#'  year you are evaluating out), and then projecting forecast anomalies onto 
-#'  this EOF (FALSE). The default value is TRUE.
+#'  reference (TRUE, default) or compute the NAO by first computing the leading
+#'  EOF of the forecast anomalies (in cross-validation mode, i.e. leave the 
+#'  evaluated year out), then projecting forecast anomalies onto this EOF
+#'  (FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE
+#'  and 'exp' will be used when obsproj is FALSE, and no cross-validation is
+#'  applied.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
 #'@return 
-#'A list which contains:
+#'A list which contains some of the following items depending on the data inputs:
 #'\item{exp}{
-#'  A numeric array of forecast NAO index in verification format with the same 
+#'  A numeric array of hindcast NAO index in verification format with the same 
 #'  dimensions as 'exp' except space_dim and ftime_dim. If ftime_avg is NULL, 
 #'  ftime_dim remains.
 #'  }
 #'\item{obs}{
-#'  A numeric array of observed NAO index in verification format with the same
+#'  A numeric array of observation NAO index in verification format with the same
 #'  dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL,
 #'  ftime_dim remains.
 #'}
+#'\item{exp_cor}{
+#'  A numeric array of forecast NAO index in verification format with the same 
+#'  dimensions as 'exp_cor' except space_dim and ftime_dim. If ftime_avg is NULL, 
+#'  ftime_dim remains.
+#'  }
 #'
 #'@references
 #'Doblas-Reyes, F.J., Pavan, V. and Stephenson, D. (2003). The skill of 
@@ -73,6 +87,8 @@
 #'lon <- seq(-80, 40, length.out = 9) 
 #'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 #'
+#'exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
+#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE)
 #'# plot the NAO index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
@@ -84,13 +100,12 @@
 #'@import multiApply
 #'@importFrom ClimProjDiags Subset
 #'@export
-NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
+NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sdate',
                 memb_dim = 'member', space_dim = c('lat', 'lon'),
                 ftime_dim = 'ftime', ftime_avg = 2:4, 
                 obsproj = TRUE, ncores = NULL) {
-
   # Check inputs 
-  ## exp and obs (1)
+  ## exp, obs, and exp_cor (1)
   if (is.null(obs) & is.null(exp)) {
     stop("Parameter 'exp' and 'obs' cannot both be NULL.")
   }
@@ -118,20 +133,45 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'obs' must have dimension names.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (!is.numeric(exp_cor)) {
+      stop("Parameter 'exp_cor' must be a numeric array.")
+    }
+    if (is.null(dim(exp_cor))) {
+      stop(paste0("Parameter 'exp_cor' must have at least dimensions ",
+                  "time_dim, memb_dim, space_dim, and ftime_dim."))
+    }
+    if (any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
+      stop("Parameter 'exp_cor' must have dimension names.")
+    }
+    if (is.null(exp) || is.null(obs)) {
+      stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not provided.")
+    }
+  }
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) > 1) {
     stop("Parameter 'time_dim' must be a character string.")
   }
   if (!is.null(exp)) {
     if (!time_dim %in% names(dim(exp))) {
-      stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+      stop("Parameter 'time_dim' is not found in 'exp' dimension.")
     }
   }
   if (!is.null(obs)) {
     if (!time_dim %in% names(dim(obs))) {
-      stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+      stop("Parameter 'time_dim' is not found in 'obs' dimension.")
+    }
+  }
+  if (!is.null(exp_cor)) {
+    if (!time_dim %in% names(dim(exp_cor))) {
+      stop("Parameter 'time_dim' is not found in 'exp_cor' dimension.")
+    }
+    if (dim(exp_cor)[time_dim] > 1) {
+      stop("Parameter 'exp_cor' is expected to have length 1 in ",
+           time_dim, "dimension.")
     }
   }
+
   ## memb_dim
   if (!is.character(memb_dim) | length(memb_dim) > 1) {
     stop("Parameter 'memb_dim' must be a character string.")
@@ -141,6 +181,7 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
     }
   }
+  add_member_back <- FALSE
   if (!is.null(obs)) {
     if (memb_dim %in% names(dim(obs))) {
       if (dim(obs)[memb_dim] != 1) {
@@ -149,8 +190,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
         add_member_back <- TRUE
         obs <- ClimProjDiags::Subset(obs, memb_dim, 1, drop = 'selected')
       }
-    } else {
-      add_member_back <- FALSE
+    }
+  }
+  if (!is.null(exp_cor)) {
+    if (!memb_dim %in% names(dim(exp_cor))) {
+      stop("Parameter 'memb_dim' is not found in 'exp_cor' dimension.")
     }
   }
   ## space_dim
@@ -167,6 +211,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (any(!space_dim %in% names(dim(exp_cor)))) {
+      stop("Parameter 'space_dim' is not found in 'exp_cor' dimensions.")
+    }
+  }
   ## ftime_dim
   if (!is.character(ftime_dim) | length(ftime_dim) > 1) {
     stop("Parameter 'ftime_dim' must be a character string.")
@@ -181,7 +230,13 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (!ftime_dim %in% names(dim(exp_cor))) {
+      stop("Parameter 'ftime_dim' is not found in 'exp_cor' dimensions.")
+    }
+  }
   ## exp and obs (2)
+  #TODO: Add checks for exp_cor
   if (!is.null(exp) & !is.null(obs)) {
     name_exp <- sort(names(dim(exp)))
     name_obs <- sort(names(dim(obs)))
@@ -208,11 +263,17 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       if (max(ftime_avg) > dim(exp)[ftime_dim] | min(ftime_avg) < 1) {
         stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
       }
-    } else {
+    } 
+    if (!is.null(obs)) {
       if (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) {
         stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
       }
     }
+    if (!is.null(exp_cor)) {
+      if (max(ftime_avg) > dim(exp_cor)[ftime_dim] | min(ftime_avg) < 1) {
+        stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
+      }
+    }
   }
   ## sdate >= 2
   if (!is.null(exp)) {
@@ -234,7 +295,8 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'lon' must be a numeric vector with the same ",
            "length as the longitude dimension of 'exp' and 'obs'.")
     }
-  } else {
+  } 
+  if (!is.null(obs)) {
     if (!is.numeric(lat) | length(lat) != dim(obs)[space_dim[1]]) {
       stop("Parameter 'lat' must be a numeric vector with the same ",
            "length as the latitude dimension of 'exp' and 'obs'.")
@@ -244,6 +306,16 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
            "length as the longitude dimension of 'exp' and 'obs'.")
     }
   } 
+  if (!is.null(exp_cor)) {
+    if (!is.numeric(lat) | length(lat) != dim(exp_cor)[space_dim[1]]) {
+      stop("Parameter 'lat' must be a numeric vector with the same ",
+           "length as the latitude dimension of 'exp_cor'.")
+    }
+    if (!is.numeric(lon) | length(lon) != dim(exp_cor)[space_dim[2]]) {
+      stop("Parameter 'lon' must be a numeric vector with the same ",
+           "length as the longitude dimension of 'exp_cor'.")
+    }
+  }
   stop_needed <- FALSE
   if (max(lat) > 80 | min(lat) < 20) {
     stop_needed <- TRUE
@@ -279,8 +351,8 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
     if (is.null(obs)) {
       stop("Parameter 'obsproj' set to TRUE but no 'obs' provided.")
     }
-    if (is.null(exp)) {
-      .warning("parameter 'obsproj' set to TRUE but no 'exp' provided.")
+    if (is.null(exp) & is.null(exp_cor)) {
+      .warning("parameter 'obsproj' set to TRUE but no 'exp' nor 'exp_cor' provided.")
     }
   }
   ## ncores
@@ -303,46 +375,63 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
       obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
     }
+    if (!is.null(exp_cor)) {
+      exp_cor_sub <- ClimProjDiags::Subset(exp_cor, ftime_dim, ftime_avg, drop = FALSE)
+      exp_cor <- MeanDims(exp_cor_sub, ftime_dim, na.rm = TRUE)
+    }
   }
 
   # wght
   wght <- array(sqrt(cos(lat * pi / 180)), dim = c(length(lat), length(lon)))
-
-  if (!is.null(exp) & !is.null(obs)) {
-    res <- Apply(list(exp, obs),
+  if (is.null(exp_cor)) {
+    if (!is.null(exp) & !is.null(obs)) {
+        res <- Apply(list(exp, obs),
+                     target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                        obs = c(time_dim, space_dim)),
+                     fun = .NAO,
+                     lat = lat, wght = wght, 
+                     obsproj = obsproj, add_member_back = add_member_back,
+                     ncores = ncores)
+    } else if (!is.null(exp)) {
+      res <- Apply(list(exp = exp),
+                   target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
+                   fun = .NAO,
+                   lat = lat, wght = wght, obs = NULL, 
+                   obsproj = obsproj, add_member_back = FALSE,
+                   ncores = ncores)
+    } else if (!is.null(obs)) {
+      if (add_member_back) {
+        output_dims <- list(obs = c(time_dim, memb_dim))
+      } else {
+        output_dims <- list(obs = time_dim)
+      }
+      res <- Apply(list(obs = obs),
+                   target_dims = list(obs = c(time_dim, space_dim)),
+                   output_dims = output_dims,
+                   fun = .NAO,
+                   lat = lat, wght = wght, exp = NULL, 
+                   obsproj = obsproj, add_member_back = add_member_back,
+                   ncores = ncores)
+    }
+  } else { # exp_cor provided
+    res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor),
                  target_dims = list(exp = c(memb_dim, time_dim, space_dim),
-                                    obs = c(time_dim, space_dim)),
+                                    obs = c(time_dim, space_dim),
+                                    exp_cor = c(memb_dim, time_dim, space_dim)),
                  fun = .NAO,
-                 lat = lat, wght = wght, 
-                 obsproj = obsproj, add_member_back = add_member_back,
-                 ncores = ncores)
-  } else if (!is.null(exp)) {
-    res <- Apply(list(exp = exp),
-                 target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
-                 fun = .NAO,
-                 lat = lat, wght = wght, obs = NULL, 
-                 obsproj = obsproj, add_member_back = FALSE,
-                 ncores = ncores)
-  } else if (!is.null(obs)) {
-    if (add_member_back) {
-      output_dims <- list(obs = c(time_dim, memb_dim))
-    } else {
-      output_dims <- list(obs = time_dim)
-    }
-    res <- Apply(list(obs = obs),
-                 target_dims = list(obs = c(time_dim, space_dim)),
-                 output_dims = output_dims,
-                 fun = .NAO,
-                 lat = lat, wght = wght, exp = NULL, 
+                 lat = lat, wght = wght,
                  obsproj = obsproj, add_member_back = add_member_back,
                  ncores = ncores)
   }
+
   return(res)
 }
 
-.NAO <- function(exp = NULL, obs = NULL, lat, wght, obsproj = TRUE, add_member_back = FALSE) {
+.NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, wght, obsproj = TRUE,
+                 add_member_back = FALSE) {
   # exp: [memb_exp, sdate, lat, lon]
   # obs: [sdate, lat, lon]
+  # exp_cor: [memb, sdate = 1, lat, lon]
   # wght: [lat, lon]
 
   if (!is.null(exp)) {
@@ -355,80 +444,130 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
     nlat <- dim(obs)[2]
     nlon <- dim(obs)[3]
   }
+  if (!is.null(exp_cor)) {
+    ntime_exp_cor <- dim(exp_cor)[2] # should be 1
+    nmemb_exp_cor <- dim(exp_cor)[1]
+  }
 
-  if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
-  if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
+  if (!is.null(obs)) nao_obs <- array(NA, dim = ntime)
+  if (!is.null(exp)) nao_exp <- array(NA, dim = c(ntime, nmemb_exp))
+  if (!is.null(exp_cor)) {
+    nao_exp_cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
+    #NOTE: The dimensions are flipped to fill in data correctly. Need to flip it back later.
+  }
 
-  for (tt in 1:ntime) {  #sdate
+  if (is.null(exp_cor)) { 
 
-    if (!is.null(obs)) {
-      ## Calculate observation EOF. Excluding one forecast start year.
-      obs_sub <- obs[(1:ntime)[-tt], , , drop = FALSE]
-      obs_EOF <- .EOF(obs_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+    for (tt in 1:ntime) {  # cross-validation
 
-      ## Correct polarity of pattern.
-      # dim(obs_EOF$EOFs): [mode, lat, lon]
-      if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-        obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
-#        obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+      if (!is.null(obs)) {
+        ## Calculate observation EOF. Excluding one forecast start year.
+        obs_sub <- obs[(1:ntime)[-tt], , , drop = FALSE]
+        EOF_obs <- .EOF(obs_sub, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+        ## Correct polarity of pattern
+        # EOF_obs: [mode = 1, lat, lon]
+        if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+          EOF_obs <- EOF_obs * (-1)
+        }
+        ## Project observed anomalies.
+        PF <- .ProjectField(obs, eof_mode = EOF_obs[1, , ], wght = wght)  # [sdate]
+        ## Keep PCs of excluded forecast start year. Fabian.
+        nao_obs[tt] <- PF[tt]
       }
-      ## Project observed anomalies.
-      PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-      ## Keep PCs of excluded forecast start year. Fabian.
-      NAOO.ver[tt] <- PF[tt]
-    }
 
-    if (!is.null(exp)) {
-      if (!obsproj) {
-        exp_sub <- exp[, (1:ntime)[-tt], , , drop = FALSE]
-        # Combine 'memb' and 'sdate' to calculate EOF
-        dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
-        exp_EOF <- .EOF(exp_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+      if (!is.null(exp)) {
+        if (!obsproj) {
+          exp_sub <- exp[, (1:ntime)[-tt], , , drop = FALSE]
+          # Combine 'memb' and 'sdate' to calculate EOF
+          dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
+          EOF_exp <- .EOF(exp_sub, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
 
-        ## Correct polarity of pattern.
-        ##NOTE: different from s2dverification, which doesn't use mean().
-#        if (0 < exp_EOF$EOFs[1, which.min(abs(lat - 65)), ]) {
-        if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-          exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
-#          exp_EOF$PCs <- exp_EOF$PCs * sign  # not used
-        }
+          ## Correct polarity of pattern
+          ##NOTE: different from s2dverification, which doesn't use mean().
+#          if (0 < EOF_exp[1, which.min(abs(lat - 65)), ]) {
+          if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+            EOF_exp <- EOF_exp * (-1)
+          }
 
-        ### Lines below could be simplified further by computing
-        ### ProjectField() only on the year of interest... (though this is
-        ### not vital). Lauriane
-        for (imemb in 1:nmemb_exp) {
-          PF <- .ProjectField(exp[imemb, , , ],
-                              eof_mode = exp_EOF$EOFs[1, , ],
-                              wght = wght) # [sdate, memb]
-          NAOF.ver[tt, imemb] <- PF[tt]
-        }    
-      } else {
-        ## Project forecast anomalies on obs EOF
-        for (imemb in 1:nmemb_exp) {
-          PF <- .ProjectField(exp[imemb, , , ],
-                              eof_mode = obs_EOF$EOFs[1, , ],
-                              wght = wght)  # [sdate]
-          NAOF.ver[tt, imemb] <- PF[tt]
+          ### Lines below could be simplified further by computing
+          ### ProjectField() only on the year of interest... (though this is
+          ### not vital). Lauriane
+          for (imemb in 1:nmemb_exp) {
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_exp[1, , ], wght = wght) # [sdate, memb]
+            nao_exp[tt, imemb] <- PF[tt]
+          }    
+        } else {
+          ## Project forecast anomalies on obs EOF
+          for (imemb in 1:nmemb_exp) {
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_obs[1, , ], wght = wght)  # [sdate]
+            nao_exp[tt, imemb] <- PF[tt]
+          }
         }
       }
+
+    }  # for loop sdate
+
+  } else { # exp_cor provided
+
+    ## Calculate observation EOF. Without cross-validation
+    EOF_obs <- .EOF(obs, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+    ## Correct polarity of pattern
+    # EOF_obs: [mode, lat, lon]
+    if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+      EOF_obs <- EOF_obs * (-1)
+    }
+    ## Project observed anomalies
+    PF <- .ProjectField(obs, eof_mode = EOF_obs, wght = wght)  # [mode = 1, sdate]
+    nao_obs[] <- PF[1, ]
+
+    if (!obsproj) {
+      # Calculate EOF_exp
+      tmp <- array(exp, dim = c(nmemb_exp * ntime, nlat, nlon))
+      EOF_exp <- .EOF(tmp, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+      ## Correct polarity of pattern
+      if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+        EOF_exp <- EOF_exp * (-1)
+      }
+      eof_mode_input <- EOF_exp[1, , ]
+    } else {
+      eof_mode_input <- EOF_obs[1, , ]
+    }
+
+    # Calculate NAO_exp
+    for (imemb in 1:dim(exp)[1]) {
+      exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, 
+                                       drop = 'selected')
+      PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght)  # [sdate]
+      nao_exp[ , imemb] <- PF
     }
 
-  }  # for loop sdate
+    # Calculate NAO_exp_cor
+    for (imemb in 1:dim(exp_cor)[1]) {
+      exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
+                                       drop = 'selected')
+      PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght)  # [sdate]
+      nao_exp_cor[, imemb] <- PF
+    }
 
+  }
   # add_member_back
   if (add_member_back) {
-    suppressWarnings({
-    NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1])
-    })
+    memb_dim_name <- ifelse(!is.null(names(dim(exp))[1]), names(dim(exp))[1], 'member')
+    nao_obs <- InsertDim(nao_obs, 2, 1, name = memb_dim_name)
   }
 
-  #NOTE: EOFs_obs is not returned because it's only the result of the last sdate
-  #     (It is returned in s2dverification.)
-  if (!is.null(exp) & !is.null(obs)) {
-    return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
-  } else if (!is.null(exp)) {
-    return(list(exp = NAOF.ver))
-  } else if (!is.null(obs)) {
-    return(list(obs = NAOO.ver))
+  # Return results
+  if (is.null(exp_cor)) {
+    res <- NULL
+    if (!is.null(exp)) {
+      res <- c(res, list(exp = nao_exp))
+    }
+    if (!is.null(obs)) {
+      res <- c(res, list(obs = nao_obs))
+    }
+    return(res)
+
+  } else {
+    return(list(exp = nao_exp, obs = nao_obs, exp_cor = nao_exp_cor))
   }
 }
diff --git a/R/PlotBoxWhisker.R b/R/PlotBoxWhisker.R
index 2ddcec0d59a7627432b388db555678dcaf4041e8..9b4c88e02836f68d32c98d4b92cc0c1d4f1f24cd 100644
--- a/R/PlotBoxWhisker.R
+++ b/R/PlotBoxWhisker.R
@@ -90,7 +90,7 @@
 #'ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 #'ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2])
 #'ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2])
-#'nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon)
+#'nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon)
 #'# Finally plot the nao index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/R/PlotStereoMap.R b/R/PlotStereoMap.R
index 4b4fbd22b93110e15e51ffe63071fac24f3b1d41..a2fdb85d27cc0faf2fd8d21160bfdfe9d1a3a6c0 100644
--- a/R/PlotStereoMap.R
+++ b/R/PlotStereoMap.R
@@ -170,8 +170,10 @@
 #'data <- matrix(rnorm(100 * 50), 100, 50)
 #'x <- seq(from = 0, to = 360, length.out = 100)
 #'y <- seq(from = -90, to = 90, length.out = 50)
+#'  \dontrun{
 #'PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50,
 #'              toptitle = "This is the title")
+#' }
 #'@import mapproj
 #'@importFrom grDevices dev.cur dev.new dev.off gray
 #'@importFrom stats median
diff --git a/R/ProjectField.R b/R/ProjectField.R
index 4fdf1892f2a8cd38a2f6492906700aa0cafb1cf2..70c8552c934becbcc2f2b365139bc1f2aa178543 100644
--- a/R/ProjectField.R
+++ b/R/ProjectField.R
@@ -229,7 +229,6 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon
   # wght: [lat, lon]
 
   ntime <- dim(ano)[1]
-
   if (length(dim(eof_mode)) == 2) {  # mode != NULL
     # Initialization of pc.ver.
     pc.ver <- array(NA, dim = ntime)  #[sdate]
diff --git a/man/NAO.Rd b/man/NAO.Rd
index 999fd75f6f45177e353006e6bd061e094f5678c2..8115fa22ff25b1a8793827ca5c9203bebb7499e1 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -7,6 +7,7 @@
 NAO(
   exp = NULL,
   obs = NULL,
+  exp_cor = NULL,
   lat,
   lon,
   time_dim = "sdate",
@@ -20,7 +21,7 @@ NAO(
 }
 \arguments{
 \item{exp}{A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
-forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+hindcast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
 dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
 If only NAO of observational data needs to be computed, this parameter can
 be left to NULL. The default value is NULL.}
@@ -31,6 +32,13 @@ dimensions 'time_dim', 'ftime_dim', and 'space_dim' at least.
 If only NAO of experimental data needs to be computed, this parameter can 
 be left to NULL. The default value is NULL.}
 
+\item{exp_cor}{A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
+forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+dimension 'time_dim' of length 1 (as in the case of an operational 
+forecast), 'memb_dim', 'ftime_dim', and 'space_dim' at least.
+If only NAO of reference period needs to be computed, this parameter can
+be left to NULL. The default value is NULL.}
+
 \item{lat}{A vector of the latitudes of 'exp' and 'obs'.}
 
 \item{lon}{A vector of the longitudes of 'exp' and 'obs'.}
@@ -55,26 +63,33 @@ value is 2:4, i.e., from 2nd to 4th forecast time steps.}
 
 \item{obsproj}{A logical value indicating whether to compute the NAO index by
 projecting the forecast anomalies onto the leading EOF of observational 
-reference (TRUE) or compute the NAO by first computing the leading 
-EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
-year you are evaluating out), and then projecting forecast anomalies onto 
-this EOF (FALSE). The default value is TRUE.}
+reference (TRUE, default) or compute the NAO by first computing the leading
+EOF of the forecast anomalies (in cross-validation mode, i.e. leave the 
+evaluated year out), then projecting forecast anomalies onto this EOF
+(FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE
+and 'exp' will be used when obsproj is FALSE, and no cross-validation is
+applied.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
 }
 \value{
-A list which contains:
+A list which contains some of the following items depending on the data inputs:
 \item{exp}{
- A numeric array of forecast NAO index in verification format with the same 
+ A numeric array of hindcast NAO index in verification format with the same 
  dimensions as 'exp' except space_dim and ftime_dim. If ftime_avg is NULL, 
  ftime_dim remains.
  }
 \item{obs}{
- A numeric array of observed NAO index in verification format with the same
+ A numeric array of observation NAO index in verification format with the same
  dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL,
  ftime_dim remains.
 }
+\item{exp_cor}{
+ A numeric array of forecast NAO index in verification format with the same 
+ dimensions as 'exp_cor' except space_dim and ftime_dim. If ftime_avg is NULL, 
+ ftime_dim remains.
+ }
 }
 \description{
 Compute the North Atlantic Oscillation (NAO) index based on the leading EOF 
@@ -84,8 +99,9 @@ observed anomalies onto the observed EOF pattern or the forecast
 anomalies onto the EOF pattern of the other years of the forecast. 
 By default (ftime_avg = 2:4), NAO() computes the NAO index for 1-month 
 lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
-cross-validated PCs of the NAO index for forecast (exp) and observations 
-(obs) based on the leading EOF pattern.
+cross-validated PCs of the NAO index for hindcast (exp) and observations 
+(obs) based on the leading EOF pattern, or, if forecast (exp_cor) is provided,
+the NAO index for forecast and the corresponding data (exp and obs).
 }
 \examples{
 # Make up synthetic data
@@ -97,6 +113,8 @@ lat <- seq(20, 80, length.out = 6)
 lon <- seq(-80, 40, length.out = 9) 
 nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 
+exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
+nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE)
 # plot the NAO index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/man/PlotBoxWhisker.Rd b/man/PlotBoxWhisker.Rd
index 9c5a3f48ab5ae30097f36dbc78bf141b0f648c9c..e8681607b4b639b9ccfe9b29651f618c2680b699 100644
--- a/man/PlotBoxWhisker.Rd
+++ b/man/PlotBoxWhisker.Rd
@@ -123,7 +123,7 @@ sampleData$lat[] <- c(20, 80)
 ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2])
 ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2])
-nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon)
+nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon)
 # Finally plot the nao index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/man/PlotStereoMap.Rd b/man/PlotStereoMap.Rd
index 1b7f166a161c3fff921d517bd4569fa54bd734bf..10e7503e3332b49deb233530a8a670118e313ddd 100644
--- a/man/PlotStereoMap.Rd
+++ b/man/PlotStereoMap.Rd
@@ -283,6 +283,8 @@ compatible with figure layouts if colour bar is disabled.
 data <- matrix(rnorm(100 * 50), 100, 50)
 x <- seq(from = 0, to = 360, length.out = 100)
 y <- seq(from = -90, to = 90, length.out = 50)
+ \dontrun{
 PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50,
              toptitle = "This is the title")
 }
+}
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index 91b994329ed8f3051dc2ce9a2150562069bfe640..c3f0f75c9d1bebf4aff8cea366591877a4a39dea 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -6,7 +6,9 @@
   obs1 <- array(rnorm(72), dim = c(dataset = 1, member = 1, sdate = 3, ftime = 4, lat = 2, lon = 3))
   lat1 <- c(20, 80)
   lon1 <- c(40, 280, 350)
-
+  set.seed(3)
+  exp1_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 3, lat = 2, lon = 3))
+ 
   # dat2
   set.seed(1)
   exp2 <- array(rnorm(216), dim = c(sdate = 3, ftime = 4, member = 2, lat = 3, lon = 3))
@@ -14,7 +16,8 @@
   obs2 <- array(rnorm(108), dim = c(sdate = 3, ftime = 4, lat = 3, lon = 3))
   lat2 <- c(80, 50, 20)
   lon2 <- c(-80, 0, 40)
-
+  set.seed(3)
+  exp2_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 2, lat = 3, lon = 3))
  
 ##############################################
 test_that("1. Input checks", {
@@ -57,7 +60,7 @@ test_that("1. Input checks", {
   )
   expect_error(
     NAO(exp1, obs1, time_dim = 'a'),
-    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+    "Parameter 'time_dim' is not found in 'exp' dimension."
   )
   # memb_dim
   expect_error(
@@ -153,7 +156,7 @@ test_that("1. Input checks", {
   )
   # ncores
   expect_error(
-    NAO(exp1, obs1, lat1, lon1, ncore = 3.5),
+    NAO(exp1, obs1, lat = lat1, lon = lon1, ncore = 3.5),
     "Parameter 'ncores' must be a positive integer."
   )
 
@@ -220,6 +223,38 @@ test_that("2. dat1", {
   c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp)
   )
 
+  # exp_cor
+  expect_equal(
+  names(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)),
+  c("exp", "obs", "exp_cor")
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor),
+  c(sdate = 1, member = 3, dataset = 1)
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp),
+  c(sdate = 3, member = 2, dataset = 1)
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs),
+  c(sdate = 3, member = 1, dataset = 1)
+  )
+  expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor),
+  c(0.3896168, 0.4384543, -0.1302738),
+  tolerance = 0.0001
+  )
+  expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp),
+  c(-0.1688756, -0.2658420, -0.8049575, -0.3022108, -0.3655258, 0.1237722),
+  tolerance = 0.0001
+  )
+    expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs),
+  c(-0.1762489, 0.1364694, -1.4581406),
+  tolerance = 0.0001
+  )
 
 })
 ##############################################
@@ -253,6 +288,39 @@ test_that("3. dat2", {
   tolerance = 0.00001
   )
 
+  # exp_cor
+  expect_equal(
+  names(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)),
+  c("exp", "obs", "exp_cor")
+  )
+  expect_equal(
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor),
+  c(sdate = 1, member = 2)
+  )
+  expect_equal(
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp),
+  c(sdate = 3, member = 2)
+  )
+  expect_equal(
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs),
+  c(sdate = 3)
+  )
+  expect_equal(
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)),
+  c(0.2121340, 0.1634516),
+  tolerance = 0.0001
+  )
+  expect_equal(
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp)),
+  c(0.01457391, 0.06668166, 0.20193275, -0.20154315, -0.49487925, -0.04181974),
+  tolerance = 0.0001
+  ) 
+  expect_equal(
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)),
+  c(0.3511294, -0.7196260, -1.5123894),
+  tolerance = 0.0001
+  )
+
 })
 
 ##############################################