model;
{

# Process equations (Catch at beginning of year)
########################

Pmed[1] <-0

for (i in 2:21) {
   	Pmed[i]<-log(max(exp(-m[i])*(rho[i-1]+alpha[i-1]/w[i-1])*(P[i-1]-Catch[i-1]/K)+exp(-m[i])*(rho[i-1]+alpha[i-1]/wr[i-1])*rp[i-1],0.01))
    }


#The following adjusts June 2004 to September 2004. Survey in June from 1983 to 2003 and then moved to August in 2004 and 2005

	Pmed[22]<-log(max(exp(-1.00*m[22])*(rho[21]+alpha[21]/w[21])*(rho2004+alpha2004/w2004)*(P[21]-Catch[21]/K)+exp(-1.00*m[22])*(rho[21]+alpha[21]/wr[21])*(rho2004+alpha2004/wr2004)*rp[21],0.01))
	
	Pmed[NY]<-log(max(exp(-m[NY])*(rho[NY-1]+alpha[NY-1]/w[NY-1])*(P[NY-1]-Catch[NY-1]/K)+exp(-m[NY])*(rho[NY-1]+alpha[NY-1]/wr[NY-1])*rp[NY-1],0.01))
	
for(i in 1:NY){
	P[i]~dlnorm(Pmed[i],isigma)I(0,5)
}


#Observation equations
#########################

#Survey Biomass (shell height >= 80mm)

for(i in 1:NY){
	Imed[i]<-log(qI[i]*K*P[i])
		I[i]~dlnorm(Imed[i],ivarepsilon)
}

#Survey numbers  (shell height >= 80mm)

for(i in 1:NY){
	Lmed[i]<-log(qI[i]*K*P[i]/(w[i]/pow(10,6)))
		L[i]~dlnorm(Lmed[i],itau)
}


#Survey recruits (65<shell height<80 mm)

for(i in 1:NY){
		rp[i]~dlnorm(-1.9,0.5)I(0,6)
}

for(i in 1:NY){
	Rmed[i]<-log(ratiolined[i]*qI[i]*K*rp[i])
		R[i]~dlnorm(Rmed[i],inu)
}


#S=dissolution rate in days/365

S~dunif(0.10,0.99)

for(i in 1:NY){
		m[i]~dlnorm(-1.9,0.5)I(0,4)
}


#Natural mortality from survey clapper index, revised for popcorn model (Smith and Lundy 2002)

cmed[1]<-log(m[1]*S*L[1])

for(i in 2:NY){
		cmed[i]<-log(m[i]*S*(S*L[i-1]+(2-S)*L[i])/2)
		}

for(i in 1:NY){
clappers[i]~dlnorm(cmed[i],iepsilon)
}


#Distribution of K
###########################


K~dlnorm(8.006,0.63390)I(500,10000)

#Distribution of qI annually, catchability coefficient for scallops with shell height >= 80 mm
#Changed to annual qI because of patterns in the predicted survey/population biomass. See Smith et al. (2005)
###########################


for (i in 1:NY){

qI[i]~dbeta(1.0,1.0)

}

#Distribution of variance terms
###########################

#Process error
isigma~dgamma(3,0.44639)
sigma<-1/isigma

#Error for survey biomass (shell height >= 80mm)
ivarepsilon~dgamma(3,0.44639)
varepsilon<-1/ivarepsilon


#Error for survey numbers (shell height >= 80mm)
itau~dgamma(3,0.44639)
tau<-1/itau

#Error for recruitment biomass
inu~dgamma(3,0.89258)
nu<-1/inu

#Error for clappers
iepsilon~dgamma(3,0.89258)
epsilon<-1/iepsilon

# Output
############################

#Predicted biomass (shell height >= 80mm)

for(t in 1:NY){
Biomass[t]<-P[t]*K
}

#Predicted survey index (shell height >= 80mm)
for(t in 1:NY){
Ipred[t]<-P[t]*K*qI[t]
}

#Predicted recruitment biomass (65< Shell height<80)
#ratiolined is the ratio of catches from the lined and unlined survey gear.  See Smith et al (2005).

for(t in 1:NY){
Rec[t]<-rp[t]*K
IRec[t]<-rp[t]*K*ratiolined[t]*qI[t]
}


#Diagnostics (Raw and standardized residuals)
#Sample from posterior to check extreme values.
#############################################

for(i in 1:NY){
	resid.p[i]<-log(P[i])-Pmed[i]
	sresid.p[i]<-resid.p[i]*sqrt(isigma)
	resid.I[i]<-log(I[i])-Imed[i]
	sresid.I[i]<-resid.I[i]*sqrt(ivarepsilon)
	resid.R[i]<-log(R[i])-Rmed[i]
	sresid.R[i]<-resid.R[i]*sqrt(inu)
	resid.c[i]<-log(clappers[i])-cmed[i]
	sresid.c[i]<-resid.c[i]*sqrt(iepsilon)
	resid.L[i]<-log(L[i])-Lmed[i]
	sresid.L[i]<-resid.L[i]*sqrt(itau)
	}


for(i in 1:NY){
	       I.rep[i]~dlnorm(Imed[i],ivarepsilon)
		p.I.smaller[i]<-step(I[i]-I.rep[i])
	       L.rep[i]~dlnorm(Lmed[i],itau)
		p.L.smaller[i]<-step(L[i]-L.rep[i])
       	R.rep[i]~dlnorm(Rmed[i],inu)
		p.R.smaller[i]<-step(R[i]-R.rep[i])
       clappers.rep[i]~dlnorm(cmed[i],iepsilon)
		p.clappers.smaller[i]<-step(clappers[i]-clappers.rep[i])
}
	


################
#Projections
################


exploit.this.year.1<-Catch[NY]/Biomass[NY]
Pr.exploit.1<-step(exploit.this.year.1-0.20)

P.nextyeara<-log(max(exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/w[NY])*(P[NY]-Catch[NY]/K)+exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/wr[NY])*rp[NY],0.01))

