diff --git a/R/NAO.R b/R/NAO.R
index 6d48dba6aaec6eb8fc41dba2462c097e514ccc3f..fb5220cf5ce02e73987ec9e85c951d0234c213b3 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -33,8 +33,8 @@
 #'@param ftime_dim A character string indicating the name of the forecast time 
 #'  dimension of 'exp' and 'obs'. The default value is 'ftime'.
 #'@param ftime_avg A numeric vector of the forecast time steps to average
-#'  across the target period. The default value is 2:4, i.e., from 2nd to 4th 
-#'  forecast time steps.
+#'  across the target period. If average is not needed, set NULL. The default 
+#'  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 
@@ -48,11 +48,13 @@
 #'A list which contains:
 #'\item{exp}{
 #'  A numeric array of forecast NAO index in verification format with the same 
-#'  dimensions as 'exp' except space_dim and ftime_dim.
+#'  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
-#'  dimensions as 'obs' except space_dim and ftime_dim.
+#'  dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL,
+#'  ftime_dim remains.
 #'}
 #'
 #'@references
@@ -198,16 +200,18 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
     }
   }
   ## ftime_avg
-  if (!is.vector(ftime_avg) | !is.integer(ftime_avg)) {
-    stop("Parameter 'ftime_avg' must be an integer vector.")
-  }
-  if (!is.null(exp)) {
-    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.")
+  if (!is.null(ftime_avg)) {
+    if (!is.vector(ftime_avg) | !is.numeric(ftime_avg)) {
+      stop("Parameter 'ftime_avg' must be an integer vector.")
     }
-  } else {
-    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)) {
+      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 (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) {
+        stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
+      }
     }
   }
   ## sdate >= 2
@@ -281,22 +285,24 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
   }
   ## ncores
   if (!is.null(ncores)) {
-    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores == 0 |
       length(ncores) > 1) {
       stop("Parameter 'ncores' must be a positive integer.")
     }
   }
 
   # Average ftime 
-  if (!is.null(exp)) {
-    exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE)
-    exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE)
-    ## Cross-validated PCs. Fabian. This should be extended to
-    ## nmod and nlt by simple loops. Virginie
-  }
-  if (!is.null(obs)) {
-    obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
-    obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
+  if (!is.null(ftime_avg)) {
+    if (!is.null(exp)) {
+      exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE)
+      exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE)
+      ## Cross-validated PCs. Fabian. This should be extended to
+      ## nmod and nlt by simple loops. Virginie
+    }
+    if (!is.null(obs)) {
+      obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
+      obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
+    }
   }
 
   # wght
diff --git a/R/ProjectField.R b/R/ProjectField.R
index 309f3efd731824e3e031a473a22a657e7d7ed845..55e7fd20e39fe8c9ec9b62f59b4e52d1b2f814cc 100644
--- a/R/ProjectField.R
+++ b/R/ProjectField.R
@@ -237,10 +237,14 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon
     # Weight
     e.1 <- eof_mode * wght
     ano <- ano * InsertDim(wght, 1, ntime)
+    #ano <- aaply(ano, 1, '*', wght)  # much heavier
 
-    na <- apply(ano, 1, mean, na.rm = TRUE)  # if [lat, lon] all NA, it's NA
+    na <- rowMeans(ano, na.rm = TRUE) # if [lat, lon] all NA, it's NA
+    #na <- apply(ano, 1, mean, na.rm = TRUE)  # much heavier
     tmp <- ano * InsertDim(e.1, 1, ntime)  # [sdate, lat, lon]
-    pc.ver <- apply(tmp, 1, sum, na.rm = TRUE) 
+    rm(ano)
+    #pc.ver <- apply(tmp, 1, sum, na.rm = TRUE)  # much heavier 
+    pc.ver <- rowSums(tmp, na.rm = TRUE)
     pc.ver[which(is.na(na))] <- NA
 
   } else {  # mode = NULL
@@ -250,7 +254,7 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon
     ano <- ano * InsertDim(wght, 1, ntime)
     dim(ano) <- c(ntime, prod(dim(ano)[2:3]))  # [sdate, lat*lon]
 
-    na <- apply(ano, 1, mean, na.rm = TRUE)  # if [lat, lon] all NA, it's NA
+    na <- rowMeans(ano, na.rm = TRUE) # if [lat, lon] all NA, it's NA
     na <- aperm(array(na, dim = c(ntime, dim(e.1)[1])), c(2, 1))
 
     # Matrix multiplication e.1 [mode, lat*lon] by ano [lat*lon, sdate]
diff --git a/man/NAO.Rd b/man/NAO.Rd
index 64b16562fb45f5ad2a84c6efd5dca8ba8e171876..999fd75f6f45177e353006e6bd061e094f5678c2 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -50,8 +50,8 @@ of 'ano'. The default value is c('lat', 'lon').}
 dimension of 'exp' and 'obs'. The default value is 'ftime'.}
 
 \item{ftime_avg}{A numeric vector of the forecast time steps to average
-across the target period. The default value is 2:4, i.e., from 2nd to 4th 
-forecast time steps.}
+across the target period. If average is not needed, set NULL. The default 
+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 
@@ -67,11 +67,13 @@ computation. The default value is NULL.}
 A list which contains:
 \item{exp}{
  A numeric array of forecast NAO index in verification format with the same 
- dimensions as 'exp' except space_dim and ftime_dim.
+ 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
- dimensions as 'obs' except space_dim and ftime_dim.
+ dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL,
+ ftime_dim remains.
 }
 }
 \description{
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index 2da0d407ea986175a9a3d2c73605bd5aead29451..05fcd22caec34d412f254f64e9819484a0269e04 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -16,6 +16,7 @@ context("s2dv::NAO tests")
   obs2 <- array(rnorm(108), dim = c(sdate = 3, ftime = 4, lat = 3, lon = 3))
   lat2 <- c(80, 50, 20)
   lon2 <- c(-80, 0, 40)
+
  
 ##############################################
 test_that("1. Input checks", {
@@ -208,6 +209,19 @@ test_that("2. dat1", {
   c(-0.1139683, 0.1056687, 0.1889449),
   tolerance = 0.0001
   )
+  expect_equal(
+  dim(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp),
+  c(sdate = 3, member = 2, dataset = 1)
+  )
+  expect_equal(
+  dim(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = NULL)$exp),
+  c(sdate = 3, member = 2, dataset = 1, ftime = 4)
+  )
+  expect_equal(
+  c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = NULL)$exp[,,1,1]),
+  c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp)
+  )
+
 
 })
 ##############################################