model { 
  # Observation model
  for (s in 1:nsurv) {
	M[s] <- n[yrsurv[s]]  # *exp(delta*season[s]) 	
    varO[s] <- max((M[s]*cvO)^2,M[s]+0.01)# variance in observed count (obs error) - if mean = var then poisson 
    R[s] <- max(.0000000001,M[s]^2/(varO[s]-M[s])) # R param for neg binomial with appropriate mean and variance 
    P[s] <- max(.0000000001,1-(varO[s]-M[s])/varO[s]) # P param for neg binomial with appropriate mean and variance 	
    ObsA[s] ~ dnegbin(P[s],R[s])   # observed count adults, reflecting negative binomial observer error
	ObsNew[s]~ dnegbin(P[s],R[s])  # New random data
	res[s] <- (ObsA[s] - M[s])/sqrt(varO[s]) # Pearson residual, observed
	res.new[s] <- (ObsNew[s] - M[s])/sqrt(varO[s]) # Pearson residual
	D[s] <- pow(res[s], 2)          # Deviance, observed data
    D.new[s] <- pow(res.new[s], 2)  # Deviance, new data
  }
  # Test statistics for posterior predictive check
  fit <- sum(D[1:nsurv])         # Summed deviance, observed data
  fit.new <- sum(D.new[1:nsurv]) # Summed deviance, new data
  # Process model
  for (t in 2:maxyr) {
	Dloglam[t-1] <- r*(1 - n[t-1]/K)
	loglam[t-1] ~ dnorm(Dloglam[t-1],tauP) # Stochastic lambda for year t, reflecting process error
	lambda[t-1] <- exp(loglam[t-1])
	n[t] <- n[t-1]*lambda[t-1]
  }
  n[1] <- n0   
  p[1] ~ dbin(pupratio,round(n[1]))
  # Prior model
  sdP ~ dt(0, 1/1^2, 1) T(0.1,) # process error, half Cauchy prior
  tauP <- 1/sdP^2
  cvO ~ dt(0, 1/1^2 ,1) T(0,) # observer error (CV), half Cauchy prior
  r ~ dnorm(.2,100) T(0.1,)  # rmax = maximum intrinsic growth rate, informed normal prior  
  K ~ dt(0, 1/100^2 ,1) T(10,500)    # Overall average K, vague cauchy prior, ~ 10-500
  n0 ~ dt(0, 1/30^2 ,1) T(5,50)      # Starting pop 2000 - vague cauchy prior ~ 5-50
}