volta.cokri2 <- function(mat.cokri, num.simu, retorna.tudo=FALSE, int.conf=0.95){
    ## Author: Wagner Bonat / Ana Beatriz Martins / Paulo Justiniano
    ## Modified for parallel and better use of RAM by: Benilton Carvalho
    require(parallel)
    nlinhas <- dim(mat.cokri[[1]])[1]
    f <- function(i){
        g <- mvrnorm(n=1, mat.cokri[[1]], mat.cokri[[2]])
        seq1 <- seq(1, nlinhas, by=2)
        gerado <- data.frame(g[seq1], g[seq1+1L])
        agl(gerado)
    }
    compos2 <- do.call(cbind, mclapply(1:num.simu, f))
    dim.vetor <- num.simu * 3
    sy1 <- seq(1, dim.vetor, by=3)
    amostra1 <- as.matrix(compos2[, sy1], ncol=num.simu)
    amostra2 <- as.matrix(compos2[, sy1+1L], ncol=num.simu)
    amostra3 <- as.matrix(compos2[, sy1+2L], ncol=num.simu)
    if (retorna.tudo) {
        retorna <- list(amostra1, amostra2, amostra3)
        return(retorna)
    }
    if (!retorna.tudo) {
        med1 <- rowMeans(amostra1)
        med2 <- rowMeans(amostra2)
        med3 <- rowMeans(amostra3)
        probs <- c(1-int.conf, int.conf)
        q1 <- t(apply(amostra1, 1, quantile, prob=probs))
        q2 <- t(apply(amostra2, 1, quantile, prob=probs))
        q3 <- t(apply(amostra3, 1, quantile, prob=probs))
        quantis <- cbind(q1, q2, q3)
        quantis <- data.frame(quantis)
        names(quantis) <- c("LI Areia", "LS Areia", "LI Silte", "LS Silte", "LI Argila", "LS Argila")
        resultado <- list(preditos=data.frame(Areia=med1, Site=med2, Argila=med3), intervalo=quantis)
        return(resultado)
    }
}