Non converging model, which should be reparameterized

I am modelling 4 variables A_i~N, which in sum give Z=sum(A_i) and Z follows a generalized normal distribution.
The data shows that the standard deviations D[A_i(Z)] are a function of Z and rho between A rho_ij(Z) as well. The code compiles and fits parameters however, stan throws a warning that the model should be reparametrized.

“WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model”

I am grateful for any advice and thank you in advance!

// MODEL 104
functions{
real func_linear(real z, real zmin, real zmax, real rho_zmin, real rho_zmax ){
// linear interpolation between zmin and zmax
real y;
y = rho_zmin + (z- zmin) * (rho_zmax - rho_zmin) / (zmax - zmin) ;
return y;
}

    real gen_normal_lpdf( vector x, real mu_, real alpha_, real beta_ ) {
        vector[rows(x)] lp;
        int N;

        N = rows(x);
        for(i in 1:N){
           // pdf according to wikipedia
           // https://en.wikipedia.org/wiki/Generalized_normal_distribution
           lp[ i ] = log(beta_) - (log(2*alpha_)+lgamma(1/beta_)) - pow( fabs(x[i]-mu_)/alpha_, beta_ );
        }
        // returns loglikelihood to maximise it (not neg. likelihood!)
        return sum( lp );
    }

}

data {
int < lower = 1 > data_n; // number of values
int < lower = 1 > axes_n; // number of axes

matrix[data_n, axes_n] A;   // data values of axe loads
vector[data_n] Z;

real < lower = 0, upper = 35 > prior_A_mean;          // mean of axe load for prior distribution
real < lower = 0, upper = 10 > prior_A_std;           // std dev of axe load for prior distribution
real < lower = 0, upper = 1 > prior_A_cov;           // CoV for the priors

real  prior_rho_c12;           //
real  prior_rho_c13;           //
real  prior_rho_c14;           //
real  prior_rho_c23;           //
real  prior_rho_c24;           //
real  prior_rho_c34;           //

real  prior_rho_d12;           //
real  prior_rho_d13;           //
real  prior_rho_d14;           //
real  prior_rho_d23;           //
real  prior_rho_d24;           //
real  prior_rho_d34;           //

real < lower = 0 > Z_0;      // beta parameter for the Z variable
real < lower = 0 > Z_1;      // beta parameter for the Z variable
real < lower = 0 > Z_beta;      // beta parameter for the Z variable

}

transformed data {
vector <lower=0, upper = 35> [data_n] A1;
vector <lower=0, upper = 35> [data_n] A2;
vector <lower=0, upper = 35> [data_n] A3;
vector <lower=0, upper = 35> [data_n] A4;

A1 = A[:, 1];
A2 = A[:, 2];
A3 = A[:, 3];
A4 = A[:, 4];

}

parameters {
real < lower = 0, upper = 30 > A1_mean; // mean of axe load
real < lower = 0, upper = 30 > A2_mean; // mean of axe load
real < lower = 0, upper = 30 > A3_mean; // mean of axe load
real < lower = 0, upper = 30 > A4_mean; // mean of axe load

real < lower = -1, upper = 1 > rho_c12;               //
real < lower = -1, upper = 1 > rho_c13;               //
real < lower = -1, upper = 1 > rho_c14;               //
real < lower = -1, upper = 1 > rho_c23;               //
real < lower = -1, upper = 1 > rho_c24;               //
real < lower = -1, upper = 1 > rho_c34;               //

real < lower = -1, upper = 1 > rho_d12;               //
real < lower = -1, upper = 1 > rho_d13;               //
real < lower = -1, upper = 1 > rho_d14;               //
real < lower = -1, upper = 1 > rho_d23;               //
real < lower = -1, upper = 1 > rho_d24;               //
real < lower = -1, upper = 1 > rho_d34;               //

real < lower = 0 > A1_std_c;               //
real < lower = 0 > A2_std_c;               //

real < lower = 0 > A1_std_d;               //
real < lower = 0 > A2_std_d;               //

real < lower = 0, upper = axes_n * 20 >     Z_alpha;  //  parameter for Z

}

