model{

  ### MAIN LOOP ###

  # loop through colonies
  for(n in 1:N)
  {
    # loop through years
    for(t in 5:(nYrs-1))
    {

      # marine density dependence determined by SST
      log(a1[n,t]) <- al0 + SST_a1*SST[n,t]

      # recruitment either marine regulated (a2>a1) or terrestrially regulated (a1>a2)
      logit(r[n,t])<-100-max(a1[n,t],a2[n])*P[n,t]

      # fecundity determined by nsAT and prec
      logit(b[n,t])<-a0 + eps[t] + nsAT_b*nsAT[n,t] - nsAT_b2*(nsAT[n,t])^2 + prec_b*prec[n,t]

      # number of young
      y[n,t]<-sj*b[n,t-4]/2*P[n,t-4]*(aH^harvest[n,t-4])

      # weight of redistribution
      w[n,t]<-P[n,t]

      # redistribution function
      W[n,t]<-w[n,t]/wtot[t]
      # number of recruits
      R[n,t]<-r[n,t]*((1-aI)*y[n,t]+aI*ytot[t]*W[n,t])

      # COMMENT-out in future scenario model running
      lambda[n,t]<-sa*P[n,t]+R[n,t]

      # COMMENT-in accounting for HPAIocc in future scenario model running
      #lambda[n,t]<-HPAIocc[n,t]*sa*P[n,t]+HPAIs[n,t]*P[n,t]+R[n,t]

      # population at next time step
      P[n,t+1]~dpois(lambda[n,t] * h[n,t])

      #overdispersion parameter
      h[n,t] ~ dgamma(theta,theta)

      # observed population at next time step
      Pd[n,t+1]~dnorm(P[n,t+1],1/(0.05*P[n,t+1]+1)^2)

      # marine carrying capacity per col per t
      K1[n,t]<-1/a1[n,t]*(100-log(rstar/(1-rstar)))

    }
  }

  for(t in 5:(nYrs-1))
  {
    # total young
    ytot[t]<-sum(y[,t])
    # total extant colonies
    wtot[t]<-sum(w[,t])
    # total available recruits
    rtot[t]<-sum(R[,t])

  }

  # introducing fec data
  # run though all fecundity observations
  for (i in 1:fec_obvs){

    # no. fleglings determined by fecundity (b) at a certain colony nf at a cetain time tf
    # and number of nests
    fledg[i]~dbin(b[nf[i],tf[i]], nests[i])

  }

  # equalise recruitment at t=4
  for(n in 1:N) {R[n,4]<-1}

  # Population time series for each colony
  P1<-P[1,1:nYrs]
  P2<-P[2,1:nYrs]
  P3<-P[3,1:nYrs]
  P4<-P[4,1:nYrs]
  P5<-P[5,1:nYrs]
  P6<-P[6,1:nYrs]
  P7<-P[7,1:nYrs]
  P8<-P[8,1:nYrs]
  P9<-P[9,1:nYrs]
  P10<-P[10,1:nYrs]
  P11<-P[11,1:nYrs]
  P12<-P[12,1:nYrs]
  P13<-P[13,1:nYrs]
  P14<-P[14,1:nYrs]
  P15<-P[15,1:nYrs]
  P16<-P[16,1:nYrs]
  P17<-P[17,1:nYrs]
  P18<-P[18,1:nYrs]
  P19<-P[19,1:nYrs]
  P20<-P[20,1:nYrs]
  P21<-P[21,1:nYrs]
  P22<-P[22,1:nYrs]
  P23<-P[23,1:nYrs]
  P24<-P[24,1:nYrs]
  P25<-P[25,1:nYrs]
  P26<-P[26,1:nYrs]
  P27<-P[27,1:nYrs]
  P28<-P[28,1:nYrs]
  P29<-P[29,1:nYrs]
  P30<-P[30,1:nYrs]
  P31<-P[31,1:nYrs]
  P32<-P[32,1:nYrs]
  P33<-P[33,1:nYrs]
  P34<-P[34,1:nYrs]
  P35<-P[35,1:nYrs]
  P36<-P[36,1:nYrs]
  P37<-P[37,1:nYrs]
  P38<-P[38,1:nYrs]
  P39<-P[39,1:nYrs]
  P40<-P[40,1:nYrs]
  P41<-P[41,1:nYrs]
  P42<-P[42,1:nYrs]
  P43<-P[43,1:nYrs]
  P44<-P[44,1:nYrs]
  P45<-P[45,1:nYrs]
  P46<-P[46,1:nYrs]
  P47<-P[47,1:nYrs]
  P48<-P[48,1:nYrs]
  P49<-P[49,1:nYrs]
  P50<-P[50,1:nYrs]
  P51<-P[51,1:nYrs]
  P52<-P[52,1:nYrs]
  P53<-P[53,1:nYrs]

  # Marine carrying capacity time series for each colony
  K11<-K1[1,5:(nYrs-1)]
  K12<-K1[2,5:(nYrs-1)]
  K13<-K1[3,5:(nYrs-1)]
  K14<-K1[4,5:(nYrs-1)]
  K15<-K1[5,5:(nYrs-1)]
  K16<-K1[6,5:(nYrs-1)]
  K17<-K1[7,5:(nYrs-1)]
  K18<-K1[8,5:(nYrs-1)]
  K19<-K1[9,5:(nYrs-1)]
  K110<-K1[10,5:(nYrs-1)]
  K111<-K1[11,5:(nYrs-1)]
  K112<-K1[12,5:(nYrs-1)]
  K113<-K1[13,5:(nYrs-1)]
  K114<-K1[14,5:(nYrs-1)]
  K115<-K1[15,5:(nYrs-1)]
  K116<-K1[16,5:(nYrs-1)]
  K117<-K1[17,5:(nYrs-1)]
  K118<-K1[18,5:(nYrs-1)]
  K119<-K1[19,5:(nYrs-1)]
  K120<-K1[20,5:(nYrs-1)]
  K121<-K1[21,5:(nYrs-1)]
  K122<-K1[22,5:(nYrs-1)]
  K123<-K1[23,5:(nYrs-1)]
  K124<-K1[24,5:(nYrs-1)]
  K125<-K1[25,5:(nYrs-1)]
  K126<-K1[26,5:(nYrs-1)]
  K127<-K1[27,5:(nYrs-1)]
  K128<-K1[28,5:(nYrs-1)]
  K129<-K1[29,5:(nYrs-1)]
  K130<-K1[30,5:(nYrs-1)]
  K131<-K1[31,5:(nYrs-1)]
  K132<-K1[32,5:(nYrs-1)]
  K133<-K1[33,5:(nYrs-1)]
  K134<-K1[34,5:(nYrs-1)]
  K135<-K1[35,5:(nYrs-1)]
  K136<-K1[36,5:(nYrs-1)]
  K137<-K1[37,5:(nYrs-1)]
  K138<-K1[38,5:(nYrs-1)]
  K139<-K1[39,5:(nYrs-1)]
  K140<-K1[40,5:(nYrs-1)]
  K141<-K1[41,5:(nYrs-1)]
  K142<-K1[42,5:(nYrs-1)]
  K143<-K1[43,5:(nYrs-1)]
  K144<-K1[44,5:(nYrs-1)]
  K145<-K1[45,5:(nYrs-1)]
  K146<-K1[46,5:(nYrs-1)]
  K147<-K1[47,5:(nYrs-1)]
  K148<-K1[48,5:(nYrs-1)]
  K149<-K1[49,5:(nYrs-1)]
  K150<-K1[50,5:(nYrs-1)]
  K151<-K1[51,5:(nYrs-1)]
  K152<-K1[52,5:(nYrs-1)]
  K153<-K1[53,5:(nYrs-1)]

  ### PRIORS ###

  # adult survival
  sa0~dbeta(12,12)
  samax<-0.918+(2*0.023)
  samin<-0.918-(2*0.023)
  sa<-samin+sa0*(samax-samin)

  # juvenile survival
  sj0~dbeta(12,12)
  sjmax<-0.279+(2*0.05)
  sjmin<-0.279-(2*0.05)
  sj<-sjmin+sj0*(sjmax-sjmin)

  # define priors for HPAI matrix (for colonies wher HPAI impact not known)
  # COMMENT-in for running future scenario model
  #HPAIs[1,123]~dunif(0.455,samax)
  #HPAIs[2,123]~dunif(0.455,samax)
  #HPAIs[12,123]~dunif(0.455,samax)
  #HPAIs[15,123]~dunif(0.455,samax)
  #HPAIs[19,123]~dunif(0.455,samax)
  #HPAIs[23,123]~dunif(0.455,samax)
  #HPAIs[38,123]~dunif(0.455,samax)
  #HPAIs[47,123]~dunif(0.455,samax)
  #HPAIs[49,123]~dunif(0.455,samax)
  #HPAIs[50,123]~dunif(0.455,samax)

  # immigration prior
  aI0~dbeta(12,12)
  aImax<-0.5
  aImin<-0.3
  aI<-aImin+aI0*(aImax-aImin)

  # harvest
  aH~dbeta(1,1)

  # overdispersion prios
  theta <- max(10000-theta0,1000)
  theta0 ~ dexp(1/200)

  # COVARIATE PRIORS
  prec_b~dnorm(0,1)     # prior for prec in b (fec)
  nsAT_b~dgamma(1,1)   # prior for linear nsAT in b (fec)
  nsAT_b2~dgamma(1,1)  # prior for non-linear nsAT in b (fec)
  SST_a1~dnorm(0,1)    # prior for SST in a1 (fec)

  # Initialise first 5 years of model
  for(n in 1:N){
    for (t in 1:4){
      y[n,t]~dunif(0,5000)
      logit(b[n,t])<-a0 + eps[t]
      w[n,t]<-min(1,P[n,t])
    }
  }

  for (n in 1:N){

    for (t in 1:5){
      P[n,t]~dunif(0,5000)
      Pd[n,t]~dnorm(P[n,t],1/(0.05*P[n,t]+1)^2)

    }
  }

  # intercept of density dependence (a1)
  # mean 0.05, sd 0.02
  a10~dgamma(6.25,125)
  al0<-log(a10)


  # intercept of fecundity
  a0~dnorm(0.84,1)
  # baseline fecundity
  logit(bstar)<-a0
  # define constant
  rstar<-(1-sa)/sj*bstar

  # priors for terrestrial carrying capacity and density dependence
  for(n in 1:N)
  {
    K0[n]~dbeta(9.2,13.8)
    Kmax[n]<-1.2*cc[n]
    Kmin[n]<-0.8*cc[n]
    K[n]<-Kmin[n]+K0[n]*(Kmax[n]-Kmin[n])

    a2[n]<-1/K[n]*(100-log(rstar/(1-rstar)))

  }

  # stochastic annual fecundity fluctuations in each year
  for(t in 1:(nYrs))
  {

    eps[t]~dnorm(0,preceps)
  }
  # precision of eps
  preceps~dgamma(10,1)


  #inits# a0

  #data# nYrs,N,Pd,harvest, cc, fledg, nests, nf, tf, fec_obvs, SST, prec, nsAT

  #extra data needed for running future scenarios model:
  # HPAIocc, HPAIs

  #monitor model coefficients historical
  #monitor# a0, al0, a2, aH, aI, sa, sj, K, theta, preceps, nsAT_b, nsAT_b2, prec_b, SST_a1

  #monitor populations sizes
  #monitor# P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27,P28,P29,P30,P31,P32,P33,P34,P35,P36,P37,P38,P39,P40,P41,P42,P43,P44,P45,P46,P47,P48,P49,P50,P51,P52,P53

  #monitor carrying capacities
  #monitor# K11,K12,K13,K14,K15,K16,K17,K18,K19,K110,K111,K112,K113,K114,K115,K116,K117,K118,K119,K120,K121,K122,K123,K124,K125,K126,K127,K128,K129,K130,K131,K132,K133,K134,K135,K136,K137,K138,K139,K140,K141,K142,K143,K144,K145,K146,K147,K148,K149,K150,K151,K152,K153

}