Skip to content

Commit

Permalink
Merge pull request #284 from nmfs-stock-synthesis/fix_seasonal_recrui…
Browse files Browse the repository at this point in the history 
…tment

Fix seasonal recruitment
  • Loading branch information
@Rick-Methot-NOAA
Rick-Methot-NOAA authored Apr 13, 2022
2 parents 9f31214 + 73f0109 commit 1486a3d
Show file tree
Hide file tree
Showing 5 changed files with 62 additions and 35 deletions.
54 changes: 36 additions & 18 deletions SS_benchfore.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -356,15 +356,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
}

if(show_MSY==1)
{
report5<<version_info<<endl<<ctime(&start);
report5<<"Bmark_relF(by_fleet_&seas) (excluding non-scaled bycatch fleets)"<<endl<<Bmark_RelF_Use<<endl<<"#"<<endl;
report5<<"Bmark_histF(by_fleet_&seas)"<<endl<<Bmark_HistF<<endl<<"#"<<endl;
report5<<"Bycatch_F"<<endl<<trans(bycatch_F)<<endl<<"#"<<endl;
report5<<"YPR_mask for including catch: "<<endl<<YPR_mask<<endl;
}

if(show_MSY==2)
{
// do not recalc the age-specific vectors
Expand Down Expand Up @@ -397,8 +388,29 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
if(use_morph(g)>0)
{
ALK_idx=(s-1)*N_subseas+mid_subseas; // for midseason
Make_FishSelex();
Make_FishSelex(); // including sel_dead_num
}
if(show_MSY==1)
{
report5<<version_info<<endl<<ctime(&start);
report5<<"Bmark_relF(by_fleet_&seas) (excluding non-scaled bycatch fleets)"<<endl<<Bmark_RelF_Use<<endl<<"#"<<endl;
report5<<"Bmark_histF(by_fleet_&seas)"<<endl<<Bmark_HistF<<endl<<"#"<<endl;
report5<<"Bycatch_F"<<endl<<trans(bycatch_F)<<endl<<"#"<<endl;
report5<<"YPR_mask for including catch: "<<endl<<YPR_mask<<endl;
report5<<"Fecundity: "<<fec(1)<<endl;
for (f=1;f<=Nfleet;f++){
if(fleet_type(f)<=2){
for(s=1;s<=nseas;s++) report5<<f<<" "<<s<<" sel_dead_bio: "<<sel_dead_bio(s,f,1)<<endl;
}}
for (f=1;f<=Nfleet;f++){
if(fleet_type(f)<=2){
for(s=1;s<=nseas;s++) report5<<f<<" "<<s<<" sel_num: "<<sel_num(s,f,1)<<endl;
}}
for (f=1;f<=Nfleet;f++){
if(fleet_type(f)<=2){
for(s=1;s<=nseas;s++) report5<<f<<" "<<s<<" sel_dead_num: "<<sel_dead_num(s,f,1)<<endl;
}}
}
}

maxpossF.initialize();
Expand Down Expand Up @@ -970,13 +982,14 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
{
report5<<j<<" "<<Fmult<<" "<<equ_F_std<<" "<<MSY_SPR<<" "<<yld1(1)<<" "<<Bmsy<<" "<<Recr_msy<<" "<<Bmsy/SSB_unf<<" "
<<dyld <<" "<<dyldp<<" "<<value(sum(equ_catch_fleet(3))*Recr_msy);
report5 << value(equ_catch_fleet(3)*Recr_msy) << " " << Cost << " " << PricePerF*YPR_val_vec*Recr_msy << " " << Profit << " ";
// report5 << " " << value(equ_catch_fleet(3)*Recr_msy) << " " << Cost << " " << PricePerF*YPR_val_vec*Recr_msy << " " << Profit << " ";
report5 << " " << Cost << " " << PricePerF*YPR_val_vec*Recr_msy << " " << Profit << " ";
for (p=1;p<=pop;p++)
for (gp=1;gp<=N_GP;gp++)
{report5<<" "<<SSB_equil_pop_gp(p,gp)*Recr_msy;}
for(int ff=1;ff<=N_catchfleets(0);ff++)
{f=fish_fleet_area(0,ff);
report5<<Hrate(f,bio_t_base+1)<<" ";}
report5<<" "<<Hrate(f,bio_t_base+1);}
report5<<endl;
}
if(j<=9)
Expand Down Expand Up @@ -1045,11 +1058,12 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
}
}
report5 <<gp<<" "<<gg<<" N "<<tempvec_a<<endl;
report5 <<gp<<" "<<gg<<" Z ";
for (a=0;a<=nages-2;a++)
{report5<<-log(tempvec_a(a+1)/tempvec_a(a))<<" ";}
report5<<" NA NA"<<endl;

}
report5 <<"natM"<<endl<<natM<<endl;
report5 <<"sel_dead_num"<<endl<<sel_dead_num<<endl;
report5 <<"equil_Z"<<endl<<equ_Z<<endl;
report5 <<"equil_catage"<<endl<<equ_catage<<endl;

Fishon=0;
Do_Equil_Calc(equ_Recr);
Expand Down Expand Up @@ -1867,7 +1881,9 @@ FUNCTION void Get_Forecast()
// if(y==endyr+1) natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle))=0.0; // to negate the additive code
natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle)) = Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g))*
mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle));
if(Fcast_Loop1==jloop && ABC_Loop==ABC_Loop_end) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
if(Fcast_Loop1==jloop && ABC_Loop==ABC_Loop_end){
if(Settle_seas(settle)==s) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
}
// the adjustment for mortality increases recruit value for elapsed time since begin of season because M will then be applied from beginning of season
}
}
Expand Down Expand Up @@ -2388,7 +2404,9 @@ FUNCTION void Get_Forecast()
// natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle)) += Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g))*
natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle)) = Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g))*
mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle));
if(Fcast_Loop1==jloop && ABC_Loop==ABC_Loop_end) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
if(Fcast_Loop1==jloop && ABC_Loop==ABC_Loop_end){
if(Settle_seas(settle)==s) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
}
}
}
}
Expand Down
6 changes: 6 additions & 0 deletions SS_biofxn.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1337,6 +1337,11 @@ FUNCTION void get_mat_fec();
GPat=GP4(g);
gg=sx(g);
gp=GPat; //
if(WTage_rd==1){
fec(g)=Wt_Age_t(t,-2,g);
}
else{ // make fecundity from biology

if(do_once==1) echoinput<<"fecundity option: "<<Fecund_Option<<" parms: "<<wtlen_p(GPat)(5,6)<<endl;

switch (Fecund_Option)
Expand Down Expand Up @@ -1470,6 +1475,7 @@ FUNCTION void get_mat_fec();
" fecundity_age: "<<fec(g)<<endl;
}
#endif
}
} // end g loop
// end maturity and fecundity in spawn_seas
}
Expand Down
4 changes: 4 additions & 0 deletions SS_global.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -782,6 +782,7 @@ FINAL_SECTION
Get_Forecast();
if(mceval_phase()==0) cout<<" finished forecast"<<endl;
}

if(write_bodywt>0)
{
bodywtout<<-9999<<" "<<1<<" "<<1<<" "<<1<<" "<<1<<" "<<0<<" "<<Wt_Age_mid(1,1)<<" #terminator "<<endl;
Expand All @@ -794,6 +795,7 @@ FINAL_SECTION
if(mceval_phase()==0) cout<<" finished StdDev quantities"<<endl;
get_posteriors();
if(mceval_phase()==0) cout<<" finished posteriors"<<endl;

// SS_Label_Info_12.4.2 #Call fxn write_summaryoutput()
if(Do_CumReport>0) write_summaryoutput();

Expand Down Expand Up @@ -840,6 +842,8 @@ FINAL_SECTION

if (pick_report_use(54)=="Y" && Do_Benchmark>0)
{
cout<<"setup_benchmark: "<<endl;
setup_Benchmark();
cout<<"SPR_profile: ";
SPR_profile();
cout<<" finished "<<endl;
Expand Down
8 changes: 4 additions & 4 deletions SS_popdyn.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ FUNCTION void get_initial_conditions()
}
Save_PopAge(t+s,p,g)=natage(t+s,p,g);
Save_PopAge(t+s,p+pop,g)=elem_prod(natage(t+s,p,g),mfexp(-Z_rate(t+s,p,g)*0.5*seasdur(s)));
Recr(p,t+1+Settle_seas_offset(settle_g(g)))+=equ_Recr*recr_dist(y,GP(g),settle_g(g),p)*platoon_distr(GP2(g));
if(Settle_seas(settle_g(g))==s) Recr(p,t+1+Settle_seas_offset(settle_g(g)))+=equ_Recr*recr_dist(y,GP(g),settle_g(g),p)*platoon_distr(GP2(g));
Save_PopBio(t+s,p,g)=elem_prod(natage(t+s,p,g),Wt_Age_beg(s,g));
Save_PopBio(t+s,p+pop,g)=elem_prod(Save_PopAge(t+s,p+pop,g),Wt_Age_mid(s,g));
}
Expand Down Expand Up @@ -516,7 +516,7 @@ FUNCTION void get_initial_conditions()
Save_PopAge(t+s,p+pop,g)=elem_prod(natage(t+s,p,g),mfexp(-Z_rate(t+s,p,g)*0.5*seasdur(s)));
Save_PopBio(t+s,p,g)=elem_prod(natage(t+s,p,g),Wt_Age_beg(s,g));
Save_PopBio(t+s,p+pop,g)=elem_prod(Save_PopAge(t+s,p+pop,g),Wt_Age_mid(s,g));
Recr(p,t+1+Settle_seas_offset(settle_g(g)))+=equ_Recr*recr_dist(y,GP(g),settle_g(g),p)*platoon_distr(GP2(g));
if(Settle_seas(settle_g(g))==s) Recr(p,t+1+Settle_seas_offset(settle_g(g)))+=equ_Recr*recr_dist(y,GP(g),settle_g(g),p)*platoon_distr(GP2(g));
}
}

Expand Down Expand Up @@ -889,7 +889,7 @@ FUNCTION void get_time_series()
natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle)) +=
Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g))*
mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle));
Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
if(Settle_seas(settle)==s) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
// the adjustment for mortality increases recruit value for elapsed time since begin of season because M will then be applied from beginning of season
if(docheckup==1) echoinput<<y<<" Recruits, dist, surv, result"<<Recruits<<" "<<recr_dist(y,GP(g),settle,p)<<" "<<mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle))<<" "<<natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle))<<endl;
}
Expand Down Expand Up @@ -1282,7 +1282,7 @@ FUNCTION void get_time_series()

natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle)) += Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g))*
mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle));
Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
if(Settle_seas(settle)==s) Recr(p,t+Settle_seas_offset(settle))+=Recruits*recr_dist(y,GP(g),settle,p)*platoon_distr(GP2(g));
if(docheckup==1) echoinput<<y<<" Recruits, dist, surv, result"<<Recruits<<" "<<recr_dist(y,GP(g),settle,p)<<" "<<mfexp(natM(s,GP3(g),Settle_age(settle))*Settle_timing_seas(settle))<<" "<<natage(t+Settle_seas_offset(settle),p,g,Settle_age(settle))<<endl;
}
}
Expand Down
25 changes: 12 additions & 13 deletions SS_write_report.tpl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3359,14 +3359,6 @@ FUNCTION void SPR_profile()
dvariable YPR_last;

SS2out<<endl<<pick_report_name(54)<<endl;
SS2out<<"SPRloop Iter Bycatch Fmult F_report SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Dead";}
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Ret";}
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Age";}
for (p=1;p<=pop;p++)
for (gp=1;gp<=N_GP;gp++)
{SS2out<<" SSB_Area:"<<p<<"_GP:"<<gp;}
SS2out<<endl;
y=styr-3;
yz=y;
bio_yr=y;
Expand All @@ -3376,9 +3368,11 @@ FUNCTION void SPR_profile()

// SPAWN-RECR: call make_fecundity for benchmark bio for SPR loop
// this code section that creates fecundity and selectivity seems antiquated; why is it different from the averages used for benchmark?

for (s=1;s<=nseas;s++)
{
t = styr-3*nseas+s-1;

if(MG_active(2)>0 || MG_active(3)>0 || save_for_report>0 || WTage_rd>0)
{
subseas=1;
Expand All @@ -3390,14 +3384,14 @@ FUNCTION void SPR_profile()
if(s==spawn_seas)
{
subseas=spawn_subseas;
if(spawn_subseas!=1 && spawn_subseas!=mid_subseas)
{
if(spawn_subseas!=1 && spawn_subseas!=mid_subseas){
//don't call get_growth3(subseas) because using an average ave_size
Make_AgeLength_Key(s, subseas); // spawn subseas
}
get_mat_fec();
}
}

for (g=1;g<=gmorph;g++)
if(use_morph(g)>0)
{
Expand All @@ -3406,7 +3400,14 @@ FUNCTION void SPR_profile()
}
}


SS2out<<"SPRloop Iter Bycatch Fmult F_report SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Dead";}
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Ret";}
for (f=1;f<=Nfleet;f++) {if(fleet_type(f)<=2) SS2out<<" "<<fleetname(f)<<"("<<f<<")Age";}
for (p=1;p<=pop;p++)
for (gp=1;gp<=N_GP;gp++)
{SS2out<<" SSB_Area:"<<p<<"_GP:"<<gp;}
SS2out<<endl;
equ_Recr=1.0;
Fishon=0;
int SPRloop1_end;
Expand Down Expand Up @@ -3513,11 +3514,9 @@ FUNCTION void SPR_profile()
Do_Equil_Calc(equ_Recr);
// SPAWN-RECR: calc equil spawn-recr in the SPR loop
SPR_temp=SSB_equil;
// Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm(2), SR_parm(3), SSB_virgin, Recr_virgin, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work(2), SR_parm_work(3), SSB_unf, Recr_unf, SPR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
Btgt_prof=Equ_SpawnRecr_Result(1);
Btgt_prof_rec=Equ_SpawnRecr_Result(2);

if(Btgt_prof<0.001 || Btgt_prof_rec<0.001)
{
Btgt_prof_rec=0.0; Btgt_prof=0.;
Expand Down

0 comments on commit 1486a3d

Please sign in to comment.