transformed parameters{
matrix[axes_n, axes_n] rho [data_n] ; // correlation coefficients
matrix < lower = -1000, upper = 1000 >[axes_n,axes_n] COVMAT [data_n] ; // covariance matrix
cholesky_factor_cov [axes_n] L_COVMAT[data_n] ; // choleski of COVMAT

vector <lower=0, upper = 35> [axes_n]  A_mean;          //  vector of mean values

real < lower = 0, upper = axes_n * 35  >  Z_mean;       //  mean of total axe load
real < lower = 0, upper = axes_n * 200 >  Z_var;        //  variance of total axe load

vector < lower = 0, upper = 20 > [data_n] A1_std;   //  std dev of axe load
vector < lower = 0, upper = 20 > [data_n] A2_std;   //  std dev of axe load
vector < lower = 0, upper = 20 > [data_n] A3_std;   //  std dev of axe load
vector < lower = 0, upper = 20 > [data_n] A4_std;   //  std dev of axe load

real   p;       //  parameter of the quadratic equation
real   q;       //  parameter of the quadratic equation

A_mean[1]= A1_mean;     //  building the vector with mean values
A_mean[2]= A2_mean;
A_mean[3]= A3_mean;
A_mean[4]= A4_mean;

//  calculate the mean of the total load Z
Z_mean = A1_mean +  A2_mean + A3_mean + A4_mean ;

//  caculate Z_var as a function of Z_alpha
Z_var = square( Z_alpha ) * tgamma( 3 / Z_beta ) / tgamma( 1 / Z_beta ) ;

for (i in 1:data_n){
    //  caculate the rhos as a function of the total axe load Z data
    rho[i,1,1] = 1;
    rho[i,1,2] = func_linear( Z[i], Z_0, Z_1, rho_c12,  rho_d12 ) ;
    rho[i,1,3] = func_linear( Z[i], Z_0, Z_1, rho_c13,  rho_d13 ) ;
    rho[i,1,4] = func_linear( Z[i], Z_0, Z_1, rho_c14,  rho_d14 ) ;
    rho[i,2,1] = rho[i,1,2];
    rho[i,2,2] = 1;
    rho[i,2,3] = func_linear( Z[i], Z_0, Z_1, rho_c23,  rho_d23 ) ;
    rho[i,2,4] = func_linear( Z[i], Z_0, Z_1, rho_c24,  rho_d24 ) ;
    rho[i,3,1] = rho[i,1,3];
    rho[i,3,2] = rho[i,2,3];
    rho[i,3,3] = 1;
    rho[i,3,4] = func_linear( Z[i], Z_0, Z_1, rho_c34,  rho_d34 ) ;
    rho[i,4,1] = rho[i,1,4];
    rho[i,4,2] = rho[i,2,4];
    rho[i,4,3] = rho[i,3,4];
    rho[i,4,4] = 1;

    //  calculate STD as a function of Z
    A2_std[i] = func_linear( Z[i], Z_0, Z_1, A2_std_c, A2_std_d) ;
    A3_std[i] = func_linear( Z[i], Z_0, Z_1, A2_std_c, A2_std_d) ;  // use A2 values
    A4_std[i] = func_linear( Z[i], Z_0, Z_1, A1_std_c, A1_std_d) ;  // use A1 values

    //  we calculate A1_std as a function of Z_var and the covariance matrix
    p = 2*( A2_std[i] *rho[i,1,2] + A3_std[i]*rho[i,1,3] + A4_std[i]*rho[i,1,4] );
    q = 0;
    q += -Z_var;
    q += A2_std[i]*A2_std[i]*rho[i,2,2] + A2_std[i]*A3_std[i]*rho[i,2,3] + A2_std[i]*A4_std[i]*rho[i,2,4];
    q += A3_std[i]*A2_std[i]*rho[i,3,2] + A3_std[i]*A3_std[i]*rho[i,3,3] + A3_std[i]*A4_std[i]*rho[i,3,4];
    q += A4_std[i]*A2_std[i]*rho[i,4,2] + A4_std[i]*A3_std[i]*rho[i,4,3] + A4_std[i]*A4_std[i]*rho[i,4,4];

    A1_std[i] = -p/2 + sqrt( square( p/2) - q );

    COVMAT[i,1,1] = rho[i,1,1] * A1_std[i] * A1_std[i] + 1e-10;
    COVMAT[i,1,2] = rho[i,1,2] * A1_std[i] * A2_std[i];
    COVMAT[i,1,3] = rho[i,1,3] * A1_std[i] * A3_std[i];
    COVMAT[i,1,4] = rho[i,1,4] * A1_std[i] * A4_std[i];
    COVMAT[i,2,1] = rho[i,2,1] * A2_std[i] * A1_std[i];
    COVMAT[i,2,2] = rho[i,2,2] * A2_std[i] * A2_std[i] + 1e-10;
    COVMAT[i,2,3] = rho[i,2,3] * A2_std[i] * A3_std[i];
    COVMAT[i,2,4] = rho[i,2,4] * A2_std[i] * A4_std[i];
    COVMAT[i,3,1] = rho[i,3,1] * A3_std[i] * A1_std[i];
    COVMAT[i,3,2] = rho[i,3,2] * A3_std[i] * A2_std[i];
    COVMAT[i,3,3] = rho[i,3,3] * A3_std[i] * A3_std[i] + 1e-10;
    COVMAT[i,3,4] = rho[i,3,4] * A3_std[i] * A4_std[i];
    COVMAT[i,4,1] = rho[i,4,1] * A4_std[i] * A1_std[i];
    COVMAT[i,4,2] = rho[i,4,2] * A4_std[i] * A2_std[i];
    COVMAT[i,4,3] = rho[i,4,3] * A4_std[i] * A3_std[i];
    COVMAT[i,4,4] = rho[i,4,4] * A4_std[i] * A4_std[i] + 1e-10;

    L_COVMAT[i,:,:] = cholesky_decompose( COVMAT[i,:,:] );

}   // end loop over data

}

