model {
for (i in 1:N) {
Y[i] ~ dnorm(mu, tau)
Y.rep[i] ~ dnorm(mu, tau)
T.data.obs[i] <- pow((Y[i] - mean(Y[]))/sd(Y[]),3)
T.data.rep[i] <- pow((Y.rep[i] - mean(Y.rep[]))/sd(Y.rep[]),3)
T.para.obs[i] <- pow((Y[i] - mu)/sigma,3)
T.para.rep[i] <- pow((Y.rep[i] - mu)/sigma,3)
}
mu ~ dunif(-100, 100)
tau ~ dgamma(0.001, 0.001)
sigma <- 1/sqrt(tau)
T.data.obs.tot <- sum(T.data.obs[])
T.data.rep.tot <- sum(T.data.rep[])
T.para.obs.tot <- sum(T.para.obs[])
T.para.rep.tot <- sum(T.para.rep[])
P.data <- step(T.data.obs.tot - T.data.rep.tot)
P.para <- step(T.para.obs.tot - T.para.rep.tot)
}

Inits:
list(mu = 0, tau = 1)

Data:
list(N=66, Y=c(
28, 26, 33, 24, 34, -44, 27, 16, 40, -2,
29, 22, 24, 21, 25, 30, 23, 29, 31, 19,
24, 20, 36, 32, 36, 28, 25, 21, 28, 29,
37, 25, 28, 26, 30, 32, 36, 26, 30, 22,
36, 23, 27, 27, 28, 27, 31, 27, 26, 33,
26, 32, 32, 24, 39, 28, 24, 25, 32, 25,
29, 27, 28, 29, 16, 23))

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   P.data   0.0   0.0   1.0E-12   0.0   0.0   0.0   1   10000
   P.para   0.0   0.0   1.0E-12   0.0   0.0   0.0   1   10000
   T.data.obs.tot   -289.8   0.0   1.0E-12   -289.8   -289.8   -289.8   1   10000
   T.data.rep.tot   0.02773   18.65   0.1727   -37.1   -0.01383   37.07   1   10000
   T.para.obs.tot   -293.3   81.25   0.8554   -475.1   -285.5   -157.4   1   10000
   T.para.rep.tot   -0.444   31.57   0.2758   -62.93   -0.4164   62.88   1   10000
   mu   26.23   1.34   0.01353   23.6   26.23   28.88   1   10000
   sigma   10.87   0.9811   0.01047   9.169   10.79   13.0   1   10000
   tau   0.00866   0.001532   1.653E-5   0.005921   0.008591   0.0119   1   10000