Occupancy model - jags to julia - NUTS sampler

I am trying to translate an occupancy model from jags to julia and arrived at the following code:

@model function model1(y, n_cells, n_years, n_traps, n_period, z, z_su, ::Type{T} = Float64) where {T}
  # -- define / preallocate
  psi = Array{T}(undef, (n_years, n_cells))
  mu_su = Array{T}(undef, (n_years, n_cells, n_traps))
  theta = Array{T}(undef, (n_years, n_cells, n_traps))
  mu_y = Array{T}(undef, (n_cells, n_years, n_traps, n_period))
  p = Array{T}(undef, (n_cells, n_years, n_traps, n_period))
  # -- priors
  psi0 ~ filldist(Beta(1, 1), n_years)
  theta0 ~ filldist(Beta(1, 1), n_years)
  p0 ~ filldist(Beta(1, 1), n_years)
  int_psi = logit.(psi0)
  int_theta = logit.(theta0)
  int_p = logit.(p0)
  # -- Likelihood (for model structure)
  for yy in 1:n_years
    for i in 1:n_cells
      # Occupancy in grid cell i for year yy
      psi[yy, i] = logistic(int_psi[yy])
      z[yy, i] ~ Bernoulli(psi[yy, i])
      for j in 1:n_traps
        # Site-use at camera location j for year yy
        theta[yy, i, j] = logistic(int_theta[yy])
        mu_su[yy, i, j] = z[yy, i] * theta[yy, i, j]
        z_su[yy, i, j] ~ Bernoulli(mu_su[yy, i, j])
        for k in 1:n_period
          # detection probability at camera location j in period k and year yy
          p[i, yy, j, k] = logistic(int_p[yy])
          mu_y[i, yy, j, k] = z_su[yy, i, j] * p[i, yy, j, k]
          y[i, yy, j, k] ~ Bernoulli(mu_y[i, yy, j, k])
        end # End period
      end # End camera trap
    end # End grid cell
  end # End year
 end # End model

The MH sampler seems to work with dummy data (although I don’t get any standard errors), but I can’t get the NUTS sampler to work… It spits a lot of warnings about rejected proposals, but the
progress bar doesn’t change at all. So I am obviously doing some noob mistakes here…
I would really appreciate if a more experienced turing user would take a (quick) look a this code and point me to what is blatantly wrong in there… Many thanks!