model
{

	for (s in 1:3) {

      q[s] ~ dlnorm(0, 100)  # Informative Priors equivalent to q=1 (1/trawlable_units) with sd=0.1

      for (d in 1:3) {

            z[s,d] ~ dunif(0, 3)    # Priors for instantaneous mortality rates with decadal variation

            } 

      sdpro[s] ~ dunif(0, 5) # Process error variance

      isigma2[s] <- pow(sdpro[s], -2)

      }

# Priors for observation error variance

sdobs[1]~dunif(0.497,2)

sdobs[2]~dunif(0.586,2)

sdobs[3]~dunif(0.513,2)

for(i in 1:3){

    itau2[i] <- pow(sdobs[i], -2)

      }

# Priors (informative) for transistion probabilities

theta[1] ~ dbeta(20, 80)

theta[2] ~ dbeta(25, 75)

# Prior for recruitment rate

recrate ~ dlnorm(1.178655, 4)I(0.01,1000)

# Informative priors for inital abundance (based on trawlable abundance in 1000s)

# sd=20

# Nad[1] is adult abundance in 1965

N69[1] ~ dlnorm(0, 6.25)

N69[2] ~ dlnorm(-1.203973, 6.25)

for (i in 1:5) {Nad[i] ~ dlnorm(1.252763, 11.11111)}

# Process equation

for (s in 1:3) {

      for (t in 1:10 ) {

            st[s,t] <- exp(-z[s,1])

            st[s,t+10] <- exp(-z[s,2])

            st[s,t+20] <- exp(-z[s,3])

            }

      for(t in 31:35){

            st[s,t] <- exp(-z[s,3])

            }

      }

nmod[1,1] <- (Nad[1] * (recrate / 2) + N69[1] * (1 - theta[1]))* st[1,1]

nmod[2,1] <- (N69[1] * theta[1] + N69[2] * (1 - theta[2]))* st[2,1]

nmod[3,1] <- (N69[2] * theta[2] + Nad[5])* st[3,1]

for (s in 1:3) {

      logn[s,1] <- log(nmod[s,1])

      N[s,1] ~ dlnorm(logn[s,1], isigma2[s])

      }

nmod[1,2] <- (Nad[2] * (recrate / 2) + N[1,1]  * (1 - theta[1]))* st[1,2]

nmod[2,2] <- (N[1,1]  * theta[1] + N[2,1] * (1 - theta[2]))* st[2,2]

nmod[3,2] <- (N[2,1]  * theta[2] + N[3,1])* st[3,2]

for (s in 1:3) {

      logn[s,2] <- log(nmod[s,2])

      N[s,2] ~ dlnorm(logn[s,2], isigma2[s])

      }

for (t in 3:5) {

      nmod[1,t] <- (Nad[t] * (recrate / 2) + N[1,t-1]  * (1 - theta[1]))* st[1,t]

      nmod[2,t] <- (N[1,t-1]  * theta[1] + N[2,t-1]  * (1 - theta[2]))* st[2,t]

      nmod[3,t] <- (N[2,t-1]  * theta[2] + N[3,t-1])* st[3,t]

      for (s in 1:3) {

            logn[s,t] <- log(nmod[s,t])

            N[s,t] ~ dlnorm(logn[s,t], isigma2[s])

            }

      }

for (t in 6:35 ) {

      nmod[1,t] <- (N[3,t-5] * (recrate / 2) + N[1,t-1]  * (1 - theta[1]))* st[1,t]

      nmod[2,t] <- (N[1,t-1]  * theta[1] + N[2,t-1]  * (1 - theta[2]))* st[2,t]

      nmod[3,t] <- (N[2,t-1]  * theta[2] + N[3,t-1])* st[3,t]

      for (s in 1:3) {

            logn[s,t] <- log(nmod[s,t])

            N[s,t] ~ dlnorm(logn[s,t], isigma2[s])

            }

      }

# Observation equation

for (t in 1:35 ) {

            for (s in 1:3) {

            Im[s,t] <- q[s] * N[s,t]

            logIm[s,t] <- log(Im[s,t])

            Iobs[s,t] ~ dlnorm(logIm[s,t], itau2[s])

            #resid[s,t] <- Iobs[s,t] - Im[s,t]

            }

      }

} # End model