Fitting GLMMs in JAGS: Non-conforming parameters in function

Question

I'm attempting to fit a Poisson GLMM in JAGS as per (AZuur, Hill & Reno 2013). I'm encountering the following error in R:

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  RUNTIME ERROR:
Non-conforming parameters in function :

This doesn't allude to the troublesome function in question, although I have an idea it may be an indexing problem with Nre - I am running models via a loop (for 217 species) but each species has a different combination of sites and years in which it has been observed. The goal is to fit the following code for three fixed effects (rain, ndvi, fire) and two random effects (site, year).

The code is question, including some dummy code:

count<- sample(0:20, 100, replace = TRUE)
rain<- rnorm(100, 0, 0.01)
ndvi<- rnorm(100, 0.2, 0.01)
burn<- rnorm(100, 0, 0.01)
site<- c(rep(1, 30), rep(2, 20), rep(3,50))
year<- c(rep(seq(2000:2009), 10))
temp<- as.data.frame(cbind(count, rain, ndvi, burn, site, year))

# Create dataframe and define bins and control variables
X<- model.matrix(~ rain + ndvi + burn, data = temp)
K <- ncol(X)
re<- as.numeric(temp$site)
Nre<- unique(temp$site)
win.data1<- list(
  Y = temp$count, # response
  X = X, # Covaraites
  K = K, # Num. betas
  N = nrow(temp), # Sample size
  Re = re, # Plot identifier
  Nre = Nre) # Number of plots

# Define JAGS model and save as .txt file 
sink("lmm.text")
cat("
        model {
        #1A. Diffuse normal priors for regression paramters
        for (i in 1:K) {beta[i] ~ dnorm(0, 0.0001) }
        
        #1B. Priors for random effect plot 
        for (i in 1:Nre){
        a[i] ~ dnorm(0, tau.plot)
        b[i] ~ dnorm(0, tau.year)
        }
        
        #1C. Priors for the two sigmas
        tau.plot<- 1/(sigma.plot*sigma.plot)
        tau.year<- 1/(sigma.year*sigma.year)
        tau.eps<- 1/(sigma.eps*sigma.eps)
        sigma.plot ~ dunif(0.001, 10)
        sigma.year ~ dunif(0.001, 10)
        sigma.eps ~ dunif(0.001, 10)
        
        #2. Likelihood of the data
        for (i in 1:N){
        Y[i] ~ dnorm(mu[i], tau.eps)
        mu[i]<- eta[i]
        eta[i]<- inprod(beta[], X[i,]) + a[Re[i]] + b[Re[i]] * X[i,2]
          }
        }
        ", fill = TRUE)
sink()

# Describe function to create initial values for JAGS
inits1<- function(){
  list(beta = rnorm(K, 0, 0.01),
       a = rnorm(Nre, 0, 1),
       b = rnorm(Nre, 0, 1), 
       sigma.eps = runif(1, 0.001, 10),
       sigma.plot = runif(1, 0.001, 10),
       sigma.year = runif(1, 0.001, 10))
}

# Define parameter list to track. Here each beta is handled later by the custom function uNames    
params1<- c("beta", "a", "b", "sigma.plot", "sigma.year", "sigma.eps")

# Run JAGS
J0<- jags(data = win.data1,
          inits = inits1,
          parameters = params1,
          model.file = "src/lmm.text",
          n.thin = 10,
          n.chains = 3,
          n.burnin = 4000,
          n.iter = 5000)

# Update chains to improve mixing
J1<- update(J0, n.iter = 10000, n.thin = 10)

# Save output
out<- J1$BUGSoutput

Once I have this up and running, the goal is to swap to a negative binomial, although I'd like to understand the error of my ways before I complicate things further.