Rick and I discussed several ways to do what you want. The simplest way to do this is with a 3dcalc command. If you have a group of censor files, one for each run, and each of those censor files is a single column of numbers with 1 for values to keep and 0 for values to remove (censor), you can then use this command:

3dcalc \
-a epi_r01+orig -b epi_r02+orig -c epi_r03+orig \
-d censor01.1D -e censor02.1D -f censor03.1D \
-expr '(a*d+b*e+c*f)/(d+e+f)' \
-prefix epi_censor_mean

This works because 1D files by default are interpreted as values that vary with time and are applied in order at each sub-brick. The example here only shows three dataset runs, but it can be expanded up to thirteen runs with another thirteen censor files to fill up the alphabet of input parameters. If you need to do more than 13, then you can do this more complicated method, but still not too bad really:

#!/bin/tcsh
# censorlist.csh
# given multiple runs of epi (multiple 3D+time datasets)
# and a censor list in a column of 0 (censor) and 1 (include) for each run,
# compute mean dataset across runs
#
# this example uses datasets with names like epi_r01+orig, epi_r02+orig, ...
# and censor files censor01.1D, censor02.1D, ...
#
set base_dset = "epi_r"
set base_censor = "censor"
set nruns = 3
set runs = ( `count -digits 2 1 $nruns` )
# make sub-bricks 0 where they are censored
foreach run ( $runs )
3dcalc -a "${base_dset}${run}+orig" -b "${base_censor}${run}.1D" \
-expr 'a*b' -prefix "${base_dset}${run}_censor+orig"
end

# find the number of times each sub-brick is included (not censored)
3dMean -prefix "${base_censor}_sum" -sum ${base_censor}??.1D

# sum up sub-bricks across runs
3dMean -prefix "${base_dset}_sum" -sum ${base_dset}??_censor+orig.HEAD

# get mean by dividing by the number of times the sub-brick is included
3dcalc -a ${base_dset}_sum+orig -b ${base_censor}_sum.1D \
-expr 'a/b' -prefix "${base_dset}_censor_mean_long"