Array function coding error for multlevel logit model with multivariate priors

I need to put a period to this thread. Below I am posting a Stan code chunk for running a multilevel logit model with multivariate priors, hopefully without any issue. Note that I silence the lines containing array, since I haven’t figured it out yet and will be working on it. Will update once I have solid answers.
stanDiscoursePost20220402.R (3.4 KB)
Attached please also find an updated R script file. Thank you all for help!

data {
  int<lower=0> N;                     // num individuals
  int<lower=1> K;                     // num ind predictors
  int<lower=1> J;                     // num groups
  int<lower=1> L;                     // num group predictors
  // array[N] int<lower=1, upper=J> jj;  // group for individual
  int<lower=1, upper=J> jj[N];  // group for individual
  matrix[N, K] x;                     // individual predictors
  row_vector[L] u[J];           // group predictors
  // vector[N] y;                        // outcomes
  int<lower=0, upper=1> y[N];
  // vector<lower=0, upper=1>[N] y;

parameters {
  corr_matrix[K] Omega;        // prior correlation
  vector<lower=0>[K] tau;      // prior scale
  matrix[L, K] gamma;          // group coeffs
  // array[J] vector[K] beta;     // indiv coeffs by group
  vector[K] beta[J];
  real<lower=0> sigma;         // prediction error scale

model {
  tau ~ cauchy(0, 2.5);
  Omega ~ lkj_corr(2);
  to_vector(gamma) ~ normal(0, 5);
    // array[J] row_vector[K] u_gamma;
    row_vector[K] u_gamma[J];
    for (j in 1:J) {
      u_gamma[j] = u[j] * gamma;
    beta ~ multi_normal(u_gamma, quad_form_diag(Omega, tau));
  for (n in 1:N) {
    real eta;
    real prob;
    eta <- x[n] * beta[jj[n]]; 
    prob <- 1/(1+exp(-1*eta));  
    y[n] ~ bernoulli(prob);
    // y[n] ~ inv_logit(x[n] * beta[jj[n]]);
    // y[n] ~ bernoulli_logit(x[n] * beta[jj[n]]);