model {

//  rho parameter priors
target += normal_lpdf( rho_c12 | prior_rho_c12, 1 );
target += normal_lpdf( rho_c13 | prior_rho_c13, 1 );
target += normal_lpdf( rho_c14 | prior_rho_c14, 1 );
target += normal_lpdf( rho_c23 | prior_rho_c23, 1 );
target += normal_lpdf( rho_c24 | prior_rho_c24, 1 );
target += normal_lpdf( rho_c34 | prior_rho_c34, 1 );

target += normal_lpdf( rho_d12 | prior_rho_d12, 1 );
target += normal_lpdf( rho_d13 | prior_rho_d13, 1 );
target += normal_lpdf( rho_d14 | prior_rho_d14, 1 );
target += normal_lpdf( rho_d23 | prior_rho_d23, 1 );
target += normal_lpdf( rho_d24 | prior_rho_d24, 1 );
target += normal_lpdf( rho_d34 | prior_rho_d34, 1 );

target += normal_lpdf( A1_std_c | prior_A_std, prior_A_std * prior_A_cov );
target += normal_lpdf( A2_std_c | prior_A_std, prior_A_std * prior_A_cov );
target += normal_lpdf( A1_std_d | prior_A_std, prior_A_std * prior_A_cov );
target += normal_lpdf( A2_std_d | prior_A_std, prior_A_std * prior_A_cov );

//  load std priors
for (i in 1: data_n) {
    target += normal_lpdf( A1_std[i] | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A2_std[i] | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A3_std[i] | prior_A_std, prior_A_std * prior_A_cov );
    target += normal_lpdf( A4_std[i] | prior_A_std, prior_A_std * prior_A_cov );
}

//  load mean priors
target += normal_lpdf( A1_mean | prior_A_mean, prior_A_mean * prior_A_cov );
target += normal_lpdf( A2_mean | prior_A_mean, prior_A_mean * prior_A_cov );
target += normal_lpdf( A3_mean | prior_A_mean, prior_A_mean * prior_A_cov );
target += normal_lpdf( A4_mean | prior_A_mean, prior_A_mean * prior_A_cov );

//  Total load priors
target += normal_lpdf( Z_mean | prior_A_mean * axes_n, prior_A_cov * prior_A_mean * axes_n );
target += normal_lpdf( Z_var  | axes_n * square( prior_A_std ), prior_A_cov * axes_n * square( prior_A_std ) );
target += normal_lpdf( Z_alpha | sqrt( axes_n ) * prior_A_std * tgamma(1/Z_beta) / tgamma(3/Z_beta),  prior_A_cov*sqrt( axes_n ) * prior_A_std * tgamma(1/Z_beta) / tgamma(3/Z_beta) );

//  Likelihood Total axe load
target += gen_normal_lpdf( Z | Z_mean, Z_alpha, Z_beta );

for (i in 1: data_n) {
    A[i,:] ~ multi_normal_cholesky_lpdf(A_mean, L_COVMAT[i,:,:] );
}

}