Building a multivariate nonlinear hierarchical model

I am trying to build up the following multivariate hierarchical model using an ODE. The data set consists of time series measurements for J subjects. Thus for each subject observed measurements (y_{j}) are modeled as:
y_j = f_{j}(\theta_{j}) + e_{j}, \quad j = 1:J, \theta_{j} = (\theta_{1}, \dots, \theta_{K}) ,
where f_{j}(\cdot ) is solution of ODE for parameter set \theta_{j} and e_{j} \sim MultiNormal(0, \sigma^{2} I{n_{j}}), with I as identity matrix and n_{j} is the number of measurements for subject j.

And log(\theta_{j}) = a +b_{j}, where b_{j} \sim MultiNormal(0, \Sigma), a \sim MultiNormal(mu_a, \Lambda) with \Lambda a diagonal matrix. (Note: I have logged parameters as some parameter values are large and some are very small).

At this stage, I am coding the priors and computing log(\theta_{j}) . Following Stans manual, I gave \Sigma an LKJ prior. The code runs okay, but I am stuck with constraining one parameter (parameter d \in [0,1]) in the parameter set, which must be between 0 and 1 on the original (non logged) scale.

Currently, taking exponential of l_theta give values for d larger than 1. How should I modify the code so that d is between 0 and 1 on the original scale?

 functions {
     real[] ode(real t,
     real[] y,
      real[] theta,
     real[] x_r, // x_r:real data, x_i: integer data
     int[] x_i) {
      real dydt[4];
     real A       =  theta[1];
    real beta    = theta[2];
      real delta   = theta[3];
      real gamma   = theta [4];
     real omega   = theta[5];
      real alpha   = theta[6];
      real kappa   = theta[7];
      real eta     = theta[8];
     real d       = theta[9];

        dydt[1] =  A - gamma*y[1]-beta * y[3] * y[1];                 
     dydt[2] = beta * y[1] * y[3] - delta * y[2];                  
    dydt[3] = (1-d)*omega * y[2]  - kappa * y[3] -alpha* y[4]*y[3] ;          
      dydt[4] = eta *y[3]* y[4];   
   return dydt;
  data {
      int<lower = 1> J;  //number of subjects
      int<lower =1> K; //number of parameters
      vector[K] mu_a;   //mean vector of a
  transformed data {
        real x_r[0];
       int x_i[0];

    matrix[J,K] a;
    vector<lower = 0>[K] tau; // sd of b
    cholesky_factor_corr [K] Omega;
    matrix[K,J] z;
 transformed parameters {
        matrix[J,K] b;
        matrix[J,K] l_theta;    //parameters on the log scale

      b = (diag_pre_multiply(tau, Omega)*z)';
      l_theta = a + b;

model {
      tau ~ cauchy(0,2.5);
      Omega ~ lkj_corr_cholesky(2); 
      to_vector(z) ~ normal(0,1);//this implies for(k in 1:K) z[k]~normal(0,1)
      to_vector(a) ~ normal(mu_a, 5);

The R code used to run the stan program is :

stan_data <- list(  J = 1, #numb. of subjects
                           K = 9, #numb. of parameters 
                         mu_a = c(log(1.4e7), log(1.72e-10), log(0.14),log(0.14), log(1e4), log(0.001),log(3.48),log(1.29e-5), log(0.5))
       #mean vector of (A, beta, delta, gamma,omega,alpha, kappa,eta, d) on log scale
fake_data <- stan("model_fit.stan",         
                  data = stan_data,
                  chains = 1, iter = 1000,
a <- extract(fake_data, pars =c("a", 'b', 'l_theta'))

Also, is the following code used to sample a \sim MultiNormal(mu_a , \Lambda), with \Lambda = diag(rep(10, K)) correct? I read from Stan manual to use a normal distribution if the co-variance matrix of multinormal is diagonal.

to_vector(a) ~ normal(mu_a, 5);

I feel that there is something incorrect with the code as the density of l_theta for each parameter is concentrated at 0 and taking exponential do not give a proper distribution for theta parameters.