model {
  for(i in 2:nday){
    dayl = 0
    pi = 2*arcsin(1)
    Tr = .5*exp(p10*.5*(mint[i] + maxt[i]))
    lai = max(.1, Cf/LMA)
    trange = (maxt[i] - mint[i])/2
    gs = pow(2, 0.789798) / (0.37836)
    pp = lai*Nit / gs *p11*exp(0.011136*maxt[i])
    qq = -204.6453
    Caqqpp.sq = pow((Ca[i] + qq - pp), 2)
    ci = .5*(Ca[i] + qq - pp + pow( (Caqqpp.sq) - 4*(Ca[i]*qq - pp*4.22273), 1/2 ))

    e0 = 7.19298*pow(lai, 2) / (pow(lai, 2) + 2.1001)

    dec = -23.4*cos(360*(yearday+10)/365 *pi/180)*pi / 180
    mult = tan(dec)*tan(lat)

    dayl = ifelse(abs(mult) < 1, 24*arccos(-mult)/pi, ifelse(mult >= 1, 24, 0) )

    cps = e0*rad[i]*gs*(Ca[i] - ci) / (e0*rad[i] + gs*(Ca[i] - ci))
    G[i] = cps*(0.0156935*dayl + 0.0453194)

    Cf.pred[i] = (1-p5)*Cf[i-1] + G[i]*(1 - p2)*p3
    Cw.pred[i] = (1 - p6)*Cw[i-1] + G[i]*(1-p2)*(1-p3)*(1-p4)
    Cr.pred[i] = (1 - p7)*Cr[i-1] + G*(1-p2)*(1-p3)*p4
    Clit.pred[i] = Cf[i-1]*p5 + Cr[i-1]*p7 + (1 - Tr*(p8+p1))*Clit[i-1]
    Csom.pred[i] = p6*Cw[i-1] + Tr*p1*Clit[i-1] + (1 - Tr*p9)*Csom[i-1]

    Cf[i] ~ dnorm(Cf.pred[i], sd_cf)
    Cw[i] ~ dnorm(Cw.pred[i], sd_cw)
    Cr[i] ~ dnorm(Cr.pred[i], sd_cr)
    Clit[i] ~ dnorm(Clit.pred[i], sd_clit)
    Csom[i] ~ dnorm(Csom.pred[i], sd_csom)
  }

  for(i in 1:nday){
    Cf.obs[i] ~ dnorm(Cf[i], sd_cf_obs)
    Cw.obs[i] ~ dnorm(Cw[i], sd_cw_obs)
    Cr.obs[i] ~ dnorm(Cr[i], sd_cr_obs)
    Clit.obs[i] ~ dnorm(Clit[i], sd_clit_obs)
    Csom.obs[i] ~ dnorm(Csom[i], sd_csom_obs)
  }

  ## Initial conditions
  Cf[1] ~ dnorm(100, 2)
  Cw[1] ~ dnorm(9200, 25)
  Cr[1] ~ dnorm(100, 2)
  Clit[1] ~ dnorm(20, 2)
  Csom[1] ~ dnorm(11000, 25)

  ## Priors on parameters
  p1 ~ dunif((10^(-6)), .01)
  p2 ~ dunif(.2, .7)
  p3 ~ dunif(.01, .5)
  p4 ~ dunif(.01, .5)
  p5 ~ dunif((10^(-4)), .1)
  p6 ~ dunif((10^(-6)), .01)
  p7 ~ dunif((10^(-4)), .01)
  p8 ~ dunif((10^(-5)), .1)
  p9 ~ dunif((10^(-6)), .01)
  p10 ~ dunif(.05, .2)
  p11 ~ dunif(5, 20)

  ## Priors on process errors
  sd_cf ~ dnorm(0, 100)T(0,)
  sd_cw ~ dnorm(0, 100)T(0,)
  sd_cr ~ dnorm(0, 100)T(0,)
  sd_clit ~ dnorm(0, 100)T(0,)
  sd_csom ~ dnorm(0, 100)T(0,)

  sd_cf_obs ~ dnorm(0, 100)T(0,)
  sd_cw_obs ~ dnorm(0, 100)T(0,)
  sd_cr_obs ~ dnorm(0, 100)T(0,)
  sd_clit_obs ~ dnorm(0, 100)T(0,)
  sd_csom_obs ~ dnorm(0, 100)T(0,)
}