diff --git a/DESCRIPTION b/DESCRIPTION
index 2f5c132250a2ae2d4498ef534675b84b1da24c27..76a298d58ef71deea74de7b70e0b6b606022d10f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -6,6 +6,7 @@ Authors@R: c(
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")),
     person("Carmen", "Alvarez-Castro", , "carmen.alvarez-castro@cmcc.it", role = "aut"),
+    person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "aut"),
     person("Jost", "von Hardenberg", , email = c("j.vonhardenberg@isac.cnr.it", "jost.hardenberg@polito.it"), role = "aut",  comment = c(ORCID = "0000-0002-5312-8070")),
     person("Llorenç", "LLedo", , "llledo@bsc.es", role = "aut"), 
     person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "aut"),
@@ -13,16 +14,17 @@ Authors@R: c(
     person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "aut"),
     person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"),
     person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"),
-    person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"),
     person("Filippo", "Cali Quaglia", , "filippo.caliquaglia@gmail.com", role = "ctb"),
     person("Chihchung", "Chou", , "chihchung.chou@bsc.es", role = "ctb"),
     person("Nicola", "Cortesi", , "nicola.cortesi@bsc.es", role = "ctb"),
     person("Susanna", "Corti", , "s.corti@isac.cnr.it", role = "ctb"),
     person("Paolo", "Davini", , "p.davini@isac.cnr.it", role = "ctb"),
+    person("Gildas", "Dayon", , "gildas.dayon@meteo.fr", role = "ctb"),
     person("Marta", "Dominguez", , "mdomingueza@aemet.es", role = "ctb"),
     person("Federico", "Fabiano", , "f.fabiano@isac.cnr.it", role = "ctb"),
     person("Ignazio", "Giuntoli", , "i.giuntoli@isac.cnr.it", role = "ctb"),
     person("Raul", "Marcos", , "raul.marcos@bsc.es", role = "ctb"),
+    person("Paola", "Marson", , "paola.marson@meteo.fr", role = "ctb"),
     person("Niti", "Mishra", , "niti.mishra@bsc.es", role = "ctb"),
     person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"),
     person("Francesc", "Roura-Adserias", , "francesc.roura@bsc.es", role = "ctb"),
@@ -43,6 +45,7 @@ Description: Exploits dynamical seasonal forecasts in order to provide
     Terzago et al. (2018) <doi:10.5194/nhess-18-2825-2018>.
     Torralba et al. (2017) <doi:10.1175/JAMC-D-16-0204.1>.
     D'Onofrio et al. (2014) <doi:10.1175/JHM-D-13-096.1>.
+    Verfaillie et al. (2017) <doi:10.5194/gmd-10-4257-2017>.
     Van Schaeybroeck et al. (2019) <doi:10.1016/B978-0-12-812372-0.00010-8>.
     Yiou et al. (2013) <doi:10.1007/s00382-012-1626-3>.
 Depends:
@@ -79,4 +82,4 @@ VignetteBuilder: knitr
 License: Apache License 2.0
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.1.1
+RoxygenNote: 7.0.2
diff --git a/R/CST_AdamontAnalog.R b/R/CST_AdamontAnalog.R
new file mode 100644
index 0000000000000000000000000000000000000000..104913acd04a8bab563a67999963dc3c7f7dd925
--- /dev/null
+++ b/R/CST_AdamontAnalog.R
@@ -0,0 +1,258 @@
+#'CST_AdamontAnalog finds analogous data in the reference dataset to experiment
+#'data based on weather types
+#'
+#'@description This function searches for analogs in a reference dataset for 
+#'experiment data, based on corresponding weather types. The experiment data is
+#'typically a hindcast, observations are typically provided by reanalysis data.
+#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+#'
+#'@param exp experiment data an object of class \code{s2dv_cube}, can be output 
+#'from quantile correction using CST_AdamontQQCorr.R
+#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data} 
+#'but lat/lon)
+#'@param obs reference data, also of class \code{s2dv_cube}. Note that lat/lon 
+#'dimensions need to be the same as \code{exp}
+#'@param wt_obs corresponding weather types (same dimensions as \code{obs$data}
+#'but lat/lon)
+#'@param nanalogs integer defining the number of analog values to return
+#'(default: 5)
+#'@param method a character string indicating the method used for analog
+#'definition
+#'       Coded are 'pattcorr': pattern correlation
+#'	           'rain1' (for precip patterns): rain occurrence consistency
+#'	           'rain01' (for precip patterns): rain occurrence/non 
+#'                        occurrence consistency
+#'@param thres real number indicating the threshold to define rain 
+#'occurrence/non occurrence in rain(0)1
+#'@param search_obsdims list of dimensions in \code{obs} along which analogs are
+#'searched for
+#'@param londim name of longitude dimension
+#'@param latdim name of latitude dimension
+#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs
+#'analog values for each value of \code{exp} input data
+#'@import multiApply
+#'@importFrom s2dverification Subset
+#'@examples
+#'\dontrun{
+#'wt_exp <- sample(1:3, 15*6*3, replace=T)
+#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+#'wt_obs <- sample(1:3, 6*3, replace=T)
+#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+# analog_vals <- CST_AdamontAnalog(exp=lonlat_data$exp, obs=lonlat_data$obs, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2)
+#'}
+CST_AdamontAnalog <- function(exp, obs, wt_exp, wt_obs, nanalogs, 
+			      method = 'pattcorr', thres = NULL, 
+			      search_obsdims = c('member', 'sdate', 'ftime'), 
+			      londim = 'lon', latdim = 'lat') {
+
+  dimnames <- names(dim(obs$data))
+  dimnamesexp <- names(dim(exp$data))
+  if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) {
+     stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ",
+	  "as output by CSTools::CST_Load.")
+  }
+  if (!(method %in% c('pattcorr','rain1','rain01'))) {
+     stop("Input parameter 'method' must be 'pattcorr', 'rain1', or 'rain01'")
+  }
+  if (is.null(nanalogs)){
+     nanalogs <- 5
+  }
+  if (!(latdim %in% dimnames) || !(londim %in% dimnames)){
+     stop("'londim' or 'latdim' input doesn't match with 'obs$data' dimension",
+	 " names")
+  }
+  if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){
+     stop("'londim' or 'latdim' input doesn't match with 'exp$data' dimension",
+	 " names")
+  }
+  if (!all(search_obsdims %in% dimnames)){
+     stop("Names in parameter 'search_obsdims' should match 'obs$data' ",
+	  "dimension names.")
+  }
+  if (!all(dim(wt_exp) %in% dim(exp$data))){
+     stop("Dimensions for 'wt_exp' should match 'exp$data' except lat/lon")
+  }
+  if (!all(dim(wt_obs) %in% dim(obs$data))){
+     stop("Dimensions for 'wt_obs' should match 'obs$data' except lat/lon")
+  }
+  plat_exp <- which(dimnamesexp==latdim)
+  plon_exp <- which(dimnamesexp==londim)
+  plat_obs <- which(dimnames==latdim)
+  plon_obs <- which(dimnames==londim)
+  if ((dim(obs$data)[plon_obs]!=dim(exp$data)[plon_exp]) ||
+      (dim(obs$data)[plat_obs]!=dim(exp$data)[plat_exp])){
+     stop("Element 'data' from parameters 'obs' and 'exp' should have",
+      "same lon / lat dimensions if working with regular grids.")
+  }
+  # End of sanity checks; call AdamontAnalog function
+  analog_vals <- AdamontAnalog(exp = exp$data, obs = obs$data, wt_exp = wt_exp,
+			      wt_obs=wt_obs, nanalogs = nanalogs,
+			      method = method, thres = thres,
+			      search_obsdims = search_obsdims, londim = londim,
+			      latdim = latdim )
+
+  return(analog_vals)
+}
+
+#'AdamontAnalog finds analogous data in the reference dataset to experiment data
+#'based on weather types
+#'
+#'@description This function searches for analogs in a reference dataset for 
+#'experiment data, based on corresponding weather types. The experiment data is
+#'typically a hindcast, observations are typically provided by reanalysis data.
+#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+#'
+#'@param exp experiment data, such as \code{$data} array of an object of class 
+#'\code{s2dv_cube}, can be output from quantile correction using 
+#'\code{CST_AdamontQQCorr.R}
+#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but 
+#'lat/lon)
+#'@param obs reference data, can also be \code{$data} array of class 
+#'\code{s2dv_cube}. Note that lat/lon dimensions need to be the same as 
+#'\code{exp}
+#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but *
+#'lat/lon)
+#'@param nanalogs integer defining the number of analog values to return 
+#'(default: 5)
+#'@param method a character string indicating the method used for analog 
+#'definition
+#'       Coded are 'pattcorr': pattern correlation
+#'	           'rain1' (for precip patterns): rain occurrence consistency
+#'	           'rain01' (for precip patterns): rain occurrence/non 
+#'                      occurrence consistency
+#'@param thres real number indicating the threshold to define rain occurrence
+#'/non occurrence in rain(0)1
+#'@param search_obsdims list of dimensions in \code{obs} along which analogs are
+#'searched for
+#'@param londim name of longitude dimension
+#'@param latdim name of latitude dimension
+#'@return analog_vals an object of class \code{s2dv_cube} containing nanalogs 
+#'analog values for each value of \code{exp} input data
+#'@import multiApply
+#'@importFrom s2dverification Subset
+#'@examples
+#'\dontrun{
+#'wt_exp <- sample(1:3, 15*6*3, replace=T)
+#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+#'wt_obs <- sample(1:3, 6*3, replace=T)
+#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+# analog_vals <- AdamontAnalog(exp=lonlat_data$exp$data, 
+#'   obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2)
+#'}
+
+AdamontAnalog <-  function(exp, obs, wt_exp, wt_obs, nanalogs=5, 
+			   method = 'pattcorr', thres = NULL,
+			   search_obsdims = c('member', 'sdate', 'ftime'), 
+			   londim = 'lon', latdim = 'lat') { 
+  # exp: lat, lon, sdate, ftime, member
+  # obs: lat, lon, dims for searching 'sdate' 'ftime'...
+  # wt_exp: sdate, ftime, member 
+  # wt_obs: the dims for searching
+  dimnames <- names(dim(obs))
+  dimnamesexp <- names(dim(exp))
+  if (method %in% c('rain1','rain01') & is.null(thres)){
+     stop("Threshold 'thres' must be defined with methods 'rain1' and 'rain01'")
+  }
+  if (method == 'pattcorr' & !is.null(thres)){
+     warning("Parameter 'thres' is not used with method 'pattcorr'.")
+  }
+  if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){
+     stop("'londim' or 'latdim' input doesn't match with 'exp' dimension names")
+  }
+  # Position of lat/lon dimensions in exp data
+  poslatexp <- which(dimnamesexp == latdim)
+  poslonexp <- which(dimnamesexp == londim)
+  poslatobs <- which(dimnames == latdim)
+  poslonobs <- which(dimnames == londim)
+  if (!all(search_obsdims %in% dimnames)){
+     stop("Names in parameter 'search_obsdims' should match 'obs' ",
+	  "dimension names.")
+  }
+  if (!all(dim(wt_exp) %in% dim(exp))){
+     stop("Dimensions for 'wt_exp' should match 'exp' except lat/lon")
+  }
+  if (!all(dim(wt_obs) %in% dim(obs))){
+     stop("Dimensions for 'wt_obs' should match 'obs' except lat/lon")
+  }
+  if ((dim(obs)[poslonobs]!=dim(exp)[poslonexp]) ||
+      (dim(obs)[poslatobs]!=dim(exp)[poslatexp])){
+     stop("Parameters 'obs' and 'exp' should have same lon / lat dimensions.")
+  }
+
+  ## Reshaping obs:
+  ## The dimensions where to search in a single dim
+  if (length(search_obsdims) > 1) {
+    for (i in 1:(length(search_obsdims) - 1)) {
+      obs <- MergeDims(obs, search_obsdims[i:(i + 1)], 
+		       rename_dim = search_obsdims[i + 1])
+      wt_obs <- MergeDims(wt_obs, search_obsdims[i:(i + 1)], 
+			  rename_dim = search_obsdims[i + 1])
+    }
+  }
+  names(dim(obs))[which(names(dim(obs)) == search_obsdims[length(search_obsdims)])] <- 'time'  
+  names(dim(wt_obs))[which(names(dim(wt_obs)) == search_obsdims[length(search_obsdims)])] <- 'time' 
+  # Split 'time' dim in weather types
+  obs <- SplitDim(obs, split_dim = 'time', indices = as.vector(wt_obs),
+		  new_dim_name='type')
+
+  analog_vals <- Apply(list(exp, obs, wt_exp), 
+		       target_dims = list(c(londim, latdim),
+					  c(londim, latdim, 'time', 'type'), 
+					  NULL),
+		       .analogs, method = method, thres = thres)$output1
+  
+  # Reshaping output:
+  analog_vals <- Subset(analog_vals,along='type',indices=1,drop='selected')
+  poslat <- which(names(dim(analog_vals)) == latdim)
+  poslon <- which(names(dim(analog_vals)) == londim)
+  postime <- which(names(dim(analog_vals)) == 'time') # Dimension with N analogs
+  pos <- 1:length(dim(analog_vals))
+  if (poslatexp > poslonexp){
+      analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],
+					 postime,poslon,poslat))
+  } else {
+      analog_vals <- aperm(analog_vals,c(pos[-c(poslon,poslat,postime)],
+					 postime,poslat,poslon))
+  }
+  # Renaming 'time' dim to 'analog'
+  names(dim(analog_vals))[which(names(dim(analog_vals)) == 'time')] <- 'analog' 
+  return(analog_vals)
+}
+
+
+.analogs <- function(exp, obs, wt_exp, nanalogs = 5, method = 'pattcorr', 
+		     thres = NULL, londimexp = 'lon', latdimexp = 'lat', 
+		     londimobs = 'lon', latdimobs = 'lat') {
+  # exp: lon, lat
+  # obs: lon, lat, time, wt
+  # wt_exp: wt single scalar
+  	
+  search_analog <- switch(method, 'rain1' = .rain1, 'rain01' = .rain01, 
+			  'pattcorr' = .pattcor,
+                          stop(paste0("Adamont Analog function only supports ",
+				      "methods 'rain1', 'rain01', 'pattcorr'")))
+
+  obs <- Subset(obs, along = 'type', indices = wt_exp)
+  accuracy <- Apply(list(exp, obs), target_dims = list(c(londimexp, latdimexp),
+                                                       c(londimobs, latdimobs)),
+                    search_analog, thres = thres)$output1
+  obs <- Subset(obs, along = 'time', indices = order(accuracy, decreasing = TRUE)[1:nanalogs])
+  return(obs)
+}
+
+.rain1 <- function(exp_day, obs_day, thres) {
+  accuracy <- sum((obs_day >= thres) * (exp_day >= thres))
+  return(accuracy)
+}
+.rain01 <- function(exp_day, obs_day, thres) {
+  accuracy <- sum(diag(table((obs_day >= thres),(exp_day >= thres))))
+  return(accuracy)
+}
+.pattcor <- function(exp_day, obs_day, thres = NULL) {
+  accuracy <- cor(as.vector(obs_day),as.vector(exp_day))
+  return(accuracy)
+}
+
+
diff --git a/R/CST_AdamontQQCorr.R b/R/CST_AdamontQQCorr.R
new file mode 100644
index 0000000000000000000000000000000000000000..1d72c2a87b788ded2a1ae8e00d80ff0f5aa0abff
--- /dev/null
+++ b/R/CST_AdamontQQCorr.R
@@ -0,0 +1,393 @@
+#'CST_AdamontQQCorr computes quantile-quantile correction of seasonal or 
+#'decadal forecast data using weather types
+#'
+#'@description This function computes a quantile mapping based on weather types 
+#'for experiment data (typically a hindcast) onto reference \code{obs},
+#'typically provided by reanalysis data.
+#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr}
+#'@author Paola Marson, \email{paola.marson@meteo.fr}
+#'@author Gildas Dayon, \email{gildas.dayon@meteo.fr}
+#'
+#'@param exp experiment data an object of class \code{s2dv_cube}
+#'@param wt_exp corresponding weather types (same dimensions as \code{exp$data}
+#'    but lat/lon)
+#'@param obs reference data, also of class \code{s2dv_cube}. lat/lon dimensions
+#'    can differ from \code{exp} if non rectilinear latlon grids are used, 
+#'    in which case regrid should be set to TRUE and .NearestNeighbors \code{NN}
+#'    output should be provided
+#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but
+#'lat/lon)
+#'@param corrdims list of dimensions in \code{exp} for which quantile mapping
+#'    correction is applied
+#'@param londim character name of longitude dimension in \code{exp} and 
+#'    \code{obs}
+#'@param latdim character name of latitude dimension in \code{exp} and 
+#'    \code{obs}
+#'
+#'@return an object of class \code{s2dv_cube} containing experiment data on the
+#'    lat/lon grid of \code{obs} input data, corrected by quantile mapping 
+#'    depending on the weather types \code{wt_exp}
+#'
+#'@import qmap
+#'@importFrom s2dverification Subset
+#'@import multiApply
+#'@import abind
+#'@examples
+#'\dontrun{
+#'wt_exp <- sample(1:3, 15*6*3, replace=T)
+#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+#'wt_obs <- sample(1:3, 6*3, replace=T)
+#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+#'exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, 
+#'    obs=lonlat_data$obs, wt_obs=wt_obs, 
+#'    corrdims = c('dataset','member','sdate','ftime'))
+#'}
+CST_AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, 
+			      corrdims = c('member','sdate','ftime'),
+			      londim='lon', latdim='lat') {
+
+   if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')){
+     stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ",
+	  "as output by CSTools::CST_Load.")
+   }
+   dimnames <- names(dim(obs$data))
+   dimnamesexp <- names(dim(exp$data))
+   if (!(latdim %in% dimnames) || !(londim %in% dimnames)){
+     stop("'londim' or 'latdim' input doesn't match with 'obs$data' dimension",
+	 " names")
+   }
+   if (!(latdim %in% dimnamesexp) || !(londim %in% dimnamesexp)){
+     stop("'londim' or 'latdim' input doesn't match with 'exp$data' dimension",
+	 " names")
+   }
+   if (!(('time' %in% corrdims) || ('ftime' %in% corrdims))){
+     warning("Forecast time should be one of the dimensions for the correction
+	  specified in corrdims input list")
+   }
+   if (!all(corrdims %in% dimnamesexp)){
+     stop("Names in parameter 'corrdims' should match input dimension names.")
+   }
+   if (!all(dim(wt_exp) %in% dim(exp$data))){
+     stop("Dimensions for 'wt_exp' should match 'exp$data' except lat/lon")
+   }
+   if (!all(dim(wt_obs) %in% dim(obs$data))){
+     stop("Dimensions for 'wt_obs' should match 'obs$data' except lat/lon")
+   }
+   if ((length(dim(exp$lon))==2) || (length(dim(obs$lon))==2)){
+      myNN <- .NearestNeighbors(exp=exp, obs=obs, method='ADA')
+      exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data,
+			       	wt_obs=wt_obs, corrdims=corrdims,
+			       	londim=londim, latdim=latdim, 
+				regrid=TRUE, NN=myNN) 
+   } else {
+   ## If not (standard case)
+      ## exp$data lat/lon dimensions should match obs$data
+      plat_exp <- which(dimnamesexp==latdim)
+      plon_exp <- which(dimnamesexp==londim)
+      plat_obs <- which(dimnames==latdim)
+      plon_obs <- which(dimnames==londim)
+      if ((dim(obs$data)[plon_obs]!=dim(exp$data)[plon_exp]) ||
+	  (dim(obs$data)[plat_obs]!=dim(exp$data)[plat_exp])){
+         stop("Element 'data' from parameters 'obs' and 'exp' should have",
+	      "same lon / lat dimensions if working with regular grids.")
+      }
+      exp_corr <- AdamontQQCorr(exp=exp$data, wt_exp=wt_exp, obs=obs$data,
+			       	wt_obs=wt_obs, corrdims=corrdims,
+			       	londim=londim, latdim=latdim, regrid=FALSE)
+   }
+   return(exp_corr)
+}
+
+
+#'AdamontQQCorr computes quantile-quantile correction of seasonal or decadal
+#'forecast data using weather types
+#'
+#'@description This function computes a quantile mapping based on weather types
+#'for experiment data (typically a hindcast) onto reference \code{obs},
+#'typically provided by reanalysis data.
+#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+#'
+#'@param exp array with named dimensions (such as \code{$data} array of 
+#'experiment data from an object of class \code{s2dv_cube})
+#'@param wt_exp corresponding weather types (same dimensions as \code{exp} but
+#'lat/lon)
+#'@param obs array with named dimensions with reference data (can also be 
+#'\code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ 
+#'from \code{exp} if non rectilinear latlon grids are used, in which case 
+#'regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be
+#'provided
+#'@param wt_obs corresponding weather types (same dimensions as \code{obs} but
+#'lat/lon)
+#'@param corrdims list of dimensions in \code{exp} for which quantile mapping 
+#'correction is applied
+#'@param londim character name of longitude dimension in \code{exp} and 
+#'\code{obs}
+#'@param latdim character name of latitude dimension in \code{exp} and 
+#'\code{obs}
+#'@param regrid (optional) boolean indicating whether .NearestNeighbors
+#'regridding is needed
+#'@param NN (optional, if regrid=TRUE) list (output from .NearestNeighbors)
+#'maps (nlat, nlon) onto (nlat_o, nlon_o)
+#'
+#'@return an array (such as \code{$data} array from an object of class 
+#'\code{s2dv_cube}) with named dimensions, containing experiment data on the 
+#'lat/lon grid of \code{obs} array, corrected by quantile mapping depending on 
+#'the weather types \code{wt_exp}
+#'
+#'@import qmap
+#'@importFrom s2dverification Subset
+#'@import multiApply
+#'@import abind
+#'@examples
+#'\dontrun{
+#'wt_exp <- sample(1:3, 15*6*3, replace=T)
+#'dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+#'wt_obs <- sample(1:3, 6*3, replace=T)
+#'dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+#'exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, 
+#'    obs=lonlat_data$obs$data, wt_obs=wt_obs, 
+#'    corrdims = c('dataset','member','sdate','ftime'))
+#'}
+AdamontQQCorr <- function(exp, wt_exp, obs, wt_obs, 
+			  corrdims = c('member', 'sdate', 'ftime'), 
+			  londim='lon', latdim='lat', regrid=FALSE, NN=NULL) {
+
+   dimnames <- names(dim(obs))
+   dimnamesexp <- names(dim(exp))
+   if (!(latdim %in% dimnames) || !(londim %in% dimnames)){
+     stop("'londim' or 'latdim' input doesn't match with 'obs' dimension names")
+   }
+   if (!(('time' %in% corrdims) || ('ftime' %in% corrdims))){
+     warning("Forecast time should be one of the dimensions for the correction",
+	     " specified in corrdims input list")
+   }
+   if (!all(corrdims %in% dimnamesexp)){
+     stop("Names in parameter 'corrdims' should match input dimension names.")
+   }
+   if (!all(dim(wt_exp) %in% dim(exp))){
+     stop("Dimensions for 'wt_exp' should match 'exp' except lat/lon")
+   }
+   if (!all(dim(wt_obs) %in% dim(obs))){
+     stop("Dimensions for 'wt_obs' should match 'obs' except lat/lon")
+   }
+   if ((regrid == 'TRUE') & is.null(NN)){
+     stop("regrid set to TRUE: provide nearest neighbors input NN")
+   }
+   # The regridding part should only be done if lat/lon dimensions of obs and
+   # exp differ.
+   if (regrid == 'TRUE'){
+      obsdims <- names(dim(obs))
+      poslat <- which(obsdims == latdim)
+      poslon <- which(obsdims == londim)
+      nlat_o <- dim(obs)[poslat]
+      nlon_o <- dim(obs)[poslon]
+      ilat_o <- array(c(1:nlat_o))
+      names(dim(ilat_o))[1] <- latdim
+      ilon_o <- array(c(1:nlon_o))
+      names(dim(ilon_o))[1] <- londim
+      ## First step if obs data is higher resolution than exp data is to use
+      ## nearest neighbor to compute downscaling of exp data	
+      exp_corr <- Apply(list(exp,ilat_o,ilon_o),
+			target_dims=list(c(latdim,londim),latdim,londim),
+			.getNN,NN=NN)$output1
+
+      ## Reorder exp_corr dimensions to match exp dimensions
+      dexpc <- match(names(dim(exp)), names(dim(exp_corr)))
+      exp_corr <- aperm(exp_corr,dexpc)
+      dimnames(exp_corr) <- dimnames(exp)[dexpc]
+      ## Keep original wt_exp for remapping data
+      wt_exp2 <- wt_exp
+      ## Both exp and obs data are now on the same grid
+   } else {
+      ## exp lat/lon dimensions should match obs
+      plat_exp <- which(dimnamesexp==latdim)
+      plon_exp <- which(dimnamesexp==londim)
+      plat_obs <- which(dimnames==latdim)
+      plon_obs <- which(dimnames==londim)
+      if ((dim(obs)[plon_obs]!=dim(exp)[plon_exp]) ||
+	  (dim(obs)[plat_obs]!=dim(exp)[plat_exp])){
+         stop("Parameters 'obs' and 'exp' should have same lon / lat",
+	     " dimensions if regrid set to 'FALSE' (regular grid case).")
+      }
+      exp_corr <- exp
+      ## Keep original wt_exp for remapping data
+      wt_exp2 <- wt_exp
+   }
+   
+   ## Use CST_QuantileMapping function for quantile mapping
+   ## depending on weather type
+   for (i in 1:(length(corrdims) - 1)) {
+      obs <- MergeDims(obs, corrdims[i:(i+1)], rename_dim=corrdims[i+1])
+      wt_obs <- MergeDims(wt_obs, corrdims[i:(i+1)], rename_dim=corrdims[i+1])
+      exp_corr <- MergeDims(exp_corr, corrdims[i:(i+1)], rename_dim=corrdims[i+1])
+      wt_exp2 <- MergeDims(wt_exp2, corrdims[i:(i+1)], rename_dim=corrdims[i+1])
+   }
+   names(dim(obs))[which(names(dim(obs)) == corrdims[length(corrdims)])] <- 'time'
+   names(dim(wt_obs))[which(names(dim(wt_obs)) == corrdims[length(corrdims)])] <- 'time'
+   names(dim(exp_corr))[which(names(dim(exp_corr)) == corrdims[length(corrdims)])] <- 'time'
+   names(dim(wt_exp2))[which(names(dim(wt_exp2)) == corrdims[length(corrdims)])] <- 'time'
+   # Split 'time' dim in weather types
+   obs <- SplitDim(obs, split_dim='time',indices=as.vector(wt_obs),
+		   new_dim_name='type')
+   exp_corr <- SplitDim(exp_corr, split_dim='time',indices=as.vector(wt_exp2),
+		       	new_dim_name='type')
+   ## Add NAs to exp_corr if needed to have compatible sample dimensions
+   numtobs <- dim(obs)[which(names(dim(obs))=='time')]
+   numtexp <- dim(exp_corr)[which(names(dim(exp_corr))=='time')]
+   if (numtexp%%numtobs > 0){
+      ## Create extra dimension and include NAs
+      ndimexp <- names(dim(exp_corr))
+      ndimobs <- names(dim(obs))
+      postime <- which(ndimexp=='time')
+      dimadd <- dim(exp_corr)
+      dimadd[postime] <- ceiling(numtexp/numtobs)*numtobs-numtexp
+      exp_corr <- abind::abind(exp_corr,array(NA,dimadd),along=postime)
+      names(dim(exp_corr)) <- ndimexp
+      exp_corr <- SplitDim(exp_corr,'time',freq=numtobs,indices=NULL)
+      dimobs <- c(dim(obs),1)
+      dim(obs) <- dimobs
+      names(dim(obs)) <- c(ndimobs,'index')
+      res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims=c('time','index'),
+			     method='RQUANT')
+      res <- MergeDims(res,c('time','index'))
+      ## Remove the extra NA values added previously
+      res <- Subset(res,along='time',indices=1:numtexp)
+   } else {
+      ## Apply QuantileMapping to exp_corr depending on weather type
+      res <- QuantileMapping(exp=exp_corr,obs=obs,sample_dims='time',
+			     samplemethod='RQUANT')
+   }
+   rm(exp_corr) # Save space in memory
+   ## Reshape exp_corr data onto time dimension before 'Split'
+   rep_pos <- array(NA,c(time=length(wt_exp2)))
+   pos_time <- which(names(dim(res)) == 'time')
+   pos_type <- which(names(dim(res)) == 'type')
+   for (x in unique(wt_exp2)){
+      rep_pos[which(wt_exp2==x)]<-1:length(which(wt_exp2==x))
+   }
+   exp_corr <- .unsplit_wtype(exp=res,wt_exp=wt_exp2,rep_pos=rep_pos,
+			      pos_time=pos_time) 
+   # Now reshape exp_corr data onto original dimensions
+   dim(exp_corr) <- c(dim(wt_exp), dim(exp_corr)[-c(pos_time,pos_type)])
+   return(exp_corr)
+}
+
+.getNN <- function(exp,ilat,ilon,NN){
+   return(exp[NN$imin_lat[ilat,ilon],NN$imin_lon[ilat,ilon]])
+}
+
+.unsplit_wtype <- function(exp=exp,dim_wt='type',wt_exp=wt_exp,
+			   dim_time='time',rep_pos=rep_pos,pos_time=1){
+   # Initiate output
+   new <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[1]), along=dim_time,
+		 indices=rep_pos[1])
+   dimnames <- names(dim(new))
+   for (x in 2:length(wt_exp)){
+     	dat <- Subset(Subset(exp, along=dim_wt, indices=wt_exp[x]), 
+			   along=dim_time, indices=rep_pos[x])
+        new <- abind::abind(new,dat,along=pos_time)
+   }
+   names(dim(new)) <- dimnames
+   return(new)
+}
+#' ADAMONT Nearest Neighbors computes the distance between reference data grid centroid and SF data grid
+#'
+#'@author Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+#'@author Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+#'@description This function computes the nearest neighbor for each reference data (lon, lat) point in the experiment dataset by computing the distance between the reference dataset grid and the experiment data. This is the first step in the ADAMONT method adapted from Verfaillie et al. (2018).
+#'
+#'@param method a string among three options ('ADA': standard ADAMONT distance, 'simple': lon/lat straight euclidian distance, 'radius': distance on the sphere)
+#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment longitudes in \code{$lon} and latitudes in \code{$lat}
+#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the reference data on a different grid, with longitudes in \code{$lon} and latitudes in \code{$lat}.
+#'
+#'@return NN a list, containing the following:
+#'          min_lon: array of dimensions \code{obs$lon} giving the longitude of closest gridpoint in exp
+#'          min_lat: array of dimensions \code{obs$lat} giving the latitude of closest gridpoint in exp
+#'          imin_lon: array of dimensions \code{obs$lon} giving the longitude index of closest gridpoint in exp
+#'          imin_lat: array of dimensions \code{obs$lat} giving the latitude index of closest gridpoint in exp
+#'
+#'@importFrom s2dverification Subset
+#'@import ncdf4
+#'@noRd
+.NearestNeighbors <- function (exp, obs, method='ADA') {
+   
+   if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) {
+     stop("Inputs 'exp' and 'obs' must be of class 's2dv_cube', ",
+          "as output by CSTools::CST_Load.")
+   }
+   exp_lon <- exp$lon
+   exp_lat <- exp$lat
+   obs_lon <- obs$lon
+   obs_lat <- obs$lat
+   dim_exp_lon <- dim(exp_lon)
+   dim_exp_lat <- dim(exp_lat)
+   dim_obs_lon <- dim(obs_lon)
+   dim_obs_lat <- dim(obs_lat)
+   # Check if one of the grids is non-regular:
+   if ((length(dim_exp_lon)==2) || (length(dim_obs_lon)==2)){
+   # Flatten longitudes and latitudes in case of 2-D longitudes and latitudes (Lambert grids, etc.)
+      if ((length(dim_exp_lon)==2) & (length(dim_exp_lat)==2)){
+         dim(exp_lon) <- c(dim_exp_lon[1]*dim_exp_lon[2])
+         dim(exp_lat) <- c(dim_exp_lat[1]*dim_exp_lat[2])
+      }
+      if ((length(dim_obs_lon)==2) & (length(dim_obs_lat)==2)){
+         dim(obs_lon) <- c(dim_obs_lon[1]*dim_obs_lon[2])
+         dim(obs_lat) <- c(dim_obs_lat[1]*dim_obs_lat[2])
+      }
+      # Now lat and lon arrays have 1 dimension, length npt (= nlat*nlon)
+      OBS_grid <- cbind(obs_lon,obs_lat)
+      EXP_grid <- cbind(exp_lon,exp_lat)
+      dist_min<-min_lon<-min_lat<-imin_lon<-imin_lat<-array(dim=nrow(OBS_grid))
+      if (method == 'ADA'){
+   	C<-cos(OBS_grid[,2]*pi/180)^2
+   	for (i in 1:nrow(OBS_grid)){
+   		dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+C[i]*(OBS_grid[i,1]-EXP_grid[,1])^2
+   		dist_min[i]<-min(dist)
+   		min_lon[i]<-EXP_grid[which.min(dist),1]
+   		min_lat[i]<-EXP_grid[which.min(dist),2]
+   		imin_lon[i]<-which(exp_lon==min_lon[i])
+   		imin_lat[i]<-which(exp_lat==min_lat[i])
+   	}
+      } else if (method == 'simple'){
+   	for (i in 1:nrow(OBS_grid)){
+   		dist<-(OBS_grid[i,2]-EXP_grid[,2])^2+(OBS_grid[i,1]-EXP_grid[,1])^2
+   		dist_min[i]<-min(dist)
+   		min_lon[i]<-EXP_grid[which.min(dist),1]
+   		min_lat[i]<-EXP_grid[which.min(dist),2]
+   		imin_lon[i]<-which(exp_lon==min_lon[i]) 
+   		imin_lat[i]<-which(exp_lat==min_lat[i])
+   	}
+      } else if (method == 'radius'){
+   	R <- 6371e3 # metres, Earth radius
+   	EXP_gridr<-EXP_grid*pi/180
+   	OBS_gridr<-OBS_grid*pi/180
+   	for (i in 1:nrow(OBS_grid)){
+   		a<-sin((OBS_gridr[i,2]-EXP_gridr[,2])/2)^2 + cos(OBS_gridr[i,2])*cos(EXP_gridr[,2])*sin((OBS_gridr[i,1]-EXP_gridr[,1])/2)^2
+   		c<-2*atan2(sqrt(a),sqrt(1-a))
+   		dist<-R*c
+   		dist_min[i]<-min(dist)
+   		min_lon[i]<-EXP_grid[which.min(dist),1]
+   		min_lat[i]<-EXP_grid[which.min(dist),2]
+   		imin_lon[i]<-which(exp_lon==min_lon[i])
+   		imin_lat[i]<-which(exp_lat==min_lat[i])
+   	}
+      } else {
+      	stop("AdamontNearestNeighbors supports method = 'ADA', 'simple' or 'radius' only.")
+      }
+   
+      # Reshape outputs to original grid
+      dim(min_lon)=dim_obs_lon
+      dim(min_lat)=dim_obs_lat
+      dim(imin_lon)=dim_obs_lon
+      dim(imin_lat)=dim_obs_lat
+
+   } else {
+   # Regular lon/lat grid case: has been handled by CST_Load()
+     stop("AdamontNearestNeighbors is meant for non-regular lat/lon grids; use e.g. CST_Load to interpolate exp onto obs grid")
+   }
+
+   NN=list(min_lon=min_lon, min_lat=min_lat, imin_lon=imin_lon, imin_lat=imin_lat)
+
+   return(NN)
+}
diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R
index deb60760a2fe294b95770bbd6824fbd026e09b16..5e85182d6d23d36b1afd60e204304bf14410bfe2 100644
--- a/R/CST_SplitDim.R
+++ b/R/CST_SplitDim.R
@@ -258,3 +258,4 @@ rebuild <- function(x, data, along, indices, max_times) {
     }
     return(a)
 }