P.nextyear~dlnorm(P.nextyeara,isigma)I(0,6)

B.nextyear<-P.nextyear*K

survey.next.year.1<-P.nextyear*K*qI[NY]



################################
#Catch in 2005/2006 = 100 t

exploit.next.year.1a<-50/B.nextyear
Pr.exploit.next.1a<-step(exploit.next.year.1a-0.20)

exploit.next.year.1<-100/B.nextyear
Pr.exploit.next.1<-step(exploit.next.year.1-0.20)


exploit.next.year.2<-150/B.nextyear
Pr.exploit.next.2<-step(exploit.next.year.2-0.20)

exploit.next.year.3<-200/B.nextyear
Pr.exploit.next.3<-step(exploit.next.year.3-0.20)

exploit.next.year.4<-250/B.nextyear
Pr.exploit.next.4<-step(exploit.next.year.4-0.20)

###################################

#Catch in 2005/2006 equals 200 t

exploit.this.year.2<-200/Biomass[NY]
Pr.exploit.2<-step(exploit.this.year.2-0.20)

P.nextyear.1b<-log(max(exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/w[NY])*(P[NY]-500/K)+exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/wr[NY])*rp[NY],0.01))

P.nextyear.1~dlnorm(P.nextyear.1b,isigma)I(0,6)

B.nextyear.1<-P.nextyear.1*K

exploit.next.year.5a<-50/B.nextyear.1
Pr.exploit.next.5a<-step(exploit.next.year.5a-0.20)

exploit.next.year.5<-100/B.nextyear.1
Pr.exploit.next.5<-step(exploit.next.year.5-0.20)


exploit.next.year.6<-150/B.nextyear.1
Pr.exploit.next.6<-step(exploit.next.year.6-0.20)

exploit.next.year.7<-200/B.nextyear.1
Pr.exploit.next.7<-step(exploit.next.year.7-0.20)

exploit.next.year.8<-250/B.nextyear.1
Pr.exploit.next.8<-step(exploit.next.year.8-0.20)

###################################

#Catch in 2005/2006 equals 250 t

exploit.this.year.3<-250/Biomass[NY]
Pr.exploit.3<-step(exploit.this.year.3-0.20)

P.nextyear.1c<-log(max(exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/w[NY])*(P[NY]-600/K)+exp(-1.0*m[NY])*(rho[NY]+alpha[NY]/wr[NY])*rp[NY],0.01))

P.nextyear.2~dlnorm(P.nextyear.1c,isigma)I(0,6)

B.nextyear.2<-P.nextyear.2*K

exploit.next.year.9a<-50/B.nextyear.2
Pr.exploit.next.9a<-step(exploit.next.year.9a-0.20)

exploit.next.year.9<-100/B.nextyear.2
Pr.exploit.next.9<-step(exploit.next.year.9-0.20)


exploit.next.year.10<-150/B.nextyear.2
Pr.exploit.next.10<-step(exploit.next.year.10-0.20)

exploit.next.year.11<-200/B.nextyear.2
Pr.exploit.next.11<-step(exploit.next.year.11-0.20)

exploit.next.year.12<-250/B.nextyear.2
Pr.exploit.next.12<-step(exploit.next.year.12-0.20)

}