diff --git a/.Rbuildignore b/.Rbuildignore
index 6a3eead..ad5ee0d 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -2,3 +2,6 @@
 ^README\.md$
 ^.github$
 ^data-raw
+^tests/last.dump.rda
+^tests/.clustersize
+^tests/.proc
diff --git a/ChangeLog b/ChangeLog
index dc6da4a..61d748a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,4 @@
-5.2-0.9xxx (09/26/2024)
+5.2-0.9xxx (10/31/2024)
 -----
 Added dataset include_2024.
 
@@ -10,11 +10,14 @@ Fixed bug in setting time index when imputation is present in e0.predict.subnat(
 
 Added argument use.wpp.data to run.e0.mcmc(). 
 
-Added the mean variant into sumary output.
+Added the mean variant into summary output.
 
 e0.trajectories.plot() got arguments to plot the means (show.mean), the medians (show.median),
 as well as selected trajectories (traj.index).
 
+Added option of keeping multiple predictions directories in one simulation directory
+(argument "subdir" added to various functions).
+
 5.2-0 (09/15/2023)
 -----
 Annual subnational projections are now possible, via the argument "annual"
diff --git a/DESCRIPTION b/DESCRIPTION
index e7a08bc..7e93381 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,10 +1,22 @@
 Package: bayesLife
 Type: Package
 Title: Bayesian Projection of Life Expectancy
-Version: 5.2-0.9011
-Date: 2024-09-26
+Version: 5.2-0.9012
+Date: 2024-10-31
 Author: Hana Sevcikova, Adrian Raftery, Jennifer Chunn 
 Maintainer: Hana Sevcikova <hanas@uw.edu>
+Authors@R: c(person(given = "Hana",
+                    family = "Sevcikova",
+                    role = c("cre", "aut"),
+                    email = "hanas@uw.edu"),
+             person(given = "Adrian",
+                      family = "Raftery",
+                      role = "aut",
+                      email = "raftery@uw.edu"),
+             person(given = "Jennifer",
+                      family = "Chunn",
+                      role = "aut")
+            )
 Description: Making probabilistic projections of life expectancy for all countries of the world, using a Bayesian hierarchical model <doi:10.1007/s13524-012-0193-x>. Subnational projections are also supported.
 Depends: 
     bayesTFR (>= 7.3-0),
diff --git a/NAMESPACE b/NAMESPACE
index ca2a8f9..848f7e0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,6 +36,7 @@ export(
 	create.thinned.e0.mcmc,
 	get.e0.prediction,
 	has.e0.prediction,
+	available.e0.predictions,
 	get.e0.jmale.prediction,
 	has.e0.jmale.prediction,
 	get.rege0.prediction,
diff --git a/R/project_subnat.R b/R/project_subnat.R
index a6091c6..20d547d 100644
--- a/R/project_subnat.R
+++ b/R/project_subnat.R
@@ -1,6 +1,6 @@
 e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(), 'bayesLife.output'),
                               method = c("ar1", "shift", "scale"), predict.jmale = FALSE, my.e0M.file = NULL,
-                               end.year=2100, start.year=NULL, output.dir = NULL, annual = NULL,
+                               end.year=2100, start.year=NULL, subdir = "predictions", output.dir = NULL, annual = NULL,
                               nr.traj=NULL, seed = NULL, ar.pars = NULL, 
                                save.as.ascii = 0, verbose = TRUE, jmale.estimates = NULL, ...) {
     # Run subnational projections, using the Scale AR(1) model applied to a national bayesLife simulation 
@@ -30,7 +30,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(),
         return(compute.alpha.ar1(...))
     
     method <- match.arg(method)
-    wpred <- get.e0.prediction(sim.dir) # contains national projections
+    wpred <- get.e0.prediction(sim.dir, subdir = subdir) # contains national projections
     wdata <- wpred$e0.matrix.reconstructed
     wmeta <- wpred$mcmc.set$meta
     if(!is.null(seed)) set.seed(seed)
@@ -153,7 +153,7 @@ e0.predict.subnat <- function(countries, my.e0.file, sim.dir=file.path(getwd(),
         PIs_cqp <- array(NA, c(nr.reg, length(quantiles.to.keep), nrow(wtrajs)),
                      dimnames=list(meta$regions$country_code, dimnames(wpred$quantiles)[[2]], dimnames(wtrajs)[[1]]))
         mean_sd <- array(NA, c(nr.reg, 2, nrow(wtrajs)))
-        meta$Tc.index <- .get.Tcindex(meta$e0.matrix, cnames = meta$regions$country_name)
+        meta$Tc.index <- .get.Tcindex(meta$e0.matrix, cnames = meta$regions$country_name, stop.if.less.than2 = FALSE) # allow just one data point
         country.char <- as.character(country.obj$code)
         e0reconstructed <- meta$e0.matrix
         
diff --git a/R/projection_fcns.R b/R/projection_fcns.R
index d39f9f5..431378e 100644
--- a/R/projection_fcns.R
+++ b/R/projection_fcns.R
@@ -21,7 +21,8 @@ e0.predict <- function(mcmc.set = NULL, end.year = 2100,
                        replace.output = FALSE, predict.jmale = TRUE, 
                        nr.traj = NULL, thin = NULL, burnin = 10000, 
                        use.diagnostics = FALSE, save.as.ascii = 0, start.year = NULL,
-                       output.dir = NULL, low.memory = TRUE, ignore.last.observed = FALSE,
+                       output.dir = NULL, subdir = "predictions", low.memory = TRUE, 
+                       ignore.last.observed = FALSE,
                        seed = NULL, verbose = TRUE, ...){
 	if(!is.null(mcmc.set)) {
 		if (!inherits(mcmc.set, 'bayesLife.mcmc.set')) {
@@ -41,7 +42,7 @@ e0.predict <- function(mcmc.set = NULL, end.year = 2100,
 					replace.output = replace.output,  
 					nr.traj = nr.traj, thin = thin, burnin = burnin, 
 					save.as.ascii = save.as.ascii, start.year = start.year,
-					output.dir = output.dir, ignore.last.observed = ignore.last.observed,
+					output.dir = output.dir, subdir = subdir, ignore.last.observed = ignore.last.observed,
 					verbose = verbose)
 	if(predict.jmale && mcmc.set$meta$sex == 'F')
 		pred <- e0.jmale.predict(pred, ..., save.as.ascii = save.as.ascii, verbose = verbose)
@@ -49,7 +50,7 @@ e0.predict <- function(mcmc.set = NULL, end.year = 2100,
 }
 
 e0.predict.extra <- function(sim.dir = file.path(getwd(), 'bayesLife.output'), 
-					prediction.dir = sim.dir, 
+					prediction.dir = sim.dir, subdir = "predictions", 
 					countries = NULL, save.as.ascii = 1000, verbose = TRUE, ...) {
 	# Run prediction for given countries/regions (as codes). If they are not given it will be set to countries 
 	# for which there are MCMC results but no prediction.
@@ -58,9 +59,9 @@ e0.predict.extra <- function(sim.dir = file.path(getwd(), 'bayesLife.output'),
 	mcmc.set <- get.e0.mcmc(sim.dir)
 	if(is.null(mcmc.set))
 		stop('Error in "sim.dir" argument.')
-	pred <- get.e0.prediction(sim.dir=prediction.dir)
+	pred <- get.e0.prediction(sim.dir=prediction.dir, subdir = subdir)
 	if(is.null(pred))
-		stop('Error in "prediction.dir" argument.')
+		stop('Error in "sim.dir", "prediction.dir" or/and "subdir" argument. Use available.e0.predictions() to check on valid predictions directories.')
 	if(length(setdiff(pred$mcmc.set$meta$regions$country_code, mcmc.set$meta$regions$country_code)) > 0)
 		stop('Prediction is inconsistent with the mcmc results. Use e0.predict.')
 	if(is.null(countries)) {
@@ -77,7 +78,7 @@ e0.predict.extra <- function(sim.dir = file.path(getwd(), 'bayesLife.output'),
 	new.pred <- make.e0.prediction(mcmc.set, start.year=pred$start.year, end.year=pred$end.year, replace.output=FALSE,
 									nr.traj=pred$nr.traj, burnin=pred$burnin,
 									countries=countries.idx, save.as.ascii=0, output.dir=prediction.dir,
-									force.creating.thinned.mcmc=TRUE,
+									subdir = subdir, force.creating.thinned.mcmc=TRUE,
 									write.summary.files=FALSE, ignore.last.observed = pred$ignore.last.observed,
 									verbose=verbose)
 									
@@ -121,7 +122,7 @@ e0.predict.extra <- function(sim.dir = file.path(getwd(), 'bayesLife.output'),
 
 e0.prediction.setup <- function(...) {
     setup <- list(...)
-    mcmc.set <- start.year <- end.year <- burnin <- replace.output <- verbose <- countries <- NULL # to avoid R check note "no visible binding ..."
+    mcmc.set <- start.year <- end.year <- burnin <- replace.output <- verbose <- countries <- subdir <- NULL # to avoid R check note "no visible binding ..."
     if(is.null(setup$thin)) setup$thin <- NA # this is because thin is a method in coda and the naming clashes within the setup expression
     setup <- within(setup, {
         meta <- mcmc.set$meta
@@ -151,10 +152,10 @@ e0.prediction.setup <- function(...) {
     
         #setup output directory
         if (!exists("output.dir") || is.null(output.dir)) output.dir <- meta$output.dir
-        outdir <- file.path(output.dir, 'predictions')
+        outdir <- file.path(output.dir, basename(subdir))
     
         if(is.null(get0("countries"))) {
-            if(!replace.output && has.e0.prediction(sim.dir = output.dir))
+            if(!replace.output && has.e0.prediction(sim.dir = output.dir, subdir = subdir))
                 stop('Prediction in ', outdir,
                     ' already exists.\nSet replace.output=TRUE if you want to overwrite existing projections.')
             unlink(outdir, recursive=TRUE)
@@ -359,11 +360,11 @@ get.friendly.variant.names.bayesLife.prediction <- function(pred, ...)
 	return(c('median', 'lower 80', 'upper 80', 'lower 95', 'upper 95', 'mean', 'constant'))	
 
 convert.e0.trajectories <- function(dir=file.path(getwd(), 'bayesLife.output'), 
-								 n=1000, output.dir=NULL, 
+								 n=1000, subdir = "predictions", output.dir=NULL, 
 								 verbose=FALSE) {
 	# Converts all trajectory rda files into UN ascii, selecting n trajectories by equal spacing.
 	if(n <= 0) return()
-	pred <- get.e0.prediction(sim.dir=dir)
+	pred <- get.e0.prediction(sim.dir=dir, subdir = subdir)
 	predsex <- pred$mcmc.set$meta$sex
 	preds <- list()
 	preds[[predsex]] <- pred
@@ -376,15 +377,15 @@ convert.e0.trajectories <- function(dir=file.path(getwd(), 'bayesLife.output'),
 			else outdir <- output.dir
 		}
 		if(!file.exists(outdir)) dir.create(outdir, recursive=TRUE)
-		cat('Converting ', list(M='Male', F='Female')[[sex]], ' trajectories from', dir, '\n')
+		cat('Converting ', list(M='Male', F='Female')[[sex]], ' trajectories from', file.path(dir, subdir), '\n')
 		bayesTFR:::do.convert.trajectories(pred=preds[[sex]], n=n, output.dir=outdir, verbose=verbose)
 	}
 }
 
 write.e0.projection.summary <- function(dir=file.path(getwd(), 'bayesLife.output'), 
-									 output.dir=NULL, revision=NULL, adjusted=FALSE) {
+                                        subdir = "predictions", output.dir=NULL, revision=NULL, adjusted=FALSE) {
 # Writes four prediction summary files, one in a user-friendly format, one in a UN-format, one for each sex.
-	pred <- get.e0.prediction(sim.dir=dir)
+	pred <- get.e0.prediction(sim.dir=dir, subdir = subdir)
 	predsex <- pred$mcmc.set$meta$sex
 	preds <- list()
 	preds[[predsex]] <- pred
@@ -429,40 +430,41 @@ get.e0.reconstructed <- function(data, meta) {
     } else store.bayesLife.prediction(new.pred)
 }
 
-e0.median.reset <- function(sim.dir, countries = NULL, joint.male=FALSE) {
+e0.median.reset <- function(sim.dir, countries = NULL, joint.male=FALSE, ...) {
     if(is.null(countries)) {
-        pred <- get.e0.prediction(sim.dir, joint.male = joint.male)
+        pred <- get.e0.prediction(sim.dir, joint.male = joint.male, ...)
         pred$median.shift <- NULL
         .e0.store.adjustment(pred, sim.dir, joint.male)
         cat('\nMedians for all countries reset.\n')
     } else
-	    for(country in countries) pred <- e0.median.shift(sim.dir, country, reset=TRUE, joint.male=joint.male)
+	    for(country in countries) pred <- e0.median.shift(sim.dir, country, reset=TRUE, joint.male=joint.male, ...)
 	invisible(pred)
 }
 
 get.e0.shift <- function(country.code, pred) return(bayesTFR::get.tfr.shift(country.code, pred))
 
-e0.median.shift <- function(sim.dir, country, reset=FALSE, shift=0, from=NULL, to=NULL, joint.male=FALSE) {
-	pred <- get.e0.prediction(sim.dir, joint.male=joint.male)
+e0.median.shift <- function(sim.dir, country, reset=FALSE, shift=0, from=NULL, to=NULL, joint.male=FALSE, ...) {
+	pred <- get.e0.prediction(sim.dir, joint.male=joint.male, ...)
 	new.pred <- bayesTFR:::.bdem.median.shift(pred, type='e0', country=country, reset=reset, 
 				shift=shift, from=from, to=to)
 	.e0.store.adjustment(new.pred, sim.dir, joint.male)
 	invisible(new.pred)
 }
 
-e0.median.set <- function(sim.dir, country, values, years=NULL, joint.male=FALSE) {
-	pred <- get.e0.prediction(sim.dir, joint.male=joint.male)
+e0.median.set <- function(sim.dir, country, values, years=NULL, joint.male=FALSE, ...) {
+	pred <- get.e0.prediction(sim.dir, joint.male=joint.male, ...)
 	new.pred <- bayesTFR:::.bdem.median.set(pred, type='e0', country=country, 
 								values=values, years=years)
 	.e0.store.adjustment(new.pred, sim.dir, joint.male)
 	invisible(new.pred)
 }
 
-e0.median.adjust.jmale <- function(sim.dir, countries, factors = c(1.2, 1.1)) {
-    pred <- get.e0.prediction(sim.dir)
-    if (is.null(pred)) stop('No valid prediction in ', sim.dir)
+e0.median.adjust.jmale <- function(sim.dir, countries, factors = c(1.2, 1.1), subdir = "predictions") {
+    pred <- get.e0.prediction(sim.dir, subdir = subdir)
+    if (is.null(pred)) stop('Prediction not found in ', file.path(sim.dir, subdir), 
+                            '. Check available.e0.predictions() and use the subdir argument to set non-standard prediction subdirectory.')
     joint.male <- get.e0.jmale.prediction(pred)
-    if (is.null(joint.male)) stop('No valid male prediction in ', sim.dir)
+    if (is.null(joint.male)) stop('No valid male prediction in ', file.path(sim.dir, subdir))
     mcmc.set <- pred$mcmc.set
     if(is.null(countries)) {
         cat('\nNo countries given. Nothing to be done.\n')
@@ -485,11 +487,11 @@ e0.median.adjust.jmale <- function(sim.dir, countries, factors = c(1.2, 1.1)) {
                       new.meds[get.country.object(countries[icountry], mcmc.set$meta)$index,], joint.male = TRUE)
     }
     # reload adjusted prediction
-    invisible(get.e0.prediction(sim.dir, joint.male = TRUE))
+    invisible(get.e0.prediction(sim.dir, joint.male = TRUE, subdir = subdir))
 }
 
-e0.shift.prediction.to.wpp <- function(sim.dir, joint.male = FALSE, ...){
-    pred <- get.e0.prediction(sim.dir, joint.male = joint.male)
+e0.shift.prediction.to.wpp <- function(sim.dir, joint.male = FALSE, subdir = "predictions", ...){
+    pred <- get.e0.prediction(sim.dir, joint.male = joint.male, subdir = subdir)
     new.pred <- bayesTFR:::.do.shift.prediction.to.wpp(pred, 
                                                        wpp.dataset = if(joint.male) "e0Mproj" else "e0Fproj", ...)
     .e0.store.adjustment(new.pred, sim.dir, joint.male)
diff --git a/man/convert.trajectories.Rd b/man/convert.trajectories.Rd
index 78d69aa..3fc82ec 100644
--- a/man/convert.trajectories.Rd
+++ b/man/convert.trajectories.Rd
@@ -8,11 +8,12 @@ Converts trajectories of the life expectancy stored in a binary format into two
 }
 \usage{
 convert.e0.trajectories(dir = file.path(getwd(), "bayesLife.output"), 
-    n = 1000, output.dir = NULL, verbose = FALSE)
+    n = 1000, subdir = "predictions", output.dir = NULL, verbose = FALSE)
 }
 \arguments{
   \item{dir}{Directory containing the prediction object. It should correspond to the \code{output.dir} argument of the \code{\link{e0.predict}} function.}
   \item{n}{Number of trajectories to be stored. It can be either a single number or the word \dQuote{all} in which case all available trajectories are converted.}
+  \item{subdir}{Name of subdirectory of \code{dir} containing the prediction.}
   \item{output.dir}{Directory in which the resulting files will be stored. If \code{NULL} the same directory is used as for the prediction. Otherwise, if the directory contains joint predictions for both sexes, the ouptuts are stored into subdirectories \file{F} and \file{M}.}
   \item{verbose}{Logical switching log messages on and off.}
 }
diff --git a/man/e0.median.set.Rd b/man/e0.median.set.Rd
index 2e380d1..cf7c873 100644
--- a/man/e0.median.set.Rd
+++ b/man/e0.median.set.Rd
@@ -12,16 +12,18 @@ Editing the Projection Medians
 These functions are to be used by expert analysts. They allow to change the projection medians either to specific values, including the WPP values, or shift the medians by a given constant or a factor.
 }
 \usage{
-e0.median.set(sim.dir, country, values, years = NULL, joint.male = FALSE)
+e0.median.set(sim.dir, country, values, years = NULL, joint.male = FALSE, \dots)
 
 e0.median.shift(sim.dir, country, reset = FALSE, shift = 0, 
-    from = NULL, to = NULL, joint.male = FALSE)
+    from = NULL, to = NULL, joint.male = FALSE, \dots)
     
-e0.median.adjust.jmale(sim.dir, countries, factors = c(1.2, 1.1))
+e0.median.adjust.jmale(sim.dir, countries, factors = c(1.2, 1.1),
+    subdir = "predictions")
 
-e0.median.reset(sim.dir, countries = NULL, joint.male = FALSE)
+e0.median.reset(sim.dir, countries = NULL, joint.male = FALSE, \dots)
 
-e0.shift.prediction.to.wpp(sim.dir, joint.male = FALSE, ...)
+e0.shift.prediction.to.wpp(sim.dir, joint.male = FALSE,
+    subdir = "predictions", \dots)
 }
 \arguments{
   \item{sim.dir}{Directory containing the prediction object.}
@@ -36,12 +38,13 @@ e0.shift.prediction.to.wpp(sim.dir, joint.male = FALSE, ...)
   \item{from}{Year from which the offset/reset should start. By default, it starts at the first prediction period.}
   \item{to}{Year until which the offset/reset should be done. By default, it is set to the last prediction period.}
   \item{factors}{It should be a vector where each element corresponds to one time period. The adjustment of male medians is done as \code{e0m(t) = e0f(t) - gap(t)*factor(t)}.}
-  \item{\dots}{Additional arguments passed to the underlying adjustment function. It can be \code{verbose} to show/hide the progress of the adjustment and \code{wpp.year} to adjust it to if it differs from the wpp year of the simulation.}
+  \item{subdir}{Subdirectory of \code{sim.dir} containing the predictions.}
+  \item{\dots}{Additional arguments passed to the underlying adjustment function. For \code{e0.shift.prediction.to.wpp} it can be \code{stat} with values \dQuote{median} (default) or \dQuote{mean} to specify which statistics should be adjusted; \code{verbose} to show/hide the progress of the adjustment and \code{wpp.year} to adjust it to if it differs from the wpp year of the simulation. For the other functions it can be \code{subdir} to specify the location of the prediction.}
 }
 \details{
 The function \code{e0.median.set} can be used to set the medians of the given country to specific values. Function \code{e0.median.shift} can be used to offset the medians by a specific constant, or to reset the medians to their original BHM values. Function \code{e0.median.adjust.jmale} adjusts male medians using factors that can expand or shrink the female-male gap. 
 
-Function\code{e0.shift.prediction.to.wpp} shifts the projected medians so that they correspond to the values found in the \code{e0Fproj} (\code{joint.male = FALSE}) or \code{e0Mproj} (\code{joint.male = TRUE}) datasets of the \pkg{wpp} package that either corresponds to the package used for the simulation itself or is given by the \code{wpp.year} argument. If using \pkg{wpp2022}, the dataset name is automatically adjusted depending if it is an annual or a 5-year simulation.
+Function\code{e0.shift.prediction.to.wpp} shifts the projected medians or means (if \code{stat} is \dQuote{mean}) so that they correspond to the values found in the \code{e0Fproj} (\code{joint.male = FALSE}) or \code{e0Mproj} (\code{joint.male = TRUE}) datasets of the \pkg{wpp} package that either corresponds to the package used for the simulation itself or is given by the \code{wpp.year} argument. If using \pkg{wpp2022} or higher, the dataset name is automatically adjusted depending if it is an annual or a 5-year simulation. Note that regardless if it is an adjustment of the median or mean, the corresponding offset is always converted to a shift of the median. 
 
 Function \code{e0.median.reset} resets medians of the given countries to the original values. By default it deletes adjustments for all countries.
 
diff --git a/man/e0.predict.Rd b/man/e0.predict.Rd
index b4aae13..62390a2 100644
--- a/man/e0.predict.Rd
+++ b/man/e0.predict.Rd
@@ -13,8 +13,8 @@ e0.predict(mcmc.set = NULL, end.year = 2100,
     sim.dir = file.path(getwd(), "bayesLife.output"), replace.output = FALSE, 
     predict.jmale = TRUE, nr.traj = NULL, thin = NULL, burnin = 10000, 
     use.diagnostics = FALSE, save.as.ascii = 0, start.year = NULL,  
-    output.dir = NULL, low.memory = TRUE, ignore.last.observed = FALSE,
-    seed = NULL, verbose = TRUE, \dots)
+    output.dir = NULL, subdir = "predictions", low.memory = TRUE, 
+    ignore.last.observed = FALSE, seed = NULL, verbose = TRUE, \dots)
 }
 
 \arguments{
@@ -31,6 +31,7 @@ e0.predict(mcmc.set = NULL, end.year = 2100,
   \item{save.as.ascii}{Either a number determining how many trajectories should be converted into an ASCII file, or \dQuote{all} in which case all trajectories are converted. It should be set to 0, if no conversion is desired (default).}
   \item{start.year}{This argument should be only used if the start year of the prediction is before or at the present year of the MCMC run (see Details below). By default the prediction starts in the next time period after the present year (passed to \code{\link{run.e0.mcmc}}).}
   \item{output.dir}{Directory into which the resulting prediction object and the trajectories are stored. If it is \code{NULL}, it is set to either \code{sim.dir}, or to \code{output.dir} of \code{mcmc.set$meta} if \code{mcmc.set} is given.}
+  \item{subdir}{Subdirectory of \code{output.dir} to store the predictions. It is defined relative to \code{output.dir} and can only have one level.}
   \item{low.memory}{Logical indicating if the prediction should run in a low-memory mode. If it is \code{FALSE}, the whole traces of all parameters, including the burnin, are loaded into memory. Otherwise, burnins are discarded and parameters are loaded as they are needed and are not kept in the memory.}
   \item{ignore.last.observed}{Logical. By default, the prediction (or imputation) for each country starts one time period after the last observed data point for that country defined by the \dQuote{last.observed} column in the data. If this argument is set to \code{TRUE}, the prediction ignores that \dQuote{last.observed} value and starts at the last data point found in the data. This allows to exclude some time periods from the estimation, but include them in the prediction.}
   \item{seed}{Seed of the random number generator. If \code{NULL} no seed is set. It can be used to generate reproducible projections.}
@@ -44,7 +45,7 @@ The projection is run for all missing values before the present year, if any. Me
 
 A special case is when the argument \code{start.year} is given that is smaller or equal the present year. In such a case, imputed missing values before present year are treated as ordinary predictions (trajectories are kept). All historical data between start year and present year are used as projections.
  
-The resulting prediction object is saved into \file{\{output.dir\}/predictions}. Trajectories for all countries are saved into the same directory in a binary format, one file per country. At the end of the projection, if \code{save.as.ascii} is larger than 0, the function converts the given number of trajectories into a CSV file of a UN-specific format. They are selected by equal spacing (see function \code{\link{convert.e0.trajectories}} for more details on the conversion). In addition, two summary files are created: one in a user-friendly format, the other using a UN-specific coding of the variants and time (see \code{\link{write.e0.projection.summary}} for more details).
+The resulting prediction object is saved into \file{\{output.dir\}/\{subdir\}}. Trajectories for all countries are saved into the same directory in a binary format, one file per country. At the end of the projection, if \code{save.as.ascii} is larger than 0, the function converts the given number of trajectories into a CSV file of a UN-specific format. They are selected by equal spacing (see function \code{\link{convert.e0.trajectories}} for more details on the conversion). In addition, two summary files are created: one in a user-friendly format, the other using a UN-specific coding of the variants and time (see \code{\link{write.e0.projection.summary}} for more details).
 
 }
 \value{
diff --git a/man/e0.predict.extra.Rd b/man/e0.predict.extra.Rd
index ed8de2a..b78b40c 100644
--- a/man/e0.predict.extra.Rd
+++ b/man/e0.predict.extra.Rd
@@ -9,12 +9,13 @@ Using the posterior parameter samples the function generates posterior trajector
 }
 \usage{
 e0.predict.extra(sim.dir = file.path(getwd(), 'bayesLife.output'), 
-    prediction.dir = sim.dir, countries = NULL, 
+    prediction.dir = sim.dir, subdir = "predictions", countries = NULL, 
     save.as.ascii = 1000, verbose = TRUE, \dots)
 	}
 \arguments{
   \item{sim.dir}{Directory with the MCMC simulation results.}
   \item{prediction.dir}{Directory where the prediction object and the trajectories are stored.}
+  \item{subdir}{Subdirectory of \code{prediction.dir} containing the predictions.}
   \item{countries}{Vector of country codes for which the prediction should be made. If it is \code{NULL}, the prediction is run for all countries that are included in the MCMC object but for which no prediction was generated.}
   \item{save.as.ascii}{Either a number determining how many trajectories should be converted into an ascii file, or \dQuote{all} in which case all trajectories are converted. It should be set to 0, if no converions is desired. Note that the convertion is done on all countries.}
   \item{verbose}{Logical switching log messages on and off.}
diff --git a/man/e0.predict.subnat.Rd b/man/e0.predict.subnat.Rd
index 8056cb2..05fff65 100644
--- a/man/e0.predict.subnat.Rd
+++ b/man/e0.predict.subnat.Rd
@@ -14,8 +14,8 @@ e0.predict.subnat(countries, my.e0.file,
     sim.dir = file.path(getwd(), "bayesLife.output"), 
     method = c("ar1", "shift", "scale"),
     predict.jmale = FALSE, my.e0M.file = NULL,
-    end.year = 2100, start.year = NULL, output.dir = NULL, 
-    annual = NULL, nr.traj = NULL, seed = NULL, 
+    end.year = 2100, start.year = NULL, subdir = "predictions", 
+    output.dir = NULL, annual = NULL, nr.traj = NULL, seed = NULL, 
     ar.pars = NULL, save.as.ascii = 0, verbose = TRUE, 
     jmale.estimates = NULL, \dots)
                 
@@ -35,6 +35,7 @@ subnat.gap.estimates(annual = FALSE)
   \item{my.e0M.file}{Tab-separated ASCII file containing the subnational male e0 data.}
  \item{end.year}{End year of the projections.}
  \item{start.year}{Start year of the projections. By default, projections start at the same time point as the national projections.} 
+ \item{subdir}{Subdirectory of \code{sim.dir} containing the national predictions.}
    \item{output.dir}{Directory into which the resulting prediction objects and the trajectories are stored. See below for details.}
    \item{annual}{Logical indicating if the subnational projection should be on an annual scale or a 5-year scale. By default, 
 the scale is matched to the national simulation given by \code{sim.dir}. If given, the scale must match to the scale of the subnational data provided in \code{my.e0.file}. 
diff --git a/man/get.e0.prediction.Rd b/man/get.e0.prediction.Rd
index a5bd7e5..b839b57 100644
--- a/man/get.e0.prediction.Rd
+++ b/man/get.e0.prediction.Rd
@@ -3,21 +3,25 @@
 \alias{get.e0.jmale.prediction}
 \alias{has.e0.prediction}
 \alias{has.e0.jmale.prediction}
+\alias{available.e0.predictions}
 
 \title{
 Accessing a Prediction Object
 }
 \description{
-Function \code{get.e0.prediction} retrieves results of a prediction and creates an object of class \code{\link{bayesLife.prediction}}. Function \code{has.e0.prediction} checks an existence of such results. Analogously, functions \code{get.e0.jmale.prediction} and \code{has.e0.jmale.prediction} retrieve and check an existence of male predictions from a given female prediction object.
+Function \code{get.e0.prediction} retrieves results of a prediction and creates an object of class \code{\link{bayesLife.prediction}}. Function \code{has.e0.prediction} checks an existence of such results. Analogously, functions \code{get.e0.jmale.prediction} and \code{has.e0.jmale.prediction} retrieve and check an existence of male predictions from a given female prediction object. Function \code{available.e0.predictions} lists predictions available in the given simulation directory.
 }
 \usage{
-get.e0.prediction(mcmc = NULL, sim.dir = NULL, joint.male = FALSE, mcmc.dir = NULL)
+get.e0.prediction(mcmc = NULL, sim.dir = NULL, joint.male = FALSE, 
+    mcmc.dir = NULL, subdir = "predictions")
 
-has.e0.prediction(mcmc = NULL, sim.dir = NULL)
+has.e0.prediction(mcmc = NULL, sim.dir = NULL, subdir = "predictions")
 
 get.e0.jmale.prediction(e0.pred)
 
 has.e0.jmale.prediction(e0.pred)
+
+available.e0.predictions(mcmc = NULL, sim.dir = NULL, full.names = FALSE)
 }
 
 \arguments{
@@ -25,9 +29,13 @@ has.e0.jmale.prediction(e0.pred)
   \item{sim.dir}{Directory where the prediction is stored.  It should correspond to the value of the \code{output.dir} argument used in the \code{\link{e0.predict}} function.  Only relevant if \code{mcmc} is \code{NULL}.}
   \item{joint.male}{Logical. If \code{TRUE}, the function is applied to a male prediction that was generated using the joint female-male model implemented in the function \code{\link{e0.jmale.predict}}.}
   \item{mcmc.dir}{Optional argument to be used only in a special case when the mcmc object contained in the prediction object was estimated in different directory than in the one to which it points to (for example due to moving or renaming the original directory). The argument causes that the mcmc is redirected to the given directory.}
+  \item{subdir}{Subdirectory of \code{sim.dir} for this particular prediction.}
+  \item{full.names}{Logical. If \code{TRUE}, the directory names are given as full paths, otherwise (default) only the base names.}
   \item{e0.pred}{Object of class \code{\link{bayesLife.prediction}}.}
 }
 \details{If \code{mcmc} is not \code{NULL}, the search directory is set to \code{mcmc$meta$output.dir}. This approach assumes that the prediction was stored in the same directory as the MCMC simulation, i.e. the \code{output.dir} argument of the \code{\link{e0.predict}} function was set to \code{NULL}. If it is not the case, the argument \code{mcmc.dir} should be used.
+
+Usually, all predictions are stored in the subdirectory \dQuote{predictions} of the simulation directory. If the subdirectory has a different name, the argument \code{subdir} should be used. This allows to keep multiple predictions in one (MCMC) simulation directory. The function \code{available.e0.predictions} can be used to view all available predictions in the simulation directory.
 	
 Function \code{get.e0.jmale.prediction} extracts male projections from the \code{e0.pred} objects (which should be a female prediction object), if the male prediction was generated using the  \code{\link{e0.jmale.predict}} function. \code{has.e0.jmale.prediction} checks if such male prediction was generated.
 }
@@ -36,7 +44,10 @@ Function \code{get.e0.jmale.prediction} extracts male projections from the \code
 	Functions \code{has.e0.prediction} and \code{has.e0.jmale.prediction} return a logical indicating if a prediction exists. 
 	
 	Functions \code{get.e0.prediction} and \code{get.e0.jmale.prediction} return an
-object of class \code{\link{bayesLife.prediction}}.}
+object of class \code{\link{bayesLife.prediction}}.
+
+    Function \code{available.e0.predictions} returns a vector of directory names containing e0 predictions. 
+}
 
 \author{
 Hana Sevcikova
diff --git a/man/include.Rd b/man/include.Rd
index 40cc078..3b010e1 100644
--- a/man/include.Rd
+++ b/man/include.Rd
@@ -6,6 +6,7 @@
 \alias{include_2010}
 \alias{include_2012}
 \alias{include_2022}
+\alias{include_2024}
 
 \docType{data}
 \title{
@@ -15,6 +16,7 @@ Inclusion Codes
 Datasets containing codes that determine which countries are to be included into a simulation or/and projections. 
 }
 \usage{
+data(include_2024)
 data(include_2022)
 data(include_2019)
 data(include_2017)
diff --git a/man/write.projection.summary.Rd b/man/write.projection.summary.Rd
index 02eac40..7d68612 100644
--- a/man/write.projection.summary.Rd
+++ b/man/write.projection.summary.Rd
@@ -8,10 +8,11 @@ The function creates two files containing projection summaries, such as the medi
 }
 \usage{
 write.e0.projection.summary(dir = file.path(getwd(), "bayesLife.output"), 
-    output.dir = NULL, revision = NULL, adjusted = FALSE)
+    subdir = "predictions", output.dir = NULL, revision = NULL, adjusted = FALSE)
 }
 \arguments{
   \item{dir}{Directory containing the prediction object. It should correspond to the \code{output.dir} argument of the \code{\link{e0.predict}} function.}
+  \item{subdir}{Subdirectory of \code{dir} containing the predictions.}
   \item{output.dir}{Directory in which the resulting file will be stored. If \code{NULL} the same directory is used as for the prediction.}
   \item{revision}{UN revision number. If \code{NULL} it is determined from the corresponding WPP year: WPP 2008 corresponds to revision 13, every subsequent WPP increases the revision number by one. Used as a constant in the second file only.}
    \item{adjusted}{Logical. By default the function writes summary using the original BHM projections. If the projection medians are adjusted (using e.g. \code{\link{e0.median.set}}), setting this argument to \code{TRUE} causes writing the adjusted projections.}
diff --git a/tests/run_tests.R b/tests/run_tests.R
index 0658f20..04677eb 100644
--- a/tests/run_tests.R
+++ b/tests/run_tests.R
@@ -26,10 +26,10 @@ if(!cran) {
 	test.subnational.predictions()
 	test.my.locations.extra()
 	test.reproduce.simulation()
-	for (wpp in c(2019, 2022)){
+	for (wpp in c(2019, 2024)){
 	    test.estimate.mcmc(wpp.year = wpp)
 	    test.run.annual.simulation(wpp.year = wpp)
 	}
-	for(wpp.year in c(2010, 2012, 2015, 2017, 2022))
+	for(wpp.year in c(2010, 2012, 2015, 2017, 2022, 2024))
 	    test.get.wpp.data(wpp.year)
 }