+
diff --git a/man/AdamontAnalog.Rd b/man/AdamontAnalog.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..728ed9bccf2ab7eb6c74b37cbc5bef377817bbd8
--- /dev/null
+++ b/man/AdamontAnalog.Rd
@@ -0,0 +1,78 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_AdamontAnalog.R
+\name{AdamontAnalog}
+\alias{AdamontAnalog}
+\title{AdamontAnalog finds analogous data in the reference dataset to experiment data
+based on weather types}
+\usage{
+AdamontAnalog(
+  exp,
+  obs,
+  wt_exp,
+  wt_obs,
+  nanalogs = 5,
+  method = "pattcorr",
+  thres = NULL,
+  search_obsdims = c("member", "sdate", "ftime"),
+  londim = "lon",
+  latdim = "lat"
+)
+}
+\arguments{
+\item{exp}{experiment data, such as \code{$data} array of an object of class 
+\code{s2dv_cube}, can be output from quantile correction using 
+\code{CST_AdamontQQCorr.R}}
+
+\item{obs}{reference data, can also be \code{$data} array of class 
+\code{s2dv_cube}. Note that lat/lon dimensions need to be the same as 
+\code{exp}}
+
+\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but 
+lat/lon)}
+
+\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but *
+lat/lon)}
+
+\item{nanalogs}{integer defining the number of analog values to return 
+(default: 5)}
+
+\item{method}{a character string indicating the method used for analog 
+definition
+      Coded are 'pattcorr': pattern correlation
+           'rain1' (for precip patterns): rain occurrence consistency
+           'rain01' (for precip patterns): rain occurrence/non 
+                     occurrence consistency}
+
+\item{thres}{real number indicating the threshold to define rain occurrence
+/non occurrence in rain(0)1}
+
+\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are
+searched for}
+
+\item{londim}{name of longitude dimension}
+
+\item{latdim}{name of latitude dimension}
+}
+\value{
+analog_vals an object of class \code{s2dv_cube} containing nanalogs 
+analog values for each value of \code{exp} input data
+}
+\description{
+This function searches for analogs in a reference dataset for 
+experiment data, based on corresponding weather types. The experiment data is
+typically a hindcast, observations are typically provided by reanalysis data.
+}
+\examples{
+\dontrun{
+wt_exp <- sample(1:3, 15*6*3, replace=T)
+dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+wt_obs <- sample(1:3, 6*3, replace=T)
+dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+  obs=lonlat_data$obs$data, wt_exp=wt_exp, wt_obs=wt_obs, nanalogs=2)
+}
+}
+\author{
+Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+
+Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+}
diff --git a/man/AdamontQQCorr.Rd b/man/AdamontQQCorr.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..ec49bad348628476df668ad136d1770197f2f430
--- /dev/null
+++ b/man/AdamontQQCorr.Rd
@@ -0,0 +1,77 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_AdamontQQCorr.R
+\name{AdamontQQCorr}
+\alias{AdamontQQCorr}
+\title{AdamontQQCorr computes quantile-quantile correction of seasonal or decadal
+forecast data using weather types}
+\usage{
+AdamontQQCorr(
+  exp,
+  wt_exp,
+  obs,
+  wt_obs,
+  corrdims = c("member", "sdate", "ftime"),
+  londim = "lon",
+  latdim = "lat",
+  regrid = FALSE,
+  NN = NULL
+)
+}
+\arguments{
+\item{exp}{array with named dimensions (such as \code{$data} array of 
+experiment data from an object of class \code{s2dv_cube})}
+
+\item{wt_exp}{corresponding weather types (same dimensions as \code{exp} but
+lat/lon)}
+
+\item{obs}{array with named dimensions with reference data (can also be 
+\code{$data} array of class \code{s2dv_cube}). lat/lon dimensions can differ 
+from \code{exp} if non rectilinear latlon grids are used, in which case 
+regrid should be set to TRUE and .NearestNeighbors \code{NN} output should be
+provided}
+
+\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but
+lat/lon)}
+
+\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping 
+correction is applied}
+
+\item{londim}{character name of longitude dimension in \code{exp} and 
+\code{obs}}
+
+\item{latdim}{character name of latitude dimension in \code{exp} and 
+\code{obs}}
+
+\item{regrid}{(optional) boolean indicating whether .NearestNeighbors
+regridding is needed}
+
+\item{NN}{(optional, if regrid=TRUE) list (output from .NearestNeighbors)
+maps (nlat, nlon) onto (nlat_o, nlon_o)}
+}
+\value{
+an array (such as \code{$data} array from an object of class 
+\code{s2dv_cube}) with named dimensions, containing experiment data on the 
+lat/lon grid of \code{obs} array, corrected by quantile mapping depending on 
+the weather types \code{wt_exp}
+}
+\description{
+This function computes a quantile mapping based on weather types
+for experiment data (typically a hindcast) onto reference \code{obs},
+typically provided by reanalysis data.
+}
+\examples{
+\dontrun{
+wt_exp <- sample(1:3, 15*6*3, replace=T)
+dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+wt_obs <- sample(1:3, 6*3, replace=T)
+dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+exp_corr <- AdamontQQCorr(exp=lonlat_data$exp$data, wt_exp=wt_exp, 
+   obs=lonlat_data$obs$data, wt_obs=wt_obs, 
+   corrdims = c('dataset','member','sdate','ftime'))
+}
+}
+\author{
+Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+
+Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+}
diff --git a/man/CST_AdamontAnalog.Rd b/man/CST_AdamontAnalog.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..7d10ae02e62255e74ee3e17553f5a6da716b8a54
--- /dev/null
+++ b/man/CST_AdamontAnalog.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_AdamontAnalog.R
+\name{CST_AdamontAnalog}
+\alias{CST_AdamontAnalog}
+\title{CST_AdamontAnalog finds analogous data in the reference dataset to experiment
+data based on weather types}
+\usage{
+CST_AdamontAnalog(
+  exp,
+  obs,
+  wt_exp,
+  wt_obs,
+  nanalogs,
+  method = "pattcorr",
+  thres = NULL,
+  search_obsdims = c("member", "sdate", "ftime"),
+  londim = "lon",
+  latdim = "lat"
+)
+}
+\arguments{
+\item{exp}{experiment data an object of class \code{s2dv_cube}, can be output 
+from quantile correction using CST_AdamontQQCorr.R}
+
+\item{obs}{reference data, also of class \code{s2dv_cube}. Note that lat/lon 
+dimensions need to be the same as \code{exp}}
+
+\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data} 
+but lat/lon)}
+
+\item{wt_obs}{corresponding weather types (same dimensions as \code{obs$data}
+but lat/lon)}
+
+\item{nanalogs}{integer defining the number of analog values to return
+(default: 5)}
+
+\item{method}{a character string indicating the method used for analog
+definition
+      Coded are 'pattcorr': pattern correlation
+           'rain1' (for precip patterns): rain occurrence consistency
+           'rain01' (for precip patterns): rain occurrence/non 
+                       occurrence consistency}
+
+\item{thres}{real number indicating the threshold to define rain 
+occurrence/non occurrence in rain(0)1}
+
+\item{search_obsdims}{list of dimensions in \code{obs} along which analogs are
+searched for}
+
+\item{londim}{name of longitude dimension}
+
+\item{latdim}{name of latitude dimension}
+}
+\value{
+analog_vals an object of class \code{s2dv_cube} containing nanalogs
+analog values for each value of \code{exp} input data
+}
+\description{
+This function searches for analogs in a reference dataset for 
+experiment data, based on corresponding weather types. The experiment data is
+typically a hindcast, observations are typically provided by reanalysis data.
+}
+\examples{
+\dontrun{
+wt_exp <- sample(1:3, 15*6*3, replace=T)
+dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+wt_obs <- sample(1:3, 6*3, replace=T)
+dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+}
+}
+\author{
+Paola Marson, \email{paola.marson@meteo.fr} for PROSNOW version
+
+Lauriane Batté, \email{lauriane.batte@meteo.fr} for CSTools adaptation
+}
diff --git a/man/CST_AdamontQQCorr.Rd b/man/CST_AdamontQQCorr.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..967ce1710e4ddccdfaed41ab378b2775d331782d
--- /dev/null
+++ b/man/CST_AdamontQQCorr.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_AdamontQQCorr.R
+\name{CST_AdamontQQCorr}
+\alias{CST_AdamontQQCorr}
+\title{CST_AdamontQQCorr computes quantile-quantile correction of seasonal or 
+decadal forecast data using weather types}
+\usage{
+CST_AdamontQQCorr(
+  exp,
+  wt_exp,
+  obs,
+  wt_obs,
+  corrdims = c("member", "sdate", "ftime"),
+  londim = "lon",
+  latdim = "lat"
+)
+}
+\arguments{
+\item{exp}{experiment data an object of class \code{s2dv_cube}}
+
+\item{wt_exp}{corresponding weather types (same dimensions as \code{exp$data}
+but lat/lon)}
+
+\item{obs}{reference data, also of class \code{s2dv_cube}. lat/lon dimensions
+can differ from \code{exp} if non rectilinear latlon grids are used, 
+in which case regrid should be set to TRUE and .NearestNeighbors \code{NN}
+output should be provided}
+
+\item{wt_obs}{corresponding weather types (same dimensions as \code{obs} but
+lat/lon)}
+
+\item{corrdims}{list of dimensions in \code{exp} for which quantile mapping
+correction is applied}
+
+\item{londim}{character name of longitude dimension in \code{exp} and 
+\code{obs}}
+
+\item{latdim}{character name of latitude dimension in \code{exp} and 
+\code{obs}}
+}
+\value{
+an object of class \code{s2dv_cube} containing experiment data on the
+   lat/lon grid of \code{obs} input data, corrected by quantile mapping 
+   depending on the weather types \code{wt_exp}
+}
+\description{
+This function computes a quantile mapping based on weather types 
+for experiment data (typically a hindcast) onto reference \code{obs},
+typically provided by reanalysis data.
+}
+\examples{
+\dontrun{
+wt_exp <- sample(1:3, 15*6*3, replace=T)
+dim(wt_exp) <- c(dataset=1, member=15, sdate=6, ftime=3)
+wt_obs <- sample(1:3, 6*3, replace=T)
+dim(wt_obs) <- c(dataset=1, member=1, sdate=6, ftime=3)
+exp_corr <- CST_AdamontQQCorr(exp=lonlat_data$exp, wt_exp=wt_exp, 
+   obs=lonlat_data$obs, wt_obs=wt_obs, 
+   corrdims = c('dataset','member','sdate','ftime'))
+}
+}
+\author{
+Lauriane Batté, \email{lauriane.batte@meteo.fr}
+
+Paola Marson, \email{paola.marson@meteo.fr}
+
+Gildas Dayon, \email{gildas.dayon@meteo.fr}
+}
diff --git a/src/CSTools.so b/src/CSTools.so
index 013c91ca1300cce79ee5a744faa387225f9fd38f..b5e81187254aa25073268fa773029085c71c0779 100755
Binary files a/src/CSTools.so and b/src/CSTools.so differ