Here is an R solution. Other questions on this site demonstrate how to read and write raster files in R, so let's get right to the solution. First, some random data to illustrate the procedure:
n.rows <- 1000
n.cols <- 1000
r <- outer(1:n.rows, 1:n.cols,
function(x,y) sin(x/100) * cos(((100+x)/(1+(y/100)^2)))^2) + rnorm(n.rows*n.cols, 0, .10)
Here's the code. It applies a user-supplied function f to each 6 by 6 block of the raster r; f finds and averages the five largest values. The code uses sapply to break the array into columns of width 6 and calls aggregate to spit out a column of statistics within each such vertical band. (I coded it in this order--split the columns first--because R stores arrays in column-major order. This will maximize the locality of reference, which is important for speed on very large arrays.)
f <- function(x) {y <- as.vector(x); mean(y[order(y, decreasing=TRUE)[1:5]])}
aggregate <- function(x, k=dim(x)[2], fun) sapply(seq(1, dim(x)[1]+1-k, k),
function(i) fun(x[1:k+i-1, ]))
agg <- sapply(seq(1, dim(r)[2]+1-6, 6), function(j) aggregate(r[, 1:6+j-1], k=6, fun=f))
This takes about two seconds for the sample grid of a million cells. Plots help us compare them:
library(raster)
r.r <- raster(r); agg.r <- raster(agg)
par(mfrow=c(1,2))
image(r.r, zlim=c(min(r), max(r)), col=terrain.colors(300), main="Original")
image(agg.r, zlim=c(min(r), max(r)), col=terrain.colors(300), main="Aggregate")

It is evident that the aggregate tends to pick out the larger values in each block.