model
{
	# N.level <- 3
	# v.diag <- rep(1.0, N.level)
	# R = diag(v.diag, N.level),
	# Deg.freedom = 3,
	# Tau.noninformative = 1.0e-2,
	# Hyper.gamma = 1.0e-3
	for (i in 1:N.site) {
		# Observation
		Fish1[i] ~ dcat(p.fish[i,])
		BiomassZP[i] ~ dnorm(bzp[i], tau.err[1]) 
		Chl.a[i] ~ dnorm(chl.a[i], tau.err[2])
		log.din[i] ~ dnorm(log.DIN[i], tau.err[3])
		log.dtp[i] ~ dnorm(log.DTP[i], tau.err[4])
		w.temp[i] ~ dnorm(W.temp[i], tau.err[5]) 
		# log.density -> mean[i]
		p.fish[i,1] <- p.fish.zero[i]
		p.fish[i,2] <- 1.0 - p.fish.zero[i]
		p.fish.zero[i] <- exp(-exp(log.fp[i,1]))
		log(bzp[i]) <- log.fp[i,2]
		log(chl.a[i]) <- log.fp[i,3]
		log.fp[i,1:N.level] ~ dmnorm(mu[i,], inv.vc[,])
		# Means
		mu[i,1] <- ( # log(fish density)
			beta.f[1]
			+ beta.f[2] * w.temp[i]
			+ beta.f[3] * log.Area[i]
			+ beta.f[4] * Altitude[i]
		)
		mu[i,2] <- ( # log(biomass of zoo plankton)
			beta.z[1]
			+ beta.z[2] * w.temp[i]
		)
		mu[i,3] <- ( # log(density of Chl-a)
			beta.c[1]
			+ beta.c[2] * w.temp[i]
			+ beta.c[3] * log.din[i]
			+ beta.c[4] * log.dtp[i]
		)
	}
	# parameters
	for (k in 1:N.beta.f) {
		beta.f[k] ~ dnorm(0, Tau.noninformative)
	}
	for (k in 1:N.beta.z) {
		beta.z[k] ~ dnorm(0, Tau.noninformative)
	}
	for (k in 1:N.beta.c) {
		beta.c[k] ~ dnorm(0, Tau.noninformative)
	}
	# hyper priors
	for (k in 1:N.err) {
		tau.err[k] ~ dgamma(Hyper.gamma, Hyper.gamma)
	}
	inv.vc[1:N.level,1:N.level] ~ dwish(R[,], Deg.freedom)
	vc[1:N.level,1:N.level] <- inverse(inv.vc[,])
	for (l1 in 1:N.level) {
		for (l2 in 1:N.level) {
			rho[l1, l2] <- vc[l1,l2] / sqrt(vc[l1, l1] * vc[l2, l2])
		}
	}
}