diff --git a/.Rbuildignore b/.Rbuildignore
index 91114bf..dc93b90 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,5 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^data-raw$
+^_config\.yml$
+^cran-comments\.md$
diff --git a/DESCRIPTION b/DESCRIPTION
index 986e4e8..3ea18ac 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,6 @@
 Package: SPEI
 Type: Package
-Version: 1.8
-Date: 2022-05-05
+Version: 1.8.0
 Title: Calculation of the Standardised Precipitation-Evapotranspiration Index
 Authors@R: c(
 	person('Santiago', 'Beguería', role=c('aut','cre'),
@@ -12,8 +11,9 @@ Description: A set of functions for computing potential evapotranspiration and s
 Depends: R (>= 3.5.0)
 Imports: lmomco, lmom, TLMoments, reshape, ggplot2, checkmate, zoo, lubridate
 License: GPL-2
-URL: http://spei.csic.es
-LazyLoad: yes
+URL: http://spei.csic.es, https://github.com/sbegueria/SPEI
+BugReports: https://github.com/sbegueria/SPEI/issues
+LazyLoad: true
 Encoding: UTF-8
 Suggests: covr, testthat
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.0
diff --git a/NAMESPACE b/NAMESPACE
index 7544995..2c99f8a 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -22,9 +22,7 @@ importFrom(graphics,par)
 importFrom(graphics,plot)
 importFrom(graphics,points)
 importFrom(graphics,polygon)
-importFrom(lmom,cdfgam)
 importFrom(lmom,cdfglo)
-importFrom(lmom,cdfpe3)
 importFrom(lmom,pelgam)
 importFrom(lmom,pelglo)
 importFrom(lmom,pelpe3)
@@ -40,12 +38,20 @@ importFrom(lmomco,pwm2lmom)
 importFrom(lubridate,days_in_month)
 importFrom(lubridate,yday)
 importFrom(reshape,melt)
+importFrom(stats,aggregate)
 importFrom(stats,cycle)
 importFrom(stats,end)
 importFrom(stats,frequency)
+importFrom(stats,is.ts)
 importFrom(stats,optim)
+importFrom(stats,pnorm)
+importFrom(stats,qnorm)
+importFrom(stats,sd)
 importFrom(stats,start)
+importFrom(stats,time)
 importFrom(stats,ts)
+importFrom(stats,window)
+importFrom(zoo,as.Date.yearmon)
 importFrom(zoo,as.yearmon)
 importFrom(zoo,na.trim)
 importFrom(zoo,rollapply)
diff --git a/NEWS.md b/NEWS.md
index 9066e98..ac176e5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,7 +6,7 @@
 
 __Version history:
 
-_Version 1.8, May 2022 (current development version on github).
+_### _Version 1.8.0, November 2022 (current version on GitHub, submitted to CRAN).
 
 1. Solving several minor bugs in `<thornthwaite>`, `<hargreaves>`, and `<penman>` functions (output difference less than 0.1%).
 2. Solving a bug in `<spei>` that resulted in bad cumulative data when using a non-rectangular kernel, resulting in incorrect SPEI values.
@@ -16,16 +16,16 @@ _Version 1.8, May 2022 (current development version on github).
 6. Implementation of different versions of the Penman-Monteith ETo calculation in function `<penman>`.
 7. Implementation of an option to include CO2 concentration data in function `<penman>`.
 8. Implementation of a new option for when no wind data are available in function `<penman>`.
-9. Funtion `<plot.spei>` completely rewritten based on `ggplot2`, solving some bugs and enabling more flexibility.
+9. Function `<plot.spei>` completely rewritten based on `ggplot2`, solving some bugs and enabling more flexibility.
 
-_Version 1.7.2, June 2019 (current stable version on github).
+_Version 1.7.2, June 2019 (only on GitHub).
 
 1. Allowing for data with frequency other than 12 in `<spei>` function.
 
 _Version 1.7.1, June 2017.
 
-1. Corrected an error in `<spei>` function, which was not working when distribution was Gamma or PeasonIII and using user provided parameters. (Fixed by Emanuele Cordano, emanuele.cordano@gmail.com -- ecor)
-2. Added probability of monthly precipitation = 0 (pze) when using user provided parameters. (Fixed by Emanuele Cordano, emanuele.cordano@gmail.com -- ecor) 
+1. Corrected an error in `<spei>` function, which was not working when distribution was Gamma or PeasonIII and using user provided parameters. (Fixed by Emanuele Cordano, emanuele.cordano@gmail.com -- 'ecor')
+2. Added probability of monthly precipitation = 0 (pze) when using user provided parameters. (Fixed by Emanuele Cordano, emanuele.cordano@gmail.com -- 'ecor') 
 
 _Version 1.7, June 2017 (current on CRAN).
 
@@ -45,7 +45,7 @@ _Version 1.5, May 2013.
 
 _Version 1.4, May 2013.
 1. Minor fixes to functions \code{\link{penman}} and \code{\link{pwm}}.
-2. Documentation of the penman function defined by mistake ed as the saturation vapour pressure, while it should read 'actual vapour pressure'.
+2. Documentation of the penman function defined by mistake ed as the saturation vapor pressure, while it should read 'actual vapor pressure'.
 3. Function zzz.R added to display basic information about the SPEI package at startup.
 4. Function \code{\link{SPEINews}} added to display the NEWS file.
 
@@ -56,7 +56,7 @@ _Version 1.3, March 2013.
 
 _Version 1.2, October 2012.
 1. Fixed a bug causing several functions to fail when a time series not belonging to matrix class was provided.
-2. Function \code{\link{plot.spei}} now distinguises between calls to spei and spi and labels the axis accordingly.
+2. Function \code{\link{plot.spei}} now distinguishes between calls to \code{\link{spei}} and \code{\link{spi}} and labels the axis accordingly.
 
 _Version 1.1, March 2012.
 1. Functions \code{\link{spei}} and \code{\link{spi}} now yield an object of class "spei".
@@ -69,6 +69,6 @@ First release of the SPEI package.
 __To do (work in progress):
 
 1. Complete documentation for pwmLC.Rd.
-2. Review method plot.spei() that produces wrong results in some cases.
+2. Review method \code{\link{plot.speispei}} that produces wrong results in some cases.
 3. Implement parallel processing.
 4. Analysis functions.
diff --git a/R/SPEI-package.r b/R/SPEI-package.r
index ae7298a..ac6aed7 100755
--- a/R/SPEI-package.r
+++ b/R/SPEI-package.r
@@ -43,5 +43,10 @@
 #' 
 #' @importFrom checkmate makeAssertCollection
 #' @importFrom stats cycle end frequency start ts optim
-#' @importFrom zoo as.yearmon rollapply
+#' @importFrom zoo as.yearmon as.Date.yearmon rollapply na.trim
+#' @importFrom stats aggregate is.ts pnorm qnorm sd time window
+#' @importFrom TLMoments PWM
+#' @importFrom lmomco are.lmom.valid are.parglo.valid cdfgam cdfpe3 cdfgam pargam parglo parpe3 pwm.pp pwm2lmom
+#' @importFrom lmom pelgam pelglo pelpe3 cdfglo
+
 NULL
diff --git a/R/data.R b/R/data.R
index 9cc1eb4..e4fc1aa 100755
--- a/R/data.R
+++ b/R/data.R
@@ -1,12 +1,7 @@
 #' @name Datasets
-#' 
-#' 
-#' @aliases balance cabinda cruts4
-#' 
-#' 
+#' @aliases wichita balance cabinda cruts4
 #' @title Data sets for illustrating the functions in the SPEI package.
-#'
-#' 
+#' @keywords datasets
 #' @description 
 #' Data used in the examples of the SPEI package: 
 #' \code{wichita} dataset: monthly climate in Wichita (Kansas, lat=37.6475, 
@@ -18,19 +13,7 @@
 #' Allen et al. (1998);
 #' \code{cruts4}: 120 years of monthly climatic water balance (precipitation 
 #' minus reference evapotranspiration) data at six grid points from CRU TS 4.05. 
-#' 
-#' 
-#'
-#'  
 #' @details See description.
-#' 
-#'
-#' @usage 
-#' data(wichita)
-#' data(balance)
-#' data(cabinda)
-#'
-#'
 #' @format 
 #' \code{wichita} dataset: a data frame with:
 #' \describe{
@@ -44,13 +27,11 @@
 #'   \item{ACSH}{ monthly mean sun hours, in h.}
 #'   \item{ACSH}{ monthly mean cloud cover, in \%.}
 #'}
-#'
 #' \code{balance} dataset: a data frame with monthly climatic water balance
 #'  (precipitation minus potential evapotranspiration) at Indore (India),
 #'  Kimberley (South Africa), Albuquerque (US), Valencia (Spain), Wien (Austria),
-#'  Abashiri (Japan), Tampa (US), Sao Paulo (Brasil), Lahore (India), Punta
+#'  Abashiri (Japan), Tampa (US), Sao Paulo (Brazil), Lahore (India), Punta
 #'  Arenas (Chile) and Helsinki (Finland), in mm.
-#' 
 #' \code{cabinda} dataset: a data frame with one year of monthly climatic data
 #' at Cabinda (Angola, -5.33S 12.11E 20 m), with:
 #'  
@@ -61,7 +42,7 @@
 #'   \item{RH}{ monthly mean relative humidity, in \%.}
 #'   \item{U2}{ monthly mean wind speed, in km h-1}
 #'   \item{tsun}{ monthly mean sunshine hours, in h.}
-#'   \item{Rs}{ monthly mean dialy incoming solar radiation, MJ m-2 d-1.}
+#'   \item{Rs}{ monthly mean diaily incoming solar radiation, MJ m-2 d-1.}
 #'   \item{ET0}{ monthly ET0 from the original publication, in mm.}
 #'}
 #'
@@ -101,50 +82,19 @@
 #' summary(wichita)
 #' data(balance)
 #' summary(balance)
+#' data(cruts4)
+#' summary(cruts4)
 #'
-#'
-#'
-"wichita"
+NULL
 
+#' @rdname wichita
+"wichita"
 
-#' 
-#' @title cabinda
-#' 
-#' 
-#' @rdname Datasets
-#' 
-#' 
-#' @description See wichita
-#' 
-#' 
-#' @details See wichita
-#' 
+#' @rdname cabinda
 "cabinda"
 
-
-#' 
-#' @title balance
-#' 
-#' 
-#' @rdname Datasets
-#' 
-#' 
-#' @description See wichita
-#' 
-#' 
-#' @details See wichita
-#' 
+#' @rdname balance
 "balance"
 
-#' 
-#' @title cruts4
-#' 
-#' 
-#' @rdname Datasets
-#' 
-#' 
-#' @description See wichita
-#' 
-#' 
-#' @details See wichita
+#' @rdname cruts4
 "cruts4"
\ No newline at end of file
diff --git a/R/hargreaves.R b/R/hargreaves.R
index 75afaeb..044a78a 100755
--- a/R/hargreaves.R
+++ b/R/hargreaves.R
@@ -1,21 +1,15 @@
 #' @name Potential-evapotranspiration
-#' 
-#' 
 #' @title Computation of potential and reference evapotranspiration.
-#' 
-#' 
 #' @aliases thornthwaite penman
-#' 
-#' 
 #' @usage 
 #' thornthwaite(Tave, lat, na.rm = FALSE, verbose=TRUE)
 #' 
-#' hargreaves(Tmin, Tmax, Ra = NA, lat = NA, Pre = NA, na.rm = FALSE, verbose=TRUE)
+#' hargreaves(Tmin, Tmax, Ra = NULL, lat = NULL, Pre = NULL, na.rm = FALSE, verbose=TRUE)
 #' 
-#' penman(Tmin, Tmax, U2, Ra = NA, lat = NA, Rs = NA, tsun = NA,
-#'        CC = NA, ed = NA, Tdew = NA, RH = NA, P = NA, P0 = NA,
-#'               z = NA, crop='short', na.rm = FALSE, method='ICID',
-#'               verbose=TRUE)
+#' penman(Tmin, Tmax, U2, Ra = NULL, lat = NULL, Rs = NULL, tsun = NULL,
+#'        CC = NULL, ed = NULL, Tdew = NULL, RH = NULL, P = NULL, P0 = NULL,
+#'        CO2 = NULL, z = NULL, crop='short', na.rm = FALSE, method='ICID',
+#'        verbose=TRUE)
 #' 
 #' 
 #' @description 
@@ -35,10 +29,10 @@
 #' @param Ra   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily external radiation, MJ m-2 d-1.
 #' @param Pre   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array  of monthly total precipitation, mm.
 #' @param U2   a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily wind speeds at 2 m height, m s-1.
-#' @param Rs   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean dialy incoming solar radiation, MJ m-2 d-1.
+#' @param Rs   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily incoming solar radiation, MJ m-2 d-1.
 #' @param tsun   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily bright sunshine hours, h.
 #' @param CC   optional, numeric a vector, matrix or time series of monthly mean cloud cover, \%.
-#' @param ed   optional, numeric a vector, matrix or time series of monthly mean actual vapour pressure at 2 m height, kPa.
+#' @param ed   optional, numeric a vector, matrix or time series of monthly mean actual vapor pressure at 2 m height, kPa.
 #' @param Tdew   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily dewpoint temperature (used for estimating ed), ºC.
 #' @param RH   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean relative humidity (used for estimating ed), \%.
 #' @param P   optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean atmospheric pressure at surface, kPa.
@@ -109,7 +103,7 @@
 #'
 #' @references 
 #' Thornthwaite, C. W., 1948. An approach toward a rational classification of climate. 
-#' \emph{Geographical Review} \bold{38}: 55–94. doi:10.2307/2107309.
+#' \emph{Geographical Review} \bold{38}: 55–94. DOI:10.2307/2107309.
 #' 
 #' Hargreaves G.H., 1994. Defining and using reference evapotranspiration. 
 #' \emph{Journal of Irrigation and Drainage Engineering} \bold{120}: 1132–1139.
@@ -131,7 +125,7 @@
 #' Reston, VA, 57 pp.
 #' 
 #' Yang, Y., Roderick, M.L., Zhang, S. McVicar, T., Donohue, R.J., 2019. Hydrologic implications of vegetation
-#' response to elevated CO2 in climate projections. \emph{Nature Clim Change} \bold{9}: 44–48.
+#' response to elevated CO2 in climate projections. \emph{Nature Climate Change} \bold{9}: 44–48.
 #'
 #'
 #' @author Santiago Beguería
@@ -308,6 +302,7 @@ hargreaves <- function(Tmin, Tmax, Ra=NULL, lat=NULL, Pre=NULL,
     # 3D array input (gridded data)
     int_dims <- tmin_dims
   } else {
+    int_dims <- tmin_dims
     check$push('Input data can not have more than three dimensions.')
   }
   n_sites <- prod(int_dims[[2]], int_dims[[3]])
@@ -341,7 +336,7 @@ hargreaves <- function(Tmin, Tmax, Ra=NULL, lat=NULL, Pre=NULL,
     }
     ym <- as.yearmon(time(Tmin))
     warn$push(paste0('Time series spanning ', ym[1], ' to ', ym[n_times], '.'))
-    date <- as.Date(ym)
+    date <- as.Date.yearmon(ym)
     mlen_array <- array(as.numeric(lubridate::days_in_month(date)), dim=int_dims)
     msum_array <- array(yday(date) + round((mlen_array/2) - 1), dim=int_dims)
   } else {
@@ -355,7 +350,7 @@ hargreaves <- function(Tmin, Tmax, Ra=NULL, lat=NULL, Pre=NULL,
   # Verify the length of each input variable
   input_len <- prod(int_dims)
   if (sum(lengths(Tmin))!=input_len || sum(lengths(Tmax))!=input_len) {
-    check$push('`Tmin` and `Tmax` should not have different lengths.')
+    check$push('`Tmin` and `Tmax` cannot have different lengths.')
   }
   if (using$Ra && sum(lengths(Ra))!=input_len) {
     check$push('`Ra` has incorrect length.')
@@ -462,4 +457,4 @@ hargreaves <- function(Tmin, Tmax, Ra=NULL, lat=NULL, Pre=NULL,
   dimnames(ET0) <- names
   
   return(ET0)
-}
\ No newline at end of file
+}
diff --git a/R/kern.R b/R/kern.R
index e2ec2f7..d83856b 100755
--- a/R/kern.R
+++ b/R/kern.R
@@ -1,25 +1,15 @@
 #' @name Kernel-functions
-#' 
-#' 
 #' @title Time kernel for computing the SPEI at different time scales.
-#'
-#'
 #' @aliases kern.plot
-#'
-#'
 #' @description Function \code{kern} is used internally by \code{\link{spei}} and
 #'  \code{\link{spi}}for computing drought indices at different time scales.
-#' 
-#' 
 #' @param scale numeric, time scale or length of the kernel.
 #' @param type character, shape of the kernel function.
 #' @param shift numeric, shifting of the kernel peak.
-#' 
-#' 
 #' @details 
 #' Drought indices, such as the SPEI or the SPI, are usually 
 #' computed at different time scales to adapt to the different response 
-#' times of systems affected by drought. This is acomplished by applying 
+#' times of systems affected by drought. This is accomplished by applying 
 #' a kernel function to the data prior to computation of the SPEI. 
 #' Application of a kernel has the effect of smoothing the temporal 
 #' variability of the resulting SPEI, allowing for the major patterns 
@@ -50,7 +40,7 @@
 #'  time lag for the four different kernel shapes so they can be compared.
 #'  
 #'  
-#' @return A vector of lenght equal to \code{scale} with weights used for computing the drought index.
+#' @return A vector of length equal to \code{scale} with weights used for computing the drought index.
 #' 
 #' 
 #' @references 
@@ -96,7 +86,7 @@ kern <- function(scale, type='rectangular', shift=0) {
 	
 	if(s == 1) type == "rectangular"
 	
-	k = switch(type,
+	k <- switch(type,
 	       rectangular = rep(1,s),
 	       triangular = s:1,
 	       circular = (s^2+(1-(1:s)^2)),
diff --git a/R/kern.plot.R b/R/kern.plot.R
index 9ae5e58..245fdbb 100755
--- a/R/kern.plot.R
+++ b/R/kern.plot.R
@@ -1,18 +1,8 @@
 #' @title Plot the kern object
-#' 
-#' 
 #' @description See kern
-#' 
-#' 
 #' @details See kern
-#'
-#'
 #' @rdname Kernel-functions
-#'  
-#'  
 #' @importFrom graphics plot par
-#'   
-#'   
 #' @export
 #' 
 kern.plot  <- function(scale=12, shift=0) {
diff --git a/R/parglo.maxlik.R b/R/parglo.maxlik.R
index 2ba0752..5461cdd 100755
--- a/R/parglo.maxlik.R
+++ b/R/parglo.maxlik.R
@@ -1,45 +1,31 @@
 #' @name Generalized-Logistic
-#' 
-#' 
 #' @title Generalized Logistic maximum likelihood function
-#' 
-#' 
 #' @description Maximum likelihood fitting function for 
 #' generalized logistic distribution.
-#' 
-#' 
 #' @details This function is used internally by \code{spei} 
 #' and \code{spi} and is supposed to never be needed by the 
 #' regular user. Initial values for maximum likelihood estimation
 #' can be provided by \code{parglo}.
-#' 
-#' 
 #' @usage 
 #' parglo.maxlik(x, ini)
 #'
-#'
 #' @param x vector of quantiles for which to evaluate the PDF.
 #' @param ini a vector of initial values of the parameters to be fit.
 #'
-#'
 #' @return a list of parameters of a generalized Logistic 
 #' distribution function
 #' 
-#' 
 #' @references 
 #' S.M. Vicente-Serrano, S. Beguería, J.I. López-Moreno. 2010. 
 #' A Multi-scalar drought index sensitive to global warming: 
 #' The Standardized Precipitation Evapotranspiration Index – SPEI. 
 #' \emph{Journal of Climate} \bold{23}: 1696, DOI: 10.1175/2009JCLI2909.1.
 #'
-#'
 #' @author Santiago Beguería
 #'
-#'
 #' @importFrom stats optim
 #' @importFrom lmomco are.parglo.valid
 #'
-#'
 #' @export
 #'
 parglo.maxlik <- function(x,ini) {
diff --git a/R/penman.R b/R/penman.R
index 704214e..553029f 100755
--- a/R/penman.R
+++ b/R/penman.R
@@ -1,12 +1,6 @@
 #' @title Computation of potential evapotranspiration.
-#' 
-#' 
 #' @description See hargreaves
-#' 
-#' 
 #' @details See hargreaves
-#' 
-#' 
 #' @return  A time series with the values of monthly potential or reference 
 #' evapotranspiration, in mm. 
 #' If the input is a matrix or a multivariate time series each column will be 
@@ -157,6 +151,7 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
     # 3D array input (gridded data)
     int_dims <- tmin_dims
   } else {
+    int_dims <- tmin_dims
     check$push('Input data can not have more than 3 dimensions')
   }
   n_sites <- prod(int_dims[[2]], int_dims[[3]])
@@ -190,7 +185,7 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
     }
     ym <- as.yearmon(time(Tmin))
     warn$push(paste0('Time series spanning ', ym[1], ' to ', ym[n_times], '.'))
-    date <- as.Date(ym)
+    date <- as.Date.yearmon(ym)
     mlen_array <- array(as.numeric(days_in_month(date)), dim=int_dims)
     msum_array <- array(yday(date) + round((mlen_array/2) - 1), dim=int_dims)
   } else {
@@ -204,7 +199,7 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
   # Verify the length of each input variable
   input_len <- prod(int_dims)
   if (sum(lengths(Tmin))!=input_len || sum(lengths(Tmax))!=input_len) {
-    check$push('`Tmin` and `Tmax`cannot have different lengths.')
+    check$push('`Tmin` and `Tmax` cannot have different lengths.')
   }
   if (using$U2 && sum(lengths(U2))!=input_len) {
     check$push('`U2` has incorrect length.')
@@ -230,7 +225,7 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
   if (using$Tdew && sum(lengths(Tdew))!=input_len) {
     check$push('`Tdew` has incorrect length.')
   }
-  if (using$RH && sum(lengths(RRHa))!=input_len) {
+  if (using$RH && sum(lengths(RH))!=input_len) {
     check$push('`RH` has incorrect length.')
   }
   if (using$P && sum(lengths(P))!=input_len) {
@@ -470,7 +465,7 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
   # Daily ET0 (eq. 2.18)
   if (crop=='short') {
     c1 <- 900; c2 <- 0.34 # short reference crop (e.g. clipped grass, 0.12 m)
-  } else if (crop=='long') {
+  } else if (crop=='tall') {
     c1 <- 1600; c2 <- 0.38 # tall reference crop (e.g. alfalfa, 0.5 m)
   } else {
     stop(paste('An error occurred while estimating the daily ET0',
@@ -509,4 +504,4 @@ penman <- function(Tmin, Tmax, U2=NULL, Ra=NULL, lat=NULL, Rs=NULL,
   dimnames(ET0) <- names
   
   return(ET0)
-}
\ No newline at end of file
+}
diff --git a/R/spei.R b/R/spei.R
index 3f611e1..d89e567 100644
--- a/R/spei.R
+++ b/R/spei.R
@@ -1,28 +1,41 @@
 spei <- function(x, y,...) UseMethod('spei')
 
-
 #' @name Drought-indices
-#' 
-#' @aliases spi
-#' 
 #' @title Calculation of the Standardized Precipitation-Evapotranspiration 
 #' Index (SPEI) and the Standardized Precipitation Index (SPI).
-#' 
-#' 
+#' @aliases spei, spi
 #' @description
 #' Given a time series of the climatic water balance (precipitation minus 
 #' potential evapotranspiration), gives a time series of the Standardized 
 #' Precipitation-Evapotranspiration Index (SPEI).
-#' 
-#' 
-#' @usage 
-#' spei(data, scale, kernel = list(type = 'rectangular', shift = 0),
-#' distribution = 'log-Logistic', fit = 'ub-pwm', na.rm = FALSE,
-#' ref.start=NULL, ref.end=NULL, keep.x=FALSE, params=NULL, 
-#' verbose=TRUE, ...)
-#'
-#' 
-#' 
+#' @usage
+#' spei(
+#'  data,
+#'  scale,
+#'  kernel = list(type = 'rectangular', shift = 0),
+#'  distribution = 'log-Logistic',
+#'  fit = 'ub-pwm',
+#'  na.rm = FALSE,
+#'  ref.start=NULL,
+#'  ref.end=NULL,
+#'  keep.x=FALSE,
+#'  params=NULL, 
+#'  verbose=TRUE,
+#'  ...)
+#' 
+#' spi(
+#'  data,
+#'  scale,
+#'  kernel = list(type = 'rectangular', shift = 0),
+#'  distribution = 'Gamma',
+#'  fit = 'ub-pwm',
+#'  na.rm = FALSE,
+#'  ref.start=NULL,
+#'  ref.end=NULL,
+#'  keep.x=FALSE,
+#'  params=NULL, 
+#'  verbose=TRUE,
+#'  ...)
 #' @param data a vector, matrix or data frame with time ordered values 
 #' of precipitation (for the SPI) or of the climatic balance 
 #' precipitation minus potential evapotranspiration (for the SPEI).
@@ -100,7 +113,7 @@ spei <- function(x, y,...) UseMethod('spei')
 #' 
 #' @section Probability distributions:
 #' Following the original definitions \code{spei} uses a log-Logistic distribution 
-#' by default, and \code{spi} uses a Gamma distribution. This behaviour can be modified, 
+#' by default, and \code{spi} uses a Gamma distribution. This behavior can be modified, 
 #' however, through parameter \code{distribution}.
 #' 
 #' 
@@ -125,7 +138,7 @@ spei <- function(x, y,...) UseMethod('spei')
 #' 
 #' 
 #' @section Reference period:
-#' The default behaviour of the functions is using all the values provided in \code{data} 
+#' The default behavior of the functions is using all the values provided in \code{data} 
 #' for parameter fitting. However, this can be modified with help of parameters \code{ref.start} 
 #' and \code{ref.end}. These parameters allow defining a subset of values that will be used 
 #' for parameter fitting, i.e. a reference period. The functions, however, will compute the 
@@ -281,6 +294,7 @@ spei <- function(x, y,...) UseMethod('spei')
 #' # Modding the plot
 #' # Since plot.spei() returns a ggplot object, it is possible to add or tweak
 #' # parts of the plot.
+#' require(ggplot2)
 #' plot(spei(wichita[,'BAL'], 12)) +
 #'  ggtitle('SPEI1 at Wichita') +
 #'  scale_fill_manual(values=c('blue','red')) +  # classic SPEI look
@@ -288,11 +302,6 @@ spei <- function(x, y,...) UseMethod('spei')
 #'  theme_classic() +
 #'  theme(legend.position='bottom')
 #' 
-#' @importFrom zoo rollapply
-#' @importFrom TLMoments PWM
-#' @importFrom lmomco are.lmom.valid are.parglo.valid cdfgam cdfpe3 pargam parglo parpe3 pwm.pp pwm2lmom
-#' @importFrom lmom pelgam pelglo pelpe3 cdfglo cdfgam cdfpe3
-#' 
 #' @export
 #' 
 spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
@@ -427,6 +436,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
     # 3D array input (gridded data)
     int_dims <- data_dims
   } else {
+    int_dims <- data_dims
     check$push('Input data can not have more than three dimensions.')
   }
   n_sites <- prod(int_dims[[2]], int_dims[[3]])
@@ -492,7 +502,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
   
   # Instantiate an object to store the distribution coefficients
   # ADD PZE TO GAMMA AND PEARSONIII
-  coef = switch(distribution,
+  coef <- switch(distribution,
                 "Gamma" = array(NA, c(2, n_sites, ts_freq),
                                 list(par=c('alpha','beta'), colnames(data), NULL)),
                 "PearsonIII" = coef <- array(NA, c(3, n_sites, ts_freq),
@@ -523,7 +533,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
   
   # Convert to time series
   if (!is.ts(acu)) {
-    acu <- ts(acu, start=ts_start, fr=ts_freq)
+    acu <- ts(acu, start=ts_start, frequency=ts_freq)
   }
   
   # Trim data set to reference period for fitting (acu.ref)
@@ -561,7 +571,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
       # Probability of zero (pze)
       if(distribution != 'log-Logistic'){
         pze <- sum(x.mon==0) / length(x.mon)
-        x.mon = x.mon[x.mon > 0]
+        x.mon <- x.mon[x.mon > 0]
       }
       
       ## Compute coefficients - - - - - - - - - - - - - -
@@ -569,7 +579,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
       # Distribution parameters (f_params)
       if (!using$params) {
         # Fit distribution parameters
-        x.mon_sd = sd(x.mon, na.rm=TRUE)
+        x.mon_sd <- sd(x.mon, na.rm=TRUE)
         
         # Early stopping
         if (is.na(x.mon_sd) || (x.mon_sd == 0)) {
@@ -577,13 +587,13 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
           next()
         }
         if(length(x.mon) < 4){
-          spei[ff,s] = NA
+          spei[ff,s] <- NA
           coef[,s,c] <- NA
           next()
         }
         
         # Calculate probability weighted moments based on `lmomco` or `TLMoments`
-        pwm = switch(fit,
+        pwm <- switch(fit,
                      'pp-pwm' = pwm.pp(x.mon, -0.35, 0, nmom=3, sort=TRUE),
                      'ub-pwm' = PWM(x.mon, order=0:2)
         )
@@ -597,10 +607,10 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
         
         # `lmom` fortran functions need specific inputs L1, L2, T3
         # This is handled internally by `lmomco` with `lmorph`
-        fortran_vec = c(lmom$lambdas[1:2], lmom$ratios[3])
+        fortran_vec <- c(lmom$lambdas[1:2], lmom$ratios[3])
         
         # Calculate parameters based on distribution with `lmom`, then `lmomco`
-        f_params = switch(distribution,
+        f_params <- switch(distribution,
                           'log-Logistic' = tryCatch(pelglo(fortran_vec),
                                                     error = function(e){ parglo(lmom)$para }),
                           'Gamma' = tryCatch(pelgam(fortran_vec),
@@ -611,12 +621,12 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
         
         # Adjust if user chose `log-Logistic` and `max-lik`
         if(distribution == 'log-Logistic' && fit=='max-lik'){
-          f_params = parglo.maxlik(x.mon, f_params)$para
+          f_params <- parglo.maxlik(x.mon, f_params)$para
         }
       } else {
         # User-provided distribution parameters
         
-        f_params = as.vector(params[,s,c])
+        f_params <- as.vector(params[,s,c])
         
       }
       
@@ -679,7 +689,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
 #' 
 #' @title Generic methods for \code{spei} objects.
 #' 
-#' @aliases print.spi plot.spei plot.spi summary.spei summary.spi
+#' @aliases plot.spei summary.spei
 #' 
 #' @description 
 #' Generic methods for extracting information and plotting \code{spei} objects.
@@ -687,11 +697,10 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
 #' @usage
 #' \method{print}{spei}(x, ...)
 #' \method{summary}{spei}(object, ...)
-#' \method{plot}{spei}(x, ttext, ...)
+#' \method{plot}{spei}(x, ...)
 #' 
 #' @param x an object of class \code{spei}.
 #' @param object an object of class \code{spei}.
-#' @param ttext text to use as part of the plot title
 #' @param ... additional parameters, not used at present.
 #' 
 #' 
@@ -712,7 +721,7 @@ spei <- function(data, scale, kernel=list(type='rectangular', shift=0),
 #' @author Santiago Beguería
 #'  
 #' @examples 
-#' See examples of use in the help page of the \code{spei()} function.
+#' # See examples of use in the help page of the spei() function.
 #' 
 #' @export
 #' 
@@ -722,17 +731,9 @@ print.spei <- function (x, ...) {
 
 #' 
 #' @title summary of spei/spi
-#' 
-#' 
 #' @description See print.spei
-#' 
-#' 
 #' @details See print.spei
-#' 
-#' 
 #' @rdname Generic-methods-for-spei-objects
-#' 
-#' 
 #' @export
 #' 
 summary.spei <- function (object, ...) {
@@ -752,17 +753,10 @@ summary.spei <- function (object, ...) {
 
 #' 
 #' @title plot spei/spi
-#' 
-#' 
 #' @description See print.spei
-#' 
-#' 
 #' @details See print.spei
-#' 
-#' 
 #' @rdname Generic-methods-for-spei-objects
 #' 
-#' 
 #' @import ggplot2
 #' @importFrom zoo na.trim
 #' @importFrom reshape melt
@@ -770,7 +764,7 @@ summary.spei <- function (object, ...) {
 #' 
 #' @export
 #' 
-plot.spei <- function (x) {
+plot.spei <- function (x, ...) {
   
   ### Argument check - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   
@@ -781,6 +775,9 @@ plot.spei <- function (x) {
   
   ### Make the plot - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   
+  # A workaround to avoid R CMD check warning about no visible binding for global variables
+  #utils::globalVariables(na, value)
+  
   # Label
   if (grepl('spei', x$call[1])) {
     label <- 'SPEI (z-values)'
@@ -834,7 +831,7 @@ plot.spei <- function (x) {
   
   # Convert to time series
   if (!is.ts(data)) {
-    data <- ts(data, start=ts_start, fr=ts_freq)
+    data <- ts(data, start=ts_start, frequency=ts_freq)
   }
   
   # Determine reference period
@@ -883,7 +880,7 @@ plot.spei <- function (x) {
   # kk$cat[w] <- '(-0.5,0]'
   
   # Plot it
-  g <- ggplot(kk, aes(time, value, fill=cat, color=cat))
+  g <- ggplot(kk, aes(.data[['time']], .data[['value']], fill='cat', color='cat'))
   # reference period (if different than whole series)
   if (!is.null(x$ref.period)) {
     g <- g +
@@ -899,7 +896,7 @@ plot.spei <- function (x) {
     scale_color_manual(values=c('cyan3','tomato')) # new look
   # add NAs
   g <- g + 
-    geom_point(aes(time, na), shape=21, fill='white', color='black')
+    geom_point(aes(.data[['time']], .data[['na']]), shape=21, fill='white', color='black')
   # add other parts and options
   g <- g +
     geom_hline(yintercept=0, color='grey') +
diff --git a/R/spi.R b/R/spi.R
index eaeaafd..78e102b 100755
--- a/R/spi.R
+++ b/R/spi.R
@@ -7,9 +7,12 @@ spi <- function(x, y,...) UseMethod('spi')
 #' 
 spi <- function(data, scale, kernel=list(type='rectangular',shift=0),
 	distribution='Gamma', fit='ub-pwm', na.rm=FALSE,
-	ref.start=NULL, ref.end=NULL, x=FALSE, params=NULL, verbose=TRUE, ...){
+	ref.start=NULL, ref.end=NULL, keep.x=FALSE, params=NULL, verbose=TRUE, ...){
+  
+  # input checks
+  
 	sol <- spei(data, scale, kernel, distribution, fit, na.rm,
-		ref.start, ref.end, x, params, verbose)
+		ref.start, ref.end, keep.x, params, verbose)
 	sol$call <- match.call(expand.dots=FALSE)
 
 	return(sol)
diff --git a/R/thornthwaite.R b/R/thornthwaite.R
index d73ad3f..73dd458 100755
--- a/R/thornthwaite.R
+++ b/R/thornthwaite.R
@@ -1,20 +1,11 @@
 #' @title Computation of potential evapotranspiration.
-#' 
-#' 
 #' @description See hargreaves
-#' 
-#' 
 #' @details See hargreaves
-#' 
-#' 
 #' @return A time series with the values of monthly potential or reference evapotranspiration, in mm. 
 #' If the input is a matrix or a multivariate time series each column will be treated as independent 
 #' data (e.g., diferent observatories), and the output will be a multivariate time series.
 #' 
-#' 
 #' @rdname Potential-evapotranspiration
-#' 
-#' 
 #' @export
 #' 
 thornthwaite <- function(Tave, lat, na.rm=FALSE, verbose=TRUE) {
@@ -61,7 +52,8 @@ thornthwaite <- function(Tave, lat, na.rm=FALSE, verbose=TRUE) {
     # 3D array input (gridded data)
     int_dims <- tmin_dims
   } else {
-    check$push('Input data can not have more than 3 dimensions')
+    int_dims <- tmin_dims
+    check$push('Input data can not have more than three dimensions')
   }
   n_sites <- prod(int_dims[[2]], int_dims[[3]])
   n_times <- int_dims[[1]]
@@ -95,13 +87,13 @@ thornthwaite <- function(Tave, lat, na.rm=FALSE, verbose=TRUE) {
     }
     ym <- as.yearmon(time(Tave))
     warn$push(paste0('Time series spanning ', ym[1], ' to ', ym[n_times], '.'))
-    date <- as.Date(ym)
+    date <- as.Date.yearmon(ym)
     mlen_array <- array(as.numeric(lubridate::days_in_month(date)), dim=int_dims)
     msum_array <- array(yday(date) + round((mlen_array/2) - 1), dim=int_dims)
   } else {
     mlen <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
     msum <- cumsum(mlen) - mlen + 15
-    cyc <- array(c(1:12), dim=int_dims)[,1,1]
+    cyc <- array(c(1:12), dim=int_dims[1])
     mlen_array <- array(mlen, dim=int_dims)
     msum_array <- array(msum, dim=int_dims)
     warn$push('Assuming the data are monthly time series starting in January, all regular (non-leap) years.')
@@ -176,7 +168,7 @@ thornthwaite <- function(Tave, lat, na.rm=FALSE, verbose=TRUE) {
   if (out_type=='tsmatrix') {
     PET <- matrix(PET, nrow=n_times)
     PET <- ts(PET, frequency=ts_freq, start=ts_start)
-    colnames(PET) <- rep('PET_tho', ncol(ET0))
+    colnames(PET) <- rep('PET_tho', ncol(PET))
   } else if (out_type=='tsvector') {
     PET <- as.vector(PET)
     PET <- ts(PET, frequency=ts_freq, start=ts_start)
diff --git a/README.md b/README.md
index cffa523..6883c70 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # SPEI
 
-An R package incorporating a set of functions for computing potential or reference evapotranspiration and several widely used drought indices, currently including the **Standardized Precipitation-Evapotranspiration Index (SPEI)** and the **Standardized Precipitation Index (SPE)**.
+An R package incorporating a set of functions for computing potential or reference evapotranspiration and several widely used drought indices, currently including the **Standardized Precipitation-Evapotranspiration Index (SPEI)** and the **Standardized Precipitation Index (SPI)**.
 
 The package is centered on the SPEI. For more information on this drought index, please visit 
 [spei.csic.es](http://spei.csic.es/).
@@ -46,17 +46,17 @@ You can cite this references if you use the SPEI library on your work:
 
 Other (possibly useful) references:
 
-* Vicente-Serrano, S.M., Beguería, S., López-Moreno, J.I., Angulo, M., El Kenawy, A. 2010. A new global 0.5° gridded dataset (1901-2006) of a multiscalar drought index: comparison with current drought index datasets based on the Palmer Drought Severity Index. *Journal of Hydrometeorology* **11**: 1033--1043.
+* Vicente-Serrano, S.M., Beguería, S., López-Moreno, J.I., Angulo, M., El Kenawy, A. 2010. A new global 0.5° gridded dataset (1901-2006) of a multi-scalar drought index: comparison with current drought index datasets based on the Palmer Drought Severity Index. *Journal of Hydrometeorology* **11**: 1033--1043.
 
 * Beguería, S., Vicente-Serrano, S.M. y Angulo, M. 2010. A multi-scalar global drought data set: the SPEIbase: A new gridded product for the analysis of drought variability and impacts. *Bulletin of the American Meteorological Society* **91**: 1351--1354.
 
 * Vicente-Serrano, S.M., Beguería, S. and Juan I. López-Moreno. 2011. Comment on “Characteristics and trends in various forms of the Palmer Drought Severity Index (PDSI) during 1900-2008” by A. Dai. *Journal of Geophysical Research-Atmosphere* **116**: D19112, doi:10.1029/2011JD016410.
 
-* Vicente-Serrano, S.M., Santiago Beguería, Jorge Lorenzo-Lacruz, Jesús Julio Camarero, Juan I. López-Moreno, Cesar Azorin-Molina, Jesús Revuelto, Enrique Morán-Tejeda and Arturo Sánchez-Lorenzo. 2012. Performance of drought indices for ecological, agricultural and hydrological applications. *Earth Interactions* **16**: 1--27.
+* Vicente-Serrano, S.M., Santiago Beguería, Jorge Lorenzo-Lacruz, Jesús Julio Camarero, Juan I. López-Moreno, Cesar Azorín-Molina, Jesús Revuelto, Enrique Morán-Tejeda and Arturo Sánchez-Lorenzo. 2012. Performance of drought indices for ecological, agricultural and hydrological applications. *Earth Interactions* **16**: 1--27.
 
 * Vicente-Serrano, S.M., Célia Gouveia, Jesús Julio Camarero, Santiago Beguería, Ricardo Trigo, Juan I. López-Moreno, César Azorín-Molina, Edmond Pasho, Jorge Lorenzo-Lacruz, Jesús Revuelto, Enrique Morán-Tejeda and Arturo Sanchez-Lorenzo. 2012. The response of vegetation to drought time-scales across global land biomes. *Proceedings of the National Academy of Sciences of the United States of America*, doi: 10.1073/pnas.1207068110.
 
-* Vicente-Serrano, S.M., Gerard Van der Schrier, Santiago Beguería, Cesar Azorin-Molina, Juan-I. Lopez-Moreno. 2015. Contribution of precipitation and reference evapotranspiration to drought indices under different climates. Journal of Hydrology 426: 42--54.
+* Vicente-Serrano, S.M., Gerard Van der Schrier, Santiago Beguería, Cesar Azorín-Molina, Juan-I. Lopez-Moreno. 2015. Contribution of precipitation and reference evapotranspiration to drought indices under different climates. Journal of Hydrology 426: 42--54.
 
 * Vicente-Serrano, S.M., Beguería, S. 2016. Comment on "Candidate Distributions for Climatological Drought Indices (SPI and SPEI)" by James H. Stagge et al. *International Journal of Climatology* **36**: 2120--213.
 
@@ -64,7 +64,7 @@ Other (possibly useful) references:
 
 ## Version history
 
-### Version 1.8, May 2022 (current version on github).
+### Version 1.8.0, November 2022 (current version on GitHub, submitted to CRAN).
 
 1.  Solving several minor bugs in `<thornthwaite>`, `<hargreaves>`, and `<penman>` functions (output difference is lower than 0.1% with respect to version 1.7).
 2.  Solving a bug in `<spei>` that resulted in bad cumulative data when using a non-rectangular kernel, resulting in incorrect SPEI values.
@@ -75,9 +75,9 @@ Other (possibly useful) references:
 7.  Implementation of an option to include CO2 concentration data in function `<penman>`.
 8.  Implementation of a new option for when no wind data are available in function `<penman>`.
 9.  Functions `<spei>` and `<spi>` now admit time series of any frequency, and not only monthly (frequency 12) data.
-10. Fucntion `<plot.spei>` completely rewritten based on `ggplot2`, solving some bugs and enabling more flexibility.
+10. Function `<plot.spei>` completely rewritten based on `ggplot2`, solving some bugs and enabling more flexibility.
 
-### Version 1.7.2, January 2018 (latest stable version on github).
+### Version 1.7.2, January 2018 (only on GitHub).
 
 1. Several code optimizations and improvements (by github user @doug-friedman).
 2. Added formal unit testing with `testthat` (@doug-friedman).
@@ -110,7 +110,7 @@ Other (possibly useful) references:
 ### Version 1.4, May 2013.
 
 1. Minor fixes to functions `<penman>` and `<pwm>`.
-2. Documentation of the penman function defined by mistake ed as the saturation vapour pressure, while it should read 'actual vapour pressure'.
+2. Documentation of the penman function defined by mistake ed as the saturation vapor pressure, while it should read 'actual vapor pressure'.
 3. Function zzz.R added to display basic information about the SPEI package at startup.
 4. Function `<SPEINews>` added to display the NEWS file.
 
diff --git a/cran-comments.md b/cran-comments.md
new file mode 100644
index 0000000..858617d
--- /dev/null
+++ b/cran-comments.md
@@ -0,0 +1,5 @@
+## R CMD check results
+
+0 errors | 0 warnings | 1 note
+
+* This is a new release.
diff --git a/man/Datasets.Rd b/man/Datasets.Rd
index dc66c1f..6fe5eda 100644
--- a/man/Datasets.Rd
+++ b/man/Datasets.Rd
@@ -1,6 +1,5 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/data.R
-\docType{data}
 \name{Datasets}
 \alias{Datasets}
 \alias{wichita}
@@ -21,13 +20,11 @@
   \item{ACSH}{ monthly mean sun hours, in h.}
   \item{ACSH}{ monthly mean cloud cover, in \%.}
 }
-
 \code{balance} dataset: a data frame with monthly climatic water balance
  (precipitation minus potential evapotranspiration) at Indore (India),
  Kimberley (South Africa), Albuquerque (US), Valencia (Spain), Wien (Austria),
- Abashiri (Japan), Tampa (US), Sao Paulo (Brasil), Lahore (India), Punta
+ Abashiri (Japan), Tampa (US), Sao Paulo (Brazil), Lahore (India), Punta
  Arenas (Chile) and Helsinki (Finland), in mm.
-
 \code{cabinda} dataset: a data frame with one year of monthly climatic data
 at Cabinda (Angola, -5.33S 12.11E 20 m), with:
  
@@ -38,7 +35,7 @@ at Cabinda (Angola, -5.33S 12.11E 20 m), with:
   \item{RH}{ monthly mean relative humidity, in \%.}
   \item{U2}{ monthly mean wind speed, in km h-1}
   \item{tsun}{ monthly mean sunshine hours, in h.}
-  \item{Rs}{ monthly mean dialy incoming solar radiation, MJ m-2 d-1.}
+  \item{Rs}{ monthly mean diaily incoming solar radiation, MJ m-2 d-1.}
   \item{ET0}{ monthly ET0 from the original publication, in mm.}
 }
 
@@ -48,12 +45,6 @@ points from CRU TS 4.05 data set. The array has dimensions [time=1440,
 longitude=2, latitude=3], with time starting in January 1900. Longitudes
 are (0.25, 0.75), and latitudes (42.25, 42.75, 43.25), corresponding to the
 Central Pyrenees between Spain and France.
-
-An object of class \code{data.frame} with 12 rows and 8 columns.
-
-An object of class \code{data.frame} with 1296 rows and 11 columns.
-
-An object of class \code{array} of dimension 1440 x 2 x 3.
 }
 \source{
 The \code{wichita} data were obtained from the Global Historical Climatology Network 
@@ -63,17 +54,6 @@ Data for the \code{balance} dataset were taken from Allen et al. (1998), page 69
 The \code{cruts4} data were obtained from the CRU (Climatic Research Unit, University
 of East Anglia \url{https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.05/}) TS V4.05 data set.
 }
-\usage{
-data(wichita)
-data(balance)
-data(cabinda)
-
-cabinda
-
-balance
-
-cruts4
-}
 \description{
 Data used in the examples of the SPEI package: 
 \code{wichita} dataset: monthly climate in Wichita (Kansas, lat=37.6475, 
@@ -85,21 +65,9 @@ since January 1900;
 Allen et al. (1998);
 \code{cruts4}: 120 years of monthly climatic water balance (precipitation 
 minus reference evapotranspiration) data at six grid points from CRU TS 4.05.
-
-See wichita
-
-See wichita
-
-See wichita
 }
 \details{
 See description.
-
-See wichita
-
-See wichita
-
-See wichita
 }
 \examples{
 data(wichita)
@@ -107,8 +75,8 @@ names(wichita)
 summary(wichita)
 data(balance)
 summary(balance)
-
-
+data(cruts4)
+summary(cruts4)
 
 }
 \references{
diff --git a/man/Drought-indices.Rd b/man/Drought-indices.Rd
index b2860f6..85a5a88 100755
--- a/man/Drought-indices.Rd
+++ b/man/Drought-indices.Rd
@@ -3,14 +3,38 @@
 \name{Drought-indices}
 \alias{Drought-indices}
 \alias{spei}
+\alias{spei,}
 \alias{spi}
 \title{Calculation of the Standardized Precipitation-Evapotranspiration 
 Index (SPEI) and the Standardized Precipitation Index (SPI).}
 \usage{
-spei(data, scale, kernel = list(type = 'rectangular', shift = 0),
-distribution = 'log-Logistic', fit = 'ub-pwm', na.rm = FALSE,
-ref.start=NULL, ref.end=NULL, keep.x=FALSE, params=NULL, 
-verbose=TRUE, ...)
+spei(
+ data,
+ scale,
+ kernel = list(type = 'rectangular', shift = 0),
+ distribution = 'log-Logistic',
+ fit = 'ub-pwm',
+ na.rm = FALSE,
+ ref.start=NULL,
+ ref.end=NULL,
+ keep.x=FALSE,
+ params=NULL, 
+ verbose=TRUE,
+ ...)
+
+spi(
+ data,
+ scale,
+ kernel = list(type = 'rectangular', shift = 0),
+ distribution = 'Gamma',
+ fit = 'ub-pwm',
+ na.rm = FALSE,
+ ref.start=NULL,
+ ref.end=NULL,
+ keep.x=FALSE,
+ params=NULL, 
+ verbose=TRUE,
+ ...)
 
 spi(
   data,
@@ -21,8 +45,9 @@ spi(
   na.rm = FALSE,
   ref.start = NULL,
   ref.end = NULL,
-  x = FALSE,
+  keep.x = FALSE,
   params = NULL,
+  verbose = TRUE,
   ...
 )
 }
@@ -150,7 +175,7 @@ These parameters are then passed to the function \code{\link{kern}}.
 \section{Probability distributions}{
 
 Following the original definitions \code{spei} uses a log-Logistic distribution 
-by default, and \code{spi} uses a Gamma distribution. This behaviour can be modified, 
+by default, and \code{spi} uses a Gamma distribution. This behavior can be modified, 
 however, through parameter \code{distribution}.
 }
 
@@ -178,7 +203,7 @@ this option in order to know what distribution function should be used.
 
 \section{Reference period}{
 
-The default behaviour of the functions is using all the values provided in \code{data} 
+The default behavior of the functions is using all the values provided in \code{data} 
 for parameter fitting. However, this can be modified with help of parameters \code{ref.start} 
 and \code{ref.end}. These parameters allow defining a subset of values that will be used 
 for parameter fitting, i.e. a reference period. The functions, however, will compute the 
@@ -282,6 +307,7 @@ dim(spei_12$fitted)
 # Modding the plot
 # Since plot.spei() returns a ggplot object, it is possible to add or tweak
 # parts of the plot.
+require(ggplot2)
 plot(spei(wichita[,'BAL'], 12)) +
  ggtitle('SPEI1 at Wichita') +
  scale_fill_manual(values=c('blue','red')) +  # classic SPEI look
diff --git a/man/Generic-methods-for-spei-objects.Rd b/man/Generic-methods-for-spei-objects.Rd
index ceb1b33..e372d71 100755
--- a/man/Generic-methods-for-spei-objects.Rd
+++ b/man/Generic-methods-for-spei-objects.Rd
@@ -3,20 +3,17 @@
 \name{Generic-methods-for-spei-objects}
 \alias{Generic-methods-for-spei-objects}
 \alias{print.spei}
-\alias{print.spi}
 \alias{plot.spei}
-\alias{plot.spi}
 \alias{summary.spei}
-\alias{summary.spi}
 \title{Generic methods for \code{spei} objects.}
 \usage{
 \method{print}{spei}(x, ...)
 \method{summary}{spei}(object, ...)
-\method{plot}{spei}(x, ttext, ...)
+\method{plot}{spei}(x, ...)
 
 \method{summary}{spei}(object, ...)
 
-\method{plot}{spei}(x)
+\method{plot}{spei}(x, ...)
 }
 \arguments{
 \item{x}{an object of class \code{spei}.}
@@ -24,8 +21,6 @@
 \item{...}{additional parameters, not used at present.}
 
 \item{object}{an object of class \code{spei}.}
-
-\item{ttext}{text to use as part of the plot title}
 }
 \description{
 Generic methods for extracting information and plotting \code{spei} objects.
@@ -48,7 +43,7 @@ See print.spei
 See print.spei
 }
 \examples{
-See examples of use in the help page of the \code{spei()} function.
+# See examples of use in the help page of the spei() function.
 
 }
 \references{
diff --git a/man/Kernel-functions.Rd b/man/Kernel-functions.Rd
index 37eba4a..de76b0d 100755
--- a/man/Kernel-functions.Rd
+++ b/man/Kernel-functions.Rd
@@ -18,7 +18,7 @@ kern.plot(scale = 12, shift = 0)
 \item{shift}{numeric, shifting of the kernel peak.}
 }
 \value{
-A vector of lenght equal to \code{scale} with weights used for computing the drought index.
+A vector of length equal to \code{scale} with weights used for computing the drought index.
 }
 \description{
 Function \code{kern} is used internally by \code{\link{spei}} and
@@ -29,7 +29,7 @@ See kern
 \details{
 Drought indices, such as the SPEI or the SPI, are usually 
 computed at different time scales to adapt to the different response 
-times of systems affected by drought. This is acomplished by applying 
+times of systems affected by drought. This is accomplished by applying 
 a kernel function to the data prior to computation of the SPEI. 
 Application of a kernel has the effect of smoothing the temporal 
 variability of the resulting SPEI, allowing for the major patterns 
diff --git a/man/Potential-evapotranspiration.Rd b/man/Potential-evapotranspiration.Rd
index bd95887..f0be80a 100755
--- a/man/Potential-evapotranspiration.Rd
+++ b/man/Potential-evapotranspiration.Rd
@@ -9,12 +9,12 @@
 \usage{
 thornthwaite(Tave, lat, na.rm = FALSE, verbose=TRUE)
 
-hargreaves(Tmin, Tmax, Ra = NA, lat = NA, Pre = NA, na.rm = FALSE, verbose=TRUE)
+hargreaves(Tmin, Tmax, Ra = NULL, lat = NULL, Pre = NULL, na.rm = FALSE, verbose=TRUE)
 
-penman(Tmin, Tmax, U2, Ra = NA, lat = NA, Rs = NA, tsun = NA,
-       CC = NA, ed = NA, Tdew = NA, RH = NA, P = NA, P0 = NA,
-              z = NA, crop='short', na.rm = FALSE, method='ICID',
-              verbose=TRUE)
+penman(Tmin, Tmax, U2, Ra = NULL, lat = NULL, Rs = NULL, tsun = NULL,
+       CC = NULL, ed = NULL, Tdew = NULL, RH = NULL, P = NULL, P0 = NULL,
+       CO2 = NULL, z = NULL, crop='short', na.rm = FALSE, method='ICID',
+       verbose=TRUE)
 
 penman(
   Tmin,
@@ -57,13 +57,13 @@ thornthwaite(Tave, lat, na.rm = FALSE, verbose = TRUE)
 
 \item{U2}{a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily wind speeds at 2 m height, m s-1.}
 
-\item{Rs}{optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean dialy incoming solar radiation, MJ m-2 d-1.}
+\item{Rs}{optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily incoming solar radiation, MJ m-2 d-1.}
 
 \item{tsun}{optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily bright sunshine hours, h.}
 
 \item{CC}{optional, numeric a vector, matrix or time series of monthly mean cloud cover, \%.}
 
-\item{ed}{optional, numeric a vector, matrix or time series of monthly mean actual vapour pressure at 2 m height, kPa.}
+\item{ed}{optional, numeric a vector, matrix or time series of monthly mean actual vapor pressure at 2 m height, kPa.}
 
 \item{Tdew}{optional, a numeric vector, tsvector, matrix, tsmatrix, or 3-d array of monthly mean daily dewpoint temperature (used for estimating ed), ºC.}
 
@@ -263,7 +263,7 @@ plot(ts(cbind(pen_300, pen_co2), fr=12))
 }
 \references{
 Thornthwaite, C. W., 1948. An approach toward a rational classification of climate. 
-\emph{Geographical Review} \bold{38}: 55–94. doi:10.2307/2107309.
+\emph{Geographical Review} \bold{38}: 55–94. DOI:10.2307/2107309.
 
 Hargreaves G.H., 1994. Defining and using reference evapotranspiration. 
 \emph{Journal of Irrigation and Drainage Engineering} \bold{120}: 1132–1139.
@@ -285,7 +285,7 @@ Rep. Task Com. on Standardized Reference Evapotranspiration July 9, 2002, EWRI
 Reston, VA, 57 pp.
 
 Yang, Y., Roderick, M.L., Zhang, S. McVicar, T., Donohue, R.J., 2019. Hydrologic implications of vegetation
-response to elevated CO2 in climate projections. \emph{Nature Clim Change} \bold{9}: 44–48.
+response to elevated CO2 in climate projections. \emph{Nature Climate Change} \bold{9}: 44–48.
 }
 \author{
 Santiago Beguería
diff --git a/tests/testthat/data/penman_out_tsun.rds b/tests/testthat/data/penman_out_tsun.rds
index 3120793..dca2830 100644
Binary files a/tests/testthat/data/penman_out_tsun.rds and b/tests/testthat/data/penman_out_tsun.rds differ
diff --git a/tests/testthat/data/penman_out_tsun_ts.rds b/tests/testthat/data/penman_out_tsun_ts.rds
index dae193b..e946171 100644
Binary files a/tests/testthat/data/penman_out_tsun_ts.rds and b/tests/testthat/data/penman_out_tsun_ts.rds differ
diff --git a/tests/testthat/data/spi_12_array.rds b/tests/testthat/data/spi_12_array.rds
new file mode 100644
index 0000000..e9d84c0
Binary files /dev/null and b/tests/testthat/data/spi_12_array.rds differ
diff --git a/tests/testthat/data/spi_12_matrix.rds b/tests/testthat/data/spi_12_matrix.rds
new file mode 100644
index 0000000..7bab2d0
Binary files /dev/null and b/tests/testthat/data/spi_12_matrix.rds differ
diff --git a/tests/testthat/data/spi_12_tsmatrix.rds b/tests/testthat/data/spi_12_tsmatrix.rds
new file mode 100644
index 0000000..5ff2b93
Binary files /dev/null and b/tests/testthat/data/spi_12_tsmatrix.rds differ
diff --git a/tests/testthat/data/spi_12_tsvector.rds b/tests/testthat/data/spi_12_tsvector.rds
new file mode 100644
index 0000000..6b11879
Binary files /dev/null and b/tests/testthat/data/spi_12_tsvector.rds differ
diff --git a/tests/testthat/data/spi_12_vector.rds b/tests/testthat/data/spi_12_vector.rds
new file mode 100644
index 0000000..6d24dea
Binary files /dev/null and b/tests/testthat/data/spi_12_vector.rds differ
diff --git a/tests/testthat/data/spi_1_array.rds b/tests/testthat/data/spi_1_array.rds
new file mode 100644
index 0000000..28902e3
Binary files /dev/null and b/tests/testthat/data/spi_1_array.rds differ
diff --git a/tests/testthat/data/spi_1_matrix.rds b/tests/testthat/data/spi_1_matrix.rds
new file mode 100644
index 0000000..b68981c
Binary files /dev/null and b/tests/testthat/data/spi_1_matrix.rds differ
diff --git a/tests/testthat/data/spi_1_tsmatrix.rds b/tests/testthat/data/spi_1_tsmatrix.rds
new file mode 100644
index 0000000..0499e2a
Binary files /dev/null and b/tests/testthat/data/spi_1_tsmatrix.rds differ
diff --git a/tests/testthat/data/spi_1_tsvector.rds b/tests/testthat/data/spi_1_tsvector.rds
new file mode 100644
index 0000000..31840d8
Binary files /dev/null and b/tests/testthat/data/spi_1_tsvector.rds differ
diff --git a/tests/testthat/data/spi_1_vector.rds b/tests/testthat/data/spi_1_vector.rds
new file mode 100644
index 0000000..07b8db9
Binary files /dev/null and b/tests/testthat/data/spi_1_vector.rds differ
diff --git a/tests/testthat/data/spi_gamma.rds b/tests/testthat/data/spi_gamma.rds
new file mode 100644
index 0000000..07b8db9
Binary files /dev/null and b/tests/testthat/data/spi_gamma.rds differ
diff --git a/tests/testthat/data/spi_pe3.rds b/tests/testthat/data/spi_pe3.rds
new file mode 100644
index 0000000..fefd8c8
Binary files /dev/null and b/tests/testthat/data/spi_pe3.rds differ
diff --git a/tests/testthat/data/spi_refend.rds b/tests/testthat/data/spi_refend.rds
new file mode 100644
index 0000000..4b82df8
Binary files /dev/null and b/tests/testthat/data/spi_refend.rds differ
diff --git a/tests/testthat/data/spi_refst.rds b/tests/testthat/data/spi_refst.rds
new file mode 100644
index 0000000..afd53ae
Binary files /dev/null and b/tests/testthat/data/spi_refst.rds differ
diff --git a/tests/testthat/data/spi_refst_refend.rds b/tests/testthat/data/spi_refst_refend.rds
new file mode 100644
index 0000000..8818cd9
Binary files /dev/null and b/tests/testthat/data/spi_refst_refend.rds differ
diff --git a/tests/testthat/test_hargreaves.R b/tests/testthat/test_hargreaves.R
index 8d7b077..ff3a25e 100755
--- a/tests/testthat/test_hargreaves.R
+++ b/tests/testthat/test_hargreaves.R
@@ -65,7 +65,7 @@ test_that('input more than three dimensions', {
 
 test_that('Tmin and Tmax different size', {
   expect_error(hargreaves(c(1,TMIN), TMAX, lat=37.6475),
-               '`Tmin` and `Tmax`cannot have different lengths.')
+               '`Tmin` and `Tmax` cannot have different lengths.')
 })
 
 test_that('Ra incorrect size', {
@@ -107,7 +107,7 @@ test_that('example with lat', {
 #out <- hargreaves(TMIN, TMAX, Ra=TMAX)
 #saveRDS(out, file='./tests/testthat/data/hargreaves_out_ra.rds')
 test_that('example with Ra', {
-  harModOut = readRDS('data/hargreaves_out_ra.rds')
+  harModOut <- readRDS('data/hargreaves_out_ra.rds')
   expect_equal(harModOut, hargreaves(TMIN, TMAX, Ra=TMAX)
   )
 })
@@ -115,7 +115,7 @@ test_that('example with Ra', {
 #out <- hargreaves(TMIN, TMAX, lat=37.6475, Pre=TMAX)
 #saveRDS(out, file='./tests/testthat/data/hargreaves_out_pre.rds')
 test_that('example with Pre', {
-  harModOut = readRDS('data/hargreaves_out_pre.rds')
+  harModOut <- readRDS('data/hargreaves_out_pre.rds')
   expect_equal(harModOut, hargreaves(TMIN, TMAX, lat=37.6475, Pre=TMAX)
   )
 })
diff --git a/tests/testthat/test_kern.R b/tests/testthat/test_kern.R
index c57d0e9..e7a23b0 100755
--- a/tests/testthat/test_kern.R
+++ b/tests/testthat/test_kern.R
@@ -2,13 +2,13 @@ context("kern")
 
 
 test_that("example w/ defaults", {
-  k12Out = readRDS("data/k12_Out.rds")
+  k12Out <- readRDS("data/k12_Out.rds")
   expect_equal(k12Out, kern(12))
 })
 
 
 test_that("example w/ defaults", {
-  k12GaussOut = readRDS("data/k12_gauss_Out.rds")
+  k12GaussOut <- readRDS("data/k12_gauss_Out.rds")
   expect_equal(k12GaussOut, kern(12, 'gaussian'))
 })
 
diff --git a/tests/testthat/test_penman.R b/tests/testthat/test_penman.R
index 915f4fc..a90c2cf 100755
--- a/tests/testthat/test_penman.R
+++ b/tests/testthat/test_penman.R
@@ -32,9 +32,9 @@ test_that('No Rs nor CC provided', {
                'One of `Rs`, the pair `tsun` and `lat`, or `CC` must be provided.')
 })
 
-test_that('No Rs nor CC provided', {
+test_that('No Tmin provided', {
   expect_error(penman(Tmax=TMAX, U2=AWND, tsun=TSUN, lat=37.6475, z=402.6, na.rm=TRUE),
-               'One of `ed`, `Tdew`, `RH` or `Tmin` must be provided.')
+               'argument \"Tmin\" is missing, with no default')
 })
 
 test_that('No P nor z provided', {
@@ -59,12 +59,12 @@ test_that('NAs in Tmin', {
 
 test_that('NAs in Tmax', {
   expect_error(penman(c(1,TMIN), c(NA,TMAX), c(1,TMIN), tsun=c(1,TMIN), lat=37.6475, z=402.6, na.rm=FALSE),
-               'Data must not contain NA values if argument `na.rm` is set to FALSE.')
+               '`Tmin`, `Tmax` and `U2` must not contain NA values if argument `na.rm` is set to FALSE.')
 })
 
 test_that('NAs in U2', {
   expect_error(penman(c(1,TMIN), c(1,TMAX), c(NA,TMIN), tsun=c(1,TMIN), lat=37.6475, z=402.6, na.rm=FALSE),
-               'Data must not contain NA values if argument `na.rm` is set to FALSE.')
+               '`Tmin`, `Tmax` and `U2` must not contain NA values if argument `na.rm` is set to FALSE.')
 })
 
 test_that('NAs in Tsun', {
@@ -89,7 +89,7 @@ tmin <- array(TMIN, dim=c(nt, 2, 2, 2))
 # tsun <- array(TSUN, dim=c(nt, 2, 2, 2))
 test_that('input more than three dimensions', {
   expect_error(penman(tmin, TMAX, AWND, tsun=TSUN, lat=37.6475, z=402.6, na.rm=TRUE),
-               'Input data can not have more than three dimensions.')
+               'Input data can not have more than 3 dimensions.')
 })
 
 test_that('Tmax incorrect length', {
diff --git a/tests/testthat/test_spei.R b/tests/testthat/test_spei.R
index 4dc95af..2c94a79 100755
--- a/tests/testthat/test_spei.R
+++ b/tests/testthat/test_spei.R
@@ -10,9 +10,9 @@ data(balance)
 data(cruts4)
 
 x_vec <- as.numeric(wichita$BAL)
-x_tsvec <- ts(bal_vec, c(1980,1), fr=12)
+x_tsvec <- ts(x_vec, c(1980,1), fr=12)
 x_mat <- as.matrix(balance)
-x_tsmat <- ts(bal_mat, c(1980,1), fr=12)
+x_tsmat <- ts(x_mat, c(1980,1), fr=12)
 x_array <- cruts4
 
 # Test input data checks and error messages
@@ -91,7 +91,7 @@ test_that('Incorrect data type: passing a data.frame', {
 #saveRDS(out, file='./tests/testthat/data/spei_1_vector.rds')
 test_that('example with vector data, scale 1', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_vector.rds'),
+    readRDS('data/spei_1_vector.rds'),
     spei(x_vec, 1)$fitted
   )
 })
@@ -100,7 +100,7 @@ test_that('example with vector data, scale 1', {
 #saveRDS(out, file='./tests/testthat/data/spei_12_vector.rds')
 test_that('example with vector data, scale 12', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_12_vector.rds'),
+    readRDS('data/spei_12_vector.rds'),
     spei(x_vec, 12)$fitted
   )
 })
@@ -109,7 +109,7 @@ test_that('example with vector data, scale 12', {
 #saveRDS(out, file='./tests/testthat/data/spei_1_tsvector.rds')
 test_that('example with tsvector data, scale 1', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_tsvector.rds'),
+    readRDS('data/spei_1_tsvector.rds'),
     spei(x_tsvec, 1)$fitted
   )
 })
@@ -118,7 +118,7 @@ test_that('example with tsvector data, scale 1', {
 #saveRDS(out, file='./tests/testthat/data/spei_12_tsvector.rds')
 test_that('example with tsvector data, scale 12', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_12_tsvector.rds'),
+    readRDS('data/spei_12_tsvector.rds'),
     spei(x_tsvec, 12)$fitted
   )
 })
@@ -127,7 +127,7 @@ test_that('example with tsvector data, scale 12', {
 #saveRDS(out, file='./tests/testthat/data/spei_1_matrix.rds')
 test_that('example with matrix data, scale 1', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_matrix.rds'),
+    readRDS('data/spei_1_matrix.rds'),
     spei(x_mat, 1)$fitted
   )
 })
@@ -136,7 +136,7 @@ test_that('example with matrix data, scale 1', {
 #saveRDS(out, file='./tests/testthat/data/spei_12_matrix.rds')
 test_that('example with matrix data, scale 12', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_12_matrix.rds'),
+    readRDS('data/spei_12_matrix.rds'),
     spei(x_mat, 12)$fitted
   )
 })
@@ -145,7 +145,7 @@ test_that('example with matrix data, scale 12', {
 #saveRDS(out, file='./tests/testthat/data/spei_1_tsmatrix.rds')
 test_that('example with tsmatrix data, scale 1', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_tsmatrix.rds'),
+    readRDS('data/spei_1_tsmatrix.rds'),
     spei(x_tsmat, 1)$fitted
   )
 })
@@ -154,7 +154,7 @@ test_that('example with tsmatrix data, scale 1', {
 #saveRDS(out, file='./tests/testthat/data/spei_12_tsmatrix.rds')
 test_that('example with tsmatrix data, scale 12', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_12_tsmatrix.rds'),
+    readRDS('data/spei_12_tsmatrix.rds'),
     spei(x_tsmat, 12)$fitted
   )
 })
@@ -163,7 +163,7 @@ test_that('example with tsmatrix data, scale 12', {
 #saveRDS(out, file='./tests/testthat/data/spei_1_array.rds')
 test_that('example with 3-d array data, scale 1', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_array.rds'),
+    readRDS('data/spei_1_array.rds'),
     spei(x_array, 1)$fitted
   )
 })
@@ -172,7 +172,7 @@ test_that('example with 3-d array data, scale 1', {
 #saveRDS(out, file='./tests/testthat/data/spei_12_array.rds')
 test_that('example with 3-d array data, scale 12', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_12_array.rds'),
+    readRDS('data/spei_12_array.rds'),
     spei(x_array, 12)$fitted
   )
 })
@@ -184,7 +184,7 @@ re <- c(25,12)
 #saveRDS(out, file='./tests/testthat/data/spei_refst_refend.rds')
 test_that('example with both ref.start and ref.end', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_refst_refend.rds'),
+    readRDS('data/spei_refst_refend.rds'),
     spei(x_vec, 12, ref.start=rs, ref.end=re)$fitted
   )
 })
@@ -193,7 +193,7 @@ test_that('example with both ref.start and ref.end', {
 #saveRDS(out, file='./tests/testthat/data/spei_refst.rds')
 test_that('example with only ref.start', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_refst.rds'),
+    readRDS('data/spei_refst.rds'),
     spei(x_vec, 12, ref.start=rs)$fitted
   )
 })
@@ -202,7 +202,7 @@ test_that('example with only ref.start', {
 #saveRDS(out, file='./tests/testthat/data/spei_refend.rds')
 test_that('example with only ref.end', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_refend.rds'),
+    readRDS('data/spei_refend.rds'),
     spei(x_vec, 12, ref.end=rs)$fitted
   )
 })
@@ -211,7 +211,7 @@ test_that('example with only ref.end', {
 #saveRDS(out, file='./tests/testthat/data/spei_gamma.rds')
 test_that('example with Gamma distribution', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_gamma.rds'),
+    readRDS('data/spei_gamma.rds'),
     spei(x_vec, 1, distribution='Gamma')$fitted
   )
 })
@@ -220,7 +220,7 @@ test_that('example with Gamma distribution', {
 #saveRDS(out, file='./tests/testthat/data/spei_pe3.rds')
 test_that('example with PearsonIII distribution', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_pe3.rds'),
+    readRDS('data/spei_pe3.rds'),
     spei(x_vec, 1, distribution='PearsonIII')$fitted # , tol = 1e-7 
   )
 })
@@ -228,7 +228,7 @@ test_that('example with PearsonIII distribution', {
 pars <- spei(x_vec, 1)$coefficients
 test_that('example with user-provided parameters', {
   expect_equal(
-    readRDS('./tests/testthat/data/spei_1_vector.rds'),
+    readRDS('data/spei_1_vector.rds'),
     spei(x_vec, 1, params=pars)$fitted
   )
 })
diff --git a/tests/testthat/test_spi.R b/tests/testthat/test_spi.R
index 67f5e40..231dd42 100755
--- a/tests/testthat/test_spi.R
+++ b/tests/testthat/test_spi.R
@@ -1,30 +1,213 @@
-context("spi")
+context('spi')
 
 data(wichita)
-wichita$PET <- thornthwaite(wichita$TMED,37.6475)
+attach(wichita)
 
-# Changes in L-Moments calculation yield differences less than 1e-7
+x_vec <- as.numeric(wichita$PRCP)
+x_tsvec <- ts(x_vec, c(1980,1), fr=12)
+x_mat <- matrix(x_vec, nrow=length(x_vec), ncol=4)
+x_tsmat <- ts(x_mat, c(1980,1), fr=12)
+x_array <- array(x_vec, dim=c(length(x_vec), 2, 2))
 
-test_that("example w/ 1 month", {
-  spi1MoOut = readRDS("data/spi_1mo_Out.rds")
-  expect_equal(spi1MoOut, 
-               spi(wichita$PRCP,1),
-               tol = 1e-7
+# Test input data checks and error messages
+
+#spi(x_vec, 1)
+
+test_that('Non-numeric scale', {
+  expect_error(spi(x_vec, scale='1'),
+               'Argument `scale` must be numeric.')
+})
+
+test_that('Bad distribution', {
+  expect_error(spi(x_vec, 1, distribution='moments'),
+               'Argument `distribution` must be one of `log-Logistic`, `Gamma` or `PearsonIII`.')
+})
+
+test_that('Bad fit', {
+  expect_error(spi(x_vec, 1, fit='logistic'),
+               'Argument `fit` must be one of `ub-pwm`, `pp-pwm` or `max-lik`.')
+})
+
+test_that('Bad ref.start', {
+  expect_error(spi(x_vec, 1, ref.start=1),
+               'Argument `ref.start` must be a numeric vector of length two.')
+})
+
+test_that('Bad ref.start 2', {
+  expect_error(spi(x_vec, 1, ref.start=c('a','b')),
+               'Argument `ref.start` must be a numeric vector of length two.')
+})
+
+test_that('Bad ref.end', {
+  expect_error(spi(x_vec, 1, ref.end=1),
+               'Argument `ref.end` must be a numeric vector of length two.')
+})
+
+test_that('Bad ref.end 2', {
+  expect_error(spi(x_vec, 1, ref.end=c('a','b')),
+               'Argument `ref.end` must be a numeric vector of length two.')
+})
+
+test_that('There are NAs when na.rm=FALSE', {
+  expect_error(spi(c(NA,x_vec), 1),
+               '`data` must not contain NA values if argument `na.rm` is set to FALSE.')
+})
+
+test_that('Incorrect input dimensions', {
+  expect_error(spi(array(x_vec, c(382,2,2,2)), 1),
+               'Input data can not have more than three dimensions.')
+})
+
+test_that('Incorrect data type: passing a data.frame', {
+  expect_error(spi(as.data.frame(x_mat), 1),
+               'Bad data type: input must be a vector, tsvector, matrix, tsmatrix, or 3-d array.')
+})
+
+
+# Test function results
+
+#out <- spi(x_vec, 1)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_1_vector.rds')
+test_that('example with vector data, scale 1', {
+  expect_equal(
+    readRDS('data/spi_1_vector.rds'),
+    spi(x_vec, 1)$fitted
+  )
+})
+
+#out <- spi(x_vec, 12)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_12_vector.rds')
+test_that('example with vector data, scale 12', {
+  expect_equal(
+    readRDS('data/spi_12_vector.rds'),
+    spi(x_vec, 12)$fitted
+  )
+})
+
+#out <- spi(x_tsvec, 1)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_1_tsvector.rds')
+test_that('example with tsvector data, scale 1', {
+  expect_equal(
+    readRDS('data/spi_1_tsvector.rds'),
+    spi(x_tsvec, 1)$fitted
+  )
+})
+
+#out <- spi(x_tsvec, 12)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_12_tsvector.rds')
+test_that('example with tsvector data, scale 12', {
+  expect_equal(
+    readRDS('data/spi_12_tsvector.rds'),
+    spi(x_tsvec, 12)$fitted
+  )
+})
+
+#out <- spi(x_mat, 1)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_1_matrix.rds')
+test_that('example with matrix data, scale 1', {
+  expect_equal(
+    readRDS('data/spi_1_matrix.rds'),
+    spi(x_mat, 1)$fitted
+  )
+})
+
+#out <- spi(x_mat, 12)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_12_matrix.rds')
+test_that('example with matrix data, scale 12', {
+  expect_equal(
+    readRDS('data/spi_12_matrix.rds'),
+    spi(x_mat, 12)$fitted
+  )
+})
+
+#out <- spi(x_tsmat, 1)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_1_tsmatrix.rds')
+test_that('example with tsmatrix data, scale 1', {
+  expect_equal(
+    readRDS('data/spi_1_tsmatrix.rds'),
+    spi(x_tsmat, 1)$fitted
+  )
+})
+
+#out <- spi(x_tsmat, 12)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_12_tsmatrix.rds')
+test_that('example with tsmatrix data, scale 12', {
+  expect_equal(
+    readRDS('data/spi_12_tsmatrix.rds'),
+    spi(x_tsmat, 12)$fitted
+  )
+})
+
+#out <- spi(x_array, 1)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_1_array.rds')
+test_that('example with 3-d array data, scale 1', {
+  expect_equal(
+    readRDS('data/spi_1_array.rds'),
+    spi(x_array, 1)$fitted
+  )
+})
+
+#out <- spi(x_array, 12)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_12_array.rds')
+test_that('example with 3-d array data, scale 12', {
+  expect_equal(
+    readRDS('data/spi_12_array.rds'),
+    spi(x_array, 12)$fitted
+  )
+})
+
+rs <- c(10,1)
+re <- c(25,12)
+
+#out <- spi(x_vec, 12, ref.start=rs, ref.end=re)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_refst_refend.rds')
+test_that('example with both ref.start and ref.end', {
+  expect_equal(
+    readRDS('data/spi_refst_refend.rds'),
+    spi(x_vec, 12, ref.start=rs, ref.end=re)$fitted
+  )
+})
+
+#out <- spi(x_vec, 12, ref.start=rs)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_refst.rds')
+test_that('example with only ref.start', {
+  expect_equal(
+    readRDS('data/spi_refst.rds'),
+    spi(x_vec, 12, ref.start=rs)$fitted
+  )
+})
+
+#out <- spi(x_vec, 12, ref.end=rs)$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_refend.rds')
+test_that('example with only ref.end', {
+  expect_equal(
+    readRDS('data/spi_refend.rds'),
+    spi(x_vec, 12, ref.end=rs)$fitted
+  )
+})
+
+#out <- spi(x_vec, 1, distribution='Gamma')$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_gamma.rds')
+test_that('example with Gamma distribution', {
+  expect_equal(
+    readRDS('data/spi_gamma.rds'),
+    spi(x_vec, 1, distribution='Gamma')$fitted
   )
 })
 
-test_that("example w/ 12 months", {
-  spi12MoOut = readRDS("data/spi_12mo_Out.rds")
-  expect_equal(spi12MoOut, 
-               spi(wichita$PRCP,12),
-               tol = 1e-7
+#out <- spi(x_vec, 1, distribution='PearsonIII')$fitted
+#saveRDS(out, file='./tests/testthat/data/spi_pe3.rds')
+test_that('example with PearsonIII distribution', {
+  expect_equal(
+    readRDS('data/spi_pe3.rds'),
+    spi(x_vec, 1, distribution='PearsonIII')$fitted # , tol = 1e-7 
   )
 })
 
-test_that("example w/ 1 month (SPEI vs SPI)", {
-  spi1MoOut = readRDS("data/spi_1mo_Out.rds")
-  expect_equal(spi1MoOut[-1], 
-               spei(wichita$PRCP,1,distribution="Gamma")[-1],
-               tol = 1e-7
+pars <- spi(x_vec, 1)$coefficients
+test_that('example with user-provided parameters', {
+  expect_equal(
+    readRDS('data/spi_1_vector.rds'),
+    spi(x_vec, 1, params=pars)$fitted
   )
 })
diff --git a/tests/testthat/test_thornthwaite.R b/tests/testthat/test_thornthwaite.R
index a8cfd4c..1523479 100755
--- a/tests/testthat/test_thornthwaite.R
+++ b/tests/testthat/test_thornthwaite.R
@@ -37,13 +37,13 @@ test_that('lat has incorrect length', {
 #out <- thornthwaite(TMED, 37.6475)
 #saveRDS(out, file='./tests/testthat/data/thornthwaite_out.rds')
 test_that('example', {
-  thoOut = readRDS('data/thornthwaite_out.rds') #./tests/testthat/
+  thoOut <- readRDS('data/thornthwaite_out.rds') #./tests/testthat/
   expect_equal(thoOut, thornthwaite(TMED, 37.6475))
 })
 
 #out <- thornthwaite(ts(TMED, c(1980,1), fr=12), 37.6475)
 #saveRDS(out, file='./tests/testthat/data/thornthwaite_out_ts.rds')
 test_that('example', {
-  thoOut = readRDS('data/thornthwaite_out_ts.rds') #./tests/testthat/
+  thoOut <- readRDS('data/thornthwaite_out_ts.rds') #./tests/testthat/
   expect_equal(thoOut, thornthwaite(ts(TMED, c(1980,1), fr=12), 37.6475))
 })