Commit 38b5097d authored by Berges, Benoit's avatar Berges, Benoit
Browse files

update scripts and add no transfer scripts

parent 6a9164b2
rm(list=ls())
library(icesTAF)
load('model/NSAS_sanoba_200.RData')
metricsPeriod <- projPeriod[(length(projPeriod)-10):(length(projPeriod)-1)]
......@@ -21,10 +23,10 @@ SSB_bool <- array(FALSE,dim=c(1,nits))
for(idxIter in 1:nits){
# store value per year
SSB_riskMat[which(SSB[,idxIter] < referencePoints$Blim),idxIter] <- TRUE
SSB_riskMat[which(SSB[,idxIter] <- referencePoints$Blim),idxIter] <- TRUE
# TRUE/FALSE for each iteration
if(length(which(SSB[,idxIter] < referencePoints$Blim)!=0))
if(length(which(SSB[,idxIter] <- referencePoints$Blim)!=0))
SSB_bool[idxIter] <- TRUE
}
......@@ -41,21 +43,7 @@ colnames(lmFrame) <- c('years','risk','SSB5per')
riskLR <- lm(risk~years,data=lmFrame)
SSBTailLR <- lm(SSB5per~years,data=lmFrame)
# plot
par(mfrow=c(2,1))
plot(metricsPeriod,SSBTail[2,],type='l',
ylim=c(0,max(SSBTail[3,])),
ylab='SSB',
xlab='years',main=paste0('ftar_',ftarget,'_btrig_',btrigger))
lines(metricsPeriod,SSBTail[1,],type='l',col='blue',lty=2)
lines(metricsPeriod,SSBTail[3,],type='l',col='blue',lty=2)
plot(metricsPeriod,SSB_prob,type='l',
ylim=c(0,0.2),
ylab='Annual risk',
xlab='years',
xlim=c(an(min(metricsPeriod)),an(max(metricsPeriod))))
lines(c(2010,2040),c(0.05,0.05),lty=2,col='green')
# Fbar
FHCR26 <- FHCR[f26,ac(2021:2041)]
......@@ -78,4 +66,37 @@ LTF26 <- mean(f26Quant['50%',])
LT_RiskTrend <- summary(riskLR)$coefficients[2,1]
LT_RiskTrendSig <- ifelse(summary(riskLR)$coefficients[2,4] > 0.05,0,1) # 0 is not significant. 1 is significant
LT_SSBTailTrend <- summary(SSBTailLR)$coefficients[2,1]
LT_SSBTailTrendSig <- ifelse(summary(SSBTailLR)$coefficients[2,4] > 0.05,0,1) # 0 is not significant. 1 is significant
\ No newline at end of file
LT_SSBTailTrendSig <- ifelse(summary(SSBTailLR)$coefficients[2,4] > 0.05,0,1) # 0 is not significant. 1 is significant
perfMetrics <- cbind(LTR1,LTR3,LTY,LTF01,LTF26)
write.csv(x=perfMetrics,file = 'temp.csv')
# plotting
taf.png("NSAS_risk_ssb")
# plot
par(mfrow=c(2,1))
plot(metricsPeriod,SSBTail[2,],type='l',
ylim=c(0,max(SSBTail[3,])),
ylab='SSB',
xlab='years',main=paste0('ftar_',ftarget,'_btrig_',btrigger))
lines(metricsPeriod,SSBTail[1,],type='l',col='blue',lty=2)
lines(metricsPeriod,SSBTail[3,],type='l',col='blue',lty=2)
plot(metricsPeriod,SSB_prob,type='l',
ylim=c(0,0.2),
ylab='Annual risk',
xlab='years',
xlim=c(an(min(metricsPeriod)),an(max(metricsPeriod))))
lines(c(2010,2040),c(0.05,0.05),lty=2,col='green')
dev.off()
taf.png("NSAS_risk_ssb")
taf.png("MSE_projections")
plot(biol)
dev.off()
taf.png("MSE_projections_zoomed")
plot(window(biol,start=2000))
dev.off()
\ No newline at end of file
#-------------------------------------------------------------------------------
# 1) load packages
# setup paths
# load functions
#-------------------------------------------------------------------------------
rm(list=ls())
HCR <- 'FMSYAR'
ftarget <- 0.31
btrigger <- 1232828
cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
library(FLSAM)
library(FLEDA)
library(FLFleet)
library(ggplotFL)
library(minpack.lm) # install.packages("minpack.lm")
library(stats)
runName <- "NSAS_sanoba"
# paths to different subfolders
dataPath <- file.path(".","data/")
modelPath <- file.path(".","model/")
scriptPath <- file.path(".","side_scripts/")
functionPath <- file.path(".","functions/")
source(file.path(functionPath,"MSE_assessment.r"))
source(file.path(functionPath,"projectNSH_FMSYAR_no_transfer.R"))
# ImY forecast functions
source(file.path(functionPath,"fleet.harvest.r"))
source(file.path(functionPath,"rescaleF.r"))
source(file.path(functionPath,"TAC2sel.r"))
# fmsyAR function
source(file.path(functionPath,"fmsyAR_fun_no_transfer.r"))
source(file.path(functionPath,"find.FCAR.r"))
#-------------------------------------------------------------------------------
# 2) Initialize
#
# Define MSE parameters,
# load objects initialized previously
# * Biology
# - biol
# - surveys
# * Fisheries
# - FCProp
# - catchD
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 10
#load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
# load object
load(file.path(modelPath,paste0(runName,'_IBP_MSE_init_',ac(nits),'.RData')))
stkAssessment.ctrl <- NSH.ctrl
# load MSE parameters
load(file.path(modelPath,paste0(runName,'_IBP_MSE_params_',ac(nits),'.RData')))
fishery@landings.sel[,projPeriod] <- sweep(fishery@landings.sel[,projPeriod],c(2:6),quantMeans(fishery@landings.sel[,projPeriod]),"/")
strFleet <- c('A','B','C','D')
nFleets <- length(strFleet)
nAges <- dim(biol@stock.n)[1]
surveyNames <- names(surveys)
escapeRuns <- numeric()
runNameSave <- paste0(runName,"_MSE_",nits)
#-------------------------------------------------------------------------------
# 2) ref points
# Ftarget = FMSY and Btrigger = MSYBtrigger for now
#
#-------------------------------------------------------------------------------
referencePoints <- list(Fmsy = 0.31,
Fsq = NA,
Flim = 0.39,
Fpa = 0.31,
Blim = 874198,
Bpa = 956483,
MSYBtrigger = 1232828,
Ftarget = ftarget,
F01 = 0.05,
Btrigger = btrigger)
managementRule <- list(HCR = HCR)
#------------------------------------------------------------------------------#
# 3) Housekeeping
#------------------------------------------------------------------------------#
CATCH <- TAC
FHCR <- FLQuant(NA, dimnames=list(age=0:8,year=ac(an(projPeriod[1]):(an(projPeriod[length(projPeriod)])+3)),unit="unique",season="all",area="unique",iter=1:nits))
SSBHCR <- FLQuant(NA, dimnames=list(age="all",year=ac(an(projPeriod[1]):(an(projPeriod[length(projPeriod)])+3)),unit="unique",season=c("FcY","CtY"),area="unique",iter=1:nits))
#------------------------------------------------------------------------------#
# 4) Start running the MSE
#------------------------------------------------------------------------------#
start.time <- Sys.time()
#an(projPeriod)
for (iYr in 2021:2042){
cat(iYr,"\n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#----------------------------------------
# define year names
#----------------------------------------
TaY <- ac(iYr-1) # terminal year in the assessment
ImY <- ac(iYr) # intermediate year in the short term forecast, ie, current year
FcY <- ac(iYr+1) # year for which the advice is given, ie, forecast two years ahead the last year in the assessment
FuY <- c(ImY,FcY) # combination of future years
#----------------------------------------
# update the biol number at age in ImY
#----------------------------------------
#- Define mortality rates for iYr-1 to calculate survivors to iYr
m <- m(biol)[,ac(iYr-1),,]
z <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,]) + m
#- Update biological model to iYr
survivors <- stock.n(biol)[,ac(iYr-1)] * exp(-z)
stock.n(biol)[ac((range(biol,"min")+1):range(biol,"max")),ac(iYr),,] <- survivors[-dim(survivors)[1],,,,,]@.Data
biol@harvest[,ac(iYr-1)] <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,])
#- Update recruitment
#biol.sr@rec <- rec(biol[,ac(iYr-1)])
biol.sr@ssb[,ac(iYr-1)] <- ssb(biol[,ac(iYr-1)])
recruitBio <- predict(biol.sr[,ac(iYr-1)]) * exp(sr.res[,ac(iYr),drop=T])
stock.n(biol)[1,ac(iYr)] <- recruitBio
#- Plusgroup
if (!is.na(range(biol,"plusgroup"))){
stock.n(biol)[ac(range(biol,"max")),ac(iYr),] <- stock.n(biol)[ac(range(biol,"max")),ac(iYr),] + survivors[ac(range(biol,"max"))]
}
#- Apply process error per cohort age 1 to 8
for(idxAge in 2:nAges){
biol@stock.n[idxAge,ac(iYr)] <- biol@stock.n[idxAge,ac(iYr)]*varProccError[idxAge-1,ac(an(ac(iYr))-(idxAge-1))]
}
cat("\n Finished biology \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#- Update fishery to year iYr-1
landings.n(fishery)[,ac(iYr-1)] <- sweep(sweep(landings.sel(fishery)[,ac(iYr-1),,,,],c(1:4,6),z,"/"),c(1:4,6),stock.n(biol)[,ac(iYr-1)]*(1-exp(-z)),"*")
catch.n(biol)[,ac(iYr-1)] <- areaSums(fishery@landings.n[,ac(iYr-1)])
#print(computeLandings(fishery)[,ac(iYr-1)]/TAC[,ac(iYr-1)])
#-------------------------------------------------------------------------------
# Assessment
#-------------------------------------------------------------------------------
# filter stock object up to intermediate year to include biological variables
stkAssessment <- window(biol,end=an(TaY))
stkAssessment@catch.n <- stkAssessment@catch.n * catchVar[,ac(histMinYr:TaY),,'residuals']
stkAssessment@landings.n <- stkAssessment@catch.n
stkAssessment@landings <- computeLandings(stkAssessment)
#stkAssessment <- window(biol,end=an(TaY))
# smooth M prior to running the assessment, median filter of order 5
require(doParallel); ncores <- detectCores()-1; ncores <- ifelse(nits<ncores,nits,ncores);cl <- makeCluster(ncores); registerDoParallel(cl)
dat <- as.data.frame(stkAssessment@m)
datS <- split(dat,as.factor(paste(dat$age,dat$iter)))
res <- foreach(i = 1:length(datS)) %dopar% fitted(loess(data ~ year,data=datS[[i]],span=0.5))
stkAssessment@m <- FLQuant(c(aperm(array(unlist(res),dim=c(length(1947:TaY),nits,nAges)),c(3,1,2))),dimnames=dimnames(stkAssessment@m))
stopCluster(cl)#; detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE); #detach("package:iterators",unload=TRUE)
# update surveys
for(idxSurvey in surveyNames){
agesSurvey <- an(rownames(surveys[[idxSurvey]]@index))
yearSurvey <- an(colnames(surveys[[idxSurvey]]@index))
surveyProp <- mean(c(surveys[[idxSurvey]]@range[6],surveys[[idxSurvey]]@range[7]))
surveys[[idxSurvey]]@index[,TaY] <- surveyVars[ac(agesSurvey),TaY,idxSurvey,'catchabilities']*
exp(-z[ac(agesSurvey),TaY]*surveyProp)*
biol@stock.n[ac(agesSurvey),TaY]*surveyVars[ac(agesSurvey),TaY,idxSurvey,'residuals']
}
stkAssessment.tun <- window(surveys,end=an(TaY)+1)
stkAssessment.tun[["HERAS"]] <- window(surveys[["HERAS"]],end=an(TaY))
stkAssessment.tun[["IBTS-Q3"]] <- window(surveys[["IBTS-Q3"]],end=an(TaY))
stkAssessment.ctrl@range[5] <- an(TaY)+1
stkAssessment.ctrl@residuals <- F
#stk0 <- stkAssessment
#idx0 <- stkAssessment.tun
#sam0.ctrl <- stkAssessment.ctrl
#escapeRuns <- escapeRuns
#resInit <- stkAssessment.init
# Run stock assessment
ret <- MSE_assessment( stkAssessment,
stkAssessment.tun,
stkAssessment.ctrl,
escapeRuns) #stkAssessment.init
stkAssessment.init <- ret$resInit
escapeRuns <- ret$escapeRuns
stkAssessment <- ret$stk
#compObj <- new('FLStocks')
#compObj[['assessment']] <- stkAssessment
#compObj[['biol']] <- biol
#windows()
#plot(compObj)
#print(escapeRuns)
#print(stkAssessment@stock.n[,ac(TaY)] / biol@stock.n[,ac(TaY)])
cat("\n Finished stock assessment \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#-------------------------------------------------------------------------------
# Forecast
#-------------------------------------------------------------------------------
projNSAS <- projectNSH_FMSYAR_no_transfer(stkAssessment,
fishery,
iYr,
TAC,
TAC_var,
histMaxYr,
referencePoints)
TAC[,FcY,,,c("A","B")] <- projNSAS$TAC[,,,,c("A","B")]
FHCR[,FcY] <- projNSAS$Fbar
SSBHCR[,FcY,,"FcY"] <- projNSAS$SSB$FcY #store HCR SSB in the forecast year
SSBHCR[,FcY,,"CtY"] <- projNSAS$SSB$CtY #store HCR SSB in the continuation year
cat("\n Finished forecast \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#-------------------------------------------------------------------------------
# TAC to real F again
#-------------------------------------------------------------------------------
#-Calculate effort accordingly (assuming constant catchability)
#- Get a good starting condition
#mults <- matrix(NA,nrow=nits,ncol=4)
#for(idxIter in 1:nits)
# mults[idxIter,] <- optim(par=runif(4),fn=TAC2sel_V2,iYr=ImY,iBiol=biol[,ImY],iFishery=fishery[,ImY],iTAC=TAC[,ImY],catchVar=catchVar,TAC_var=TAC_var,iTer=idxIter,control=list(maxit=1000),lower=rep(1e-8,4),method="L-BFGS-B")$par
CATCH[,ImY,,,"A"] <- TAC[,ImY,,,"A"] + TAC_var[ImY,,'Ctransfer'] * TAC[,ImY,,,"C"]
CATCH[,ImY,,,"B"] <- TAC[,ImY,,,"B",drop=T] * TAC_var[ImY,,'Buptake',drop=T]
require(doParallel); ncores <- detectCores()-1; ncores <- ifelse(nits<ncores,nits,ncores);cl <- makeCluster(ncores); clusterEvalQ(cl,library(FLCore)); clusterEvalQ(cl,library(minpack.lm)); registerDoParallel(cl)
multso <- do.call(rbind,foreach(idxIter = 1:nits) %dopar% nls.lm(par=runif(4),
fn=TAC2sel,
iYr=ImY,
iBiol=biol[,ImY],
iFishery=fishery[,ImY],
iTAC=CATCH[,ImY],
catchVar=catchVar,
TAC_var=TAC_var,
iTer=idxIter,
control=nls.lm.control(maxiter=1000),
lower=rep(1e-8,4))$par)
stopCluster(cl); detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE)
#catchVar[1,ac(iYr),,"FCprop",,]
#Check for very high F
idx <- which(quantMeans(sweep(landings.sel(fishery)[,ac(ImY)],3:6,t(multso),"*")[ac(2:6),,,,"A"]) > 5)
if(length(idx)>0){
print(idx)
fishery@landings.sel[,ac(ImY),,,,-idx] <- sweep(landings.sel(fishery[,ac(ImY),,,,-idx]),3:6,t(multso)[,-idx],"*")
fishery@landings.sel[,ac(ImY),,,,idx] <- landings.sel(fishery[,ac(an(ImY)-1),,,,idx])
} else {
#When setting a very high F this indicate s that the optimiser doesn't work well, so replace with last year F
fishery@landings.sel[,ac(ImY)] <- sweep(landings.sel(fishery)[,ac(ImY)],3:6,t(multso),"*")
}
cat("\n Finished effort calc \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
}
save.image(file=file.path(modelPath,paste0(runName,'_',nits,".RData")))
fmsyAR_fun_no_transfer <- function( stf,
FuY,
TACS,
RECS,
referencePoints,
TAC_var,
f01,
f26){
FuY,
TACS,
RECS,
referencePoints,
TAC_var,
f01,
f26){
#TACS <- iTAC
#referencePoints <- mpPoints
#TAC_var <- iTAC_var
ImY <- FuY[1]
FcY <- FuY[2]
......@@ -25,9 +29,10 @@ fmsyAR_fun_no_transfer <- function( stf,
CATCH <- TACS
# assume no transfer of the C fleet TAC, i.e. catches only in IIIa
CATCH[,c(FcY,CtY),'C'] <- TACS[,FcY,'C']*TAC_var$Csplit
CATCH[,c(FcY,CtY),'D'] <- TACS[,FcY,'D']*TAC_var$Dsplit*TAC_var$Duptake
# assume a transfer of the C fleet TAC, i.e. catches in IIIa and IVa. The catches in IIIa is CATCH for C
# while catches in IVa are now transferred in the A fleet
CATCH[,c(FcY,CtY),,,'C'] <- TACS[,FcY,,,'C']*TAC_var[FcY,,'Csplit']#*(1-TAC_var$Ctransfer)
CATCH[,c(FcY,CtY),,,'D'] <- TACS[,FcY,,,'D']*TAC_var[FcY,,'Dsplit']*TAC_var[FcY,,'Duptake']
for(iTer in 1:dims(stf)$iter) #stf.=stf,rec.=rec,f.=fmsy,f26.=f26,f01.=f01,TACS.=TACS
res[,iTer] <- nls.lm(par=rep(1,3), # A scalor = B scalor, therefore 3 scalors
......@@ -42,11 +47,53 @@ fmsyAR_fun_no_transfer <- function( stf,
jac=NULL,nls.lm.control(ftol = (.Machine$double.eps),
maxiter = 1000))$par
res <- c(res[1],res[1],res[2],res[3])
# same scalor for B and A fleets
res <- rbind(res[1,],res[1,],res[2,],res[3,])
#iTer <- 1
#find.FCAR.verbose(res[,iTer],
# iter(stf[,FcY],iTer),
# iter(CATCH[,FcY],iTer),
# referencePoints,
# f01,
# f26)
#print(res)
# update F with scalors
stf@harvest[,FcY,c('A','B','C','D')] <- sweep(stf@harvest[,FcY,c('A','B','C','D')],
c(3,6),res,"*")
c(3,6),res,"*")
#print(fbar(stf[,FcY,'A']))
#stf@harvest[,FcY,c('A','C','D')] <- sweep(stf@harvest[,FcY,c('A','C','D')],
# c(3,6),res,"*")
# update catch and landings for each fleet in Forecast year
for(i in dms$unit){
stf@catch.n[,FcY,i] <- stf@stock.n[,FcY,i]*(1-exp(-unitSums(stf@harvest[,FcY])-stf@m[,FcY,i]))*(stf@harvest[,FcY,i]/(unitSums(stf@harvest[,FcY])+stf@m[,FcY,i]))
stf@catch[,FcY,i] <- computeCatch(stf[,FcY,i])
stf@landings.n[,FcY,i] <- stf@catch.n[,FcY,i]
stf@landings[,FcY,i] <- computeLandings(stf[,FcY,i])
}
# apply transfer
CATCH[,FcY,,,'A'] <- stf@catch[,FcY,'A']# + TAC_var[FcY,,'Ctransfer']*TACS[,FcY,,,'C']
CATCH[,FcY,,,'B'] <- stf@catch[,FcY,'B']
CATCH[,FcY,,,'C'] <- stf@catch[,FcY,'C']#*(1-TAC_var[FcY,,'Ctransfer'])
CATCH[,FcY,,,'D'] <- stf@catch[,FcY,'D']
stf@harvest[,FcY] <- fleet.harvest(stk=stf,
iYr=FcY,
CATCH=CATCH[,FcY])
#print(stf@harvest[,ImY])
#print(stf@stock.n[,ImY,,,,1,drop=T])
#print(ssb(stf[,ImY]))
# update catch and landings for each fleet in Forecast year
for(i in dms$unit){
......
projectNSH_FMSYAR_no_transfer <- function(iStocks,
iFishery,
cYr,
iTAC,
iTAC_var,
iHistMaxYr,
mpPoints){
stk <- iStocks
#===============================================================================
# Setup control file
#===============================================================================
DtY <- ac(range(stk)["maxyear"]-1) #Data year
ImY <- ac(an(DtY)+1) #Intermediate year
FcY <- ac(an(DtY)+2) #Forecast year
CtY <- ac(an(DtY)+3) #Continuation year
CtY1 <- ac(an(DtY)+4)
FuY <- c(ImY,FcY,CtY,CtY1)#Future years
pyears <- an(DtY) - iHistMaxYr #Years away from original historic max year
#- We substract pyears from the starting point, because with every year forward, we move the year ranges one forward too.
#RECS <- FLQuants("ImY"=stk@stock.n[1,ac(ImY)],"FcY"=exp(apply(log(rec(stk)[,ac((range(stk)["maxyear"]-8-pyears):(range(stk)["maxyear"]))]),3:6,mean,na.rm=T)),
# "CtY"=exp(apply(log(rec(stk)[,ac((range(stk)["maxyear"]-8-pyears):(range(stk)["maxyear"]))]),3:6,mean,na.rm=T)))
RECS <- FLQuants( "ImY" =stk@stock.n[1,ac(ImY)],
"FcY" =exp(apply(log(rec(stk)[,ac((an(DtY)-9):DtY)]),3:6,mean,na.rm=T)),
"CtY" =exp(apply(log(rec(stk)[,ac((an(DtY)-9):DtY)]),3:6,mean,na.rm=T)))
yrs1 <- list("m.spwn","harvest.spwn","stock.wt")
yrs3 <- list("mat")
yrs5 <- list("m")
dsc <- "North Sea Herring"
nam <- "NSH"
dms <- dimnames(stk@m)
dms$year <- c(rev(rev(dms$year)[1:3]),FcY,CtY,CtY1)
dms$unit <- c("A","B","C","D")
f01 <- ac(0:1)
f26 <- ac(2:6)
#===============================================================================
# Setup stock file
#===============================================================================
stf <- FLStock(FLQuant(NA,dimnames=dms),name=nam,desc=dsc)
for(i in dms$unit) stf[,,i] <- window(stk,start=an(dms$year)[1],end=rev(an(dms$year))[1])
units(stf) <- units(stk)
# Fill slots that are the same for all fleets
for(i in c(unlist(yrs1),unlist(yrs3),unlist(yrs5))){
if(i %in% unlist(yrs1)){ for(j in dms$unit){ slot(stf,i)[,FuY,j] <- slot(stk,i)[,DtY]}}
if(i %in% unlist(yrs3)){ for(j in dms$unit){ slot(stf,i)[,FuY,j] <- yearMeans(slot(stk,i)[,ac((an(DtY)-2):an(DtY))])}}
if(i %in% unlist(yrs5)){ for(j in dms$unit){ slot(stf,i)[,FuY,j] <- yearMeans(slot(stk,i)[,ac((an(DtY)-4):an(DtY))])}}
}
# Fill slots that are unique for the fleets
for(iFuY in FuY){
stf@harvest[,iFuY] <- sweep(sweep(iFishery@landings.n[,DtY],c(1,6),areaSums(iFishery@landings.n[,DtY]),"/"),c(1,6),stk@harvest[,DtY],"*")
stf@catch.wt[,iFuY] <- iFishery@landings.wt[,ac(DtY)]
stf@landings.wt[,iFuY] <- iFishery@landings.wt[,ac(DtY)]
}
# Fill slots that have no meaning for NSAS
stf@discards.n[] <- 0
stf@discards[] <- 0
stf@discards.wt[] <- 0
#===============================================================================
# Intermediate year
#===============================================================================
iCATCH <- iTAC
iCATCH[,ImY,,,"A"] <- iTAC[,ImY,,,"A"]# + iTAC_var[ImY,,'Ctransfer'] * iTAC[,ImY,,,"C"]
iCATCH[,ImY,,,"B"] <- iTAC[,ImY,,,"B",drop=T]*iTAC_var[ImY,,'Buptake',drop=T]
iCATCH[,ImY,,,"C"] <- (1-iTAC_var[ImY,,'Ctransfer'])#*iTAC[,ImY,,,"C"]*iTAC_var[ImY,,'Csplit']
iCATCH[,ImY,,,"D"] <- iTAC[,ImY,,,"D",drop=T]*iTAC_var[ImY,,'Duptake',drop=T]*iTAC_var[ImY,,'Dsplit',drop=T]
# computing F
stf@harvest[,ImY] <- fleet.harvest( stk=stf,
iYr=ImY,
CATCH=iCATCH[,ImY])
for(i in dms$unit){
stf@catch.n[,ImY,i] <- stf@stock.n[,ImY,i]*(1-exp(-unitSums(stf@harvest[,ImY])-stf@m[,ImY,i]))*(stf@harvest[,ImY,i]/(unitSums(stf@harvest[,ImY])+stf@m[,ImY,i]))
stf@catch[,ImY,i] <- computeCatch(stf[,ImY,i])
stf@landings.n[,ImY,i] <- stf@catch.n[,ImY,i]
stf@landings[,ImY,i] <- computeLandings(stf[,ImY,i])
}
#print(fbar(stf[,,'A']))
#print(iCATCH[,ImY,,,"A"])
#===============================================================================
# Forecast year
#===============================================================================
for(i in dms$unit) stf@stock.n[1,FcY,i] <- RECS$FcY
for(i in dms$unit) stf@stock.n[2:(dims(stf)$age-1),FcY,i] <- (stf@stock.n[,ImY,1]*exp(-unitSums(stf@harvest[,ImY])-stf@m[,ImY,1]))[ac(range(stf)["min"]:(range(stf)["max"]-2)),]
for(i in dms$unit) stf@stock.n[dims(stf)$age,FcY,i] <- apply((stf@stock.n[,ImY,1]*1/exp(unitSums(stf@harvest[,ImY])+stf@m[,ImY,1]))[ac((range(stf)["max"]-1):range(stf)["max"]),],2:6,sum,na.rm=T)
stf@harvest[,FcY] <- stf@harvest[,ImY]
#print(fbar(stf[,,'A']))
###--- Management options ---###