Question about Gumbel Copula estimation in Rstan

[Please include Stan program and accompanying data if possible]
Hello, everyone,
I am a beginner for Rstan, right now I encounter a problem regarding the Copula Bayesian estimation using Stan.
I try to estimate the joint probability of rainfall and streamflow and the parameter of Gumbel Copula. I assume the rainfall data comes from a Gamma distribution and the streamflow data comes from a GEV distribution. I found the Gumbel copula and GEV distribution are not available by default in Rstan, so I define the functions first.

      real gev_pdf_log(real xx, real mu, real sigma, real xi){ 
         real ystar;
         real xy_p1;
         xy_p1 = xi* ystar +1;
        return -log(sigma)-(1+1/xi)*log(xy_p1)-(xy_p1)^(-1/xi);
      real gev_cdf(real xx, real mu, real sigma, real xi){
         real ystar;
         real xy_p1;
         xy_p1 = xi* ystar +1;
        return exp(-(xy_p1)^(-1/xi));
      real gumbel_copula(real u, real v, real theta) {
         real neg_log_u = -log(u);
         real neg_log_v = -log(v);
         real theta_m1 = theta - 1;
         real temp = neg_log_u ^ theta + neg_log_v ^ theta;
        return theta_m1 * log(neg_log_u) + theta_m1 * log(neg_log_v) 
         + neg_log_u + neg_log_v - pow(temp, inv(theta)) 
         + (-2+2/theta)*log(temp)
         + log1p(theta_m1*pow(temp,-1/theta));
    data {
       int<lower=0> N;  // number of years
       vector[N] x;         // rainfall data series
       vector[N] y;         // streamflow data series
    parameters {
       real mu;
       real <lower=0>sigma;
       real xi;
       real<lower=0> alpha;
       real<lower=0> beta;
       real<lower=0,upper=1> tau; // kendall’s tau
    transformed parameters {
       real<lower=1> theta = 1 / (1 - tau);
    model {
      vector[N] summands; 
       // priors
       mu ~ normal(1000,100);
       sigma ~ normal(750,100);
       xi ~ uniform(-1,1);
       alpha ~ normal(5,5);
       beta ~ uniform(0,0.1);
       //tau ~ gamma(0.01,0.01);
       //tau ~ cauchy(0, 2.5);
        // model
       x ~ gamma(alpha,beta);
       for (n in 1:N)
         y[n] ~ gev_pdf(mu,sigma,xi);
       for (n in 1:N)
         summands[n] = gumbel_copula(gamma_cdf(x[n], alpha, beta),
                                     gev_cdf(y[n], mu, sigma,xi), theta);
       target += summands;

then I run it:

 dat <- read.csv("C:/Users/Administrator/Desktop/Rstan/downloaded/data.csv")
 dat$N <- NULL
 post <- stan("C:/Users/Administrator/Desktop/Rstan/downloaded/gumbel.stan", data = list(N = nrow(dat), x= dat$x, y = dat$y),iter=12000,warmup=2000,thin=10,chains=3)

I got the following error:

 Rejecting initial value:
   Log probability evaluates to log(0), i.e. negative infinity.
   Stan can't start sampling from this initial value.
   Initialization between (-2, 2) failed after 100 attempts.
   Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

The priors I set for parameters first are non-informative distributions, but not work. Then I tried to use the MLE method to estimate the parameters of each distribution and set the prior distributions based on the MLE results which are shown in the Code. I don’t know the reason and have been confused by days`.
the Code and data are attached.Could you help me have a look? Thanks a lot!
data.csv (583 Bytes)

You may just need to specify the init_r argument to some number closer to 0 than to its default value of 2.