Skip to content

Commit 6ebb074

Browse files
Rick-Methot-NOAARick-Methot-NOAA
committed
assure align of SRR curve and SPR/YPR profile
1 parent fcd42a0 commit 6ebb074

File tree

5 files changed

+38
-37
lines changed

5 files changed

+38
-37
lines changed

‎SS_benchfore.tpl‎

Lines changed: 22 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
// SPR refers to spawner potential ratio which is the ratio of SSBpR at some level of F to SSBpR with F = 0
1111

1212
// SSBpR_virgin is calculated in popdyn using the start year biology
13-
// SSBpR_virgin_adj used in equil_spawn_recr. _adj means could be updated if SR_update_SSBpR0_timeseries == 1.
13+
// SSBpR_virgin_4SRR is passed to equil_spawn_recr. _adj means could be updated if SR_update_SSBpR0_timeseries == 1.
1414
// Also used to get alpha in equil_spawn_recr B-H
1515

1616
FUNCTION void setup_Benchmark() // and forecast
@@ -749,12 +749,19 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
749749
SR_parm_work(j) = temp / (Bmark_Yr(10) - Bmark_Yr(9) + 1.);
750750
}
751751
}
752+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin_4SRR, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
753+
report5 << endl << " virgin - SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl<<endl;
752754
if (SR_update_SSBpR0_bmark == 0) // use virgin biology for the spawner-recruitment R0,h calculations in bmark
753755
{
754-
Recr_unf = Recr_virgin;
755-
SSB_unf = SSB_virgin;
756-
SSBpR_unf = SSB_unf / Recr_unf;
757-
SSBpR_virgin_adj = SSB_unf / Recr_unf; // so same as SSBpR_virgin, but is the version that will be used in equil_spawn_recr()
756+
Recr_unf = Recr_virgin;
757+
SSB_unf = SSB_virgin;
758+
Fishon = 0;
759+
Recr_unf = Recr_virgin;
760+
Do_Equil_Calc(Recr_unf); // this calcs SSB using benchmark biology
761+
temp1 = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
762+
temp2 = temp1 / Recr_unf;
763+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin_4SRR, temp2); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
764+
report5 << " virgSRR_&_bmarkbio - SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " " << SSBpR_virgin_4SRR << " equil: " << Equ_SpawnRecr_Result << endl<<endl;
758765
}
759766
else // use benchmark biology in the spawner-recruitment R0,h calculations
760767
{
@@ -763,23 +770,15 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
763770
Do_Equil_Calc(Recr_unf); // this calcs SSB using benchmark biology
764771
SSB_unf = SSB_equil; // equilibrium unfished SSB using the benchmark averaged Recr_unf and benchmark averaged biology
765772
SSBpR_unf = SSB_unf / Recr_unf; // this corresponds to the biology for benchmark average years, not the virgin SSB_virgin
766-
SSBpR_virgin_adj = SSB_unf / Recr_unf; // update
773+
SSBpR_virgin_4SRR = SSB_unf / Recr_unf; // update
774+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_unf); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
775+
report5 << " bmarkSRR_&_bmarkbio - SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_unf << " " << SSBpR_virgin_4SRR << " equil: " << Equ_SpawnRecr_Result << endl<<endl;
767776
}
768777
if (show_MSY == 1)
769778
{
770779
report5 << "SR_parms for benchmark: " << SR_parm_work << endl
771780
<< "Benchmark biology averaged over years: " << Bmark_Yr(1) << " " << Bmark_Yr(2) << endl << endl <<
772781
"input_SR_update_SSBpR0_rd: " << SR_update_SSBpR0_rd << "; flag for updating SSBpR0_Bmark: " << SR_update_SSBpR0_bmark << endl;
773-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
774-
report5 << " Virgin SSB, R0, SPR0: " << SSB_virgin << " " << Recr_virgin << " " << SSBpR_virgin << " " << SSB_virgin/Recr_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
775-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin, SSBpR_virgin); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
776-
report5 << " Benchmark SSB, R0, virgin_SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_virgin << " equil: " << Equ_SpawnRecr_Result << endl;
777-
if ( SR_update_SSBpR0_bmark == 1)
778-
{
779-
report5 << "SPR0 for equilibrium spawner-recruit based on benchmark biology, not virgin biology" << endl;
780-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_virgin_adj); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
781-
report5 << " Benchmark SSB, R0, bench_SPR0: " << SSB_unf << " " << Recr_unf << " " << SSBpR_virgin_adj << " equil: " << Equ_SpawnRecr_Result << endl;
782-
}
783782
Mgmt_quant(19) = Recr_unf;
784783
Mgmt_quant(20) = SSB_unf;
785784
Mgmt_quant(21) = SSB_unf; // placeholder to be replaced by SSB_HCR_infl
@@ -934,7 +933,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
934933
935934
// SPAWN-RECR: calc equil spawn-recr in YPR; need to make this area-specific
936935
SSBpR_temp = SSB_equil; // based on most recent call to Do_Equil_Calc
937-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
936+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
938937
939938
Bspr = Equ_SpawnRecr_Result(1);
940939
Bspr_rec = Equ_SpawnRecr_Result(2);
@@ -1049,7 +1048,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
10491048
if (rundetail > 0 && mceval_counter == 0 && show_MSY == 1)
10501049
echoinput << "Calculated F0.1: " << Btgt_Fmult << endl;
10511050
SSBpR_temp = SSB_equil;
1052-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
1051+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
10531052
Btgt = Equ_SpawnRecr_Result(1);
10541053
Btgt_Rec = Equ_SpawnRecr_Result(2);
10551054
YPR_Btgt_enc = YPR_enc; // total encountered yield per recruit
@@ -1145,7 +1144,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
11451144
SSBpR_temp = SSB_equil;
11461145
SPR_Btgt = SSBpR_temp / SSBpR_unf; // where SSBpR_unf = SSB_unf / Recr_unf so units of SSB/R; so result is SPR_Btgt = (fished SSB/R) / (unfished SSB/R)
11471146
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
1148-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
1147+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
11491148
yld1(ii) = Equ_SpawnRecr_Result(1);
11501149
}
11511150
@@ -1185,7 +1184,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
11851184
report5 << " " << SSB_equil_pop_gp(p, gp) * Equ_SpawnRecr_Result(2);
11861185
}
11871186
report5 << endl;
1188-
// " SSB_unf " << SSB_unf << " recr " <<Recr_unf<< " SSB_equil " << SSB_equil << " ssbpr_virginadj: " << SSBpR_virgin_adj << " Btgt " << Btgt << endl;
1187+
// " SSB_unf " << SSB_unf << " recr " <<Recr_unf<< " SSB_equil " << SSB_equil << " ssbpr_virginadj: " << SSBpR_virgin_4SRR << " Btgt " << Btgt << endl;
11891188
}
11901189
} // end search loop
11911190
@@ -1366,7 +1365,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
13661365
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
13671366
MSY_SPR = SSB_equil / SSBpR_unf;
13681367
SSBpR_temp = SSB_equil;
1369-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
1368+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
13701369
Bmsy = Equ_SpawnRecr_Result(1); // with MSY set to SPR, not directly estimated
13711370
Recr_msy = Equ_SpawnRecr_Result(2);
13721371
yld1(1) = YPR_opt * Recr_msy;
@@ -1460,7 +1459,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
14601459
// SPAWN-RECR: calc spawn-recr for MSY calcs; need to make area-specific
14611460
MSY_SPR = SSB_equil / SSBpR_unf;
14621461
SSBpR_temp = SSB_equil;
1463-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
1462+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
14641463
Bmsy = Equ_SpawnRecr_Result(1); // MSY is directly estimated
14651464
Recr_msy = Equ_SpawnRecr_Result(2);
14661465
Profit = (PricePerF * YPR_val_vec) * Recr_msy - Cost;
@@ -1756,7 +1755,7 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
17561755
SPR_Btgt2 = SSB_equil / SSBpR_unf;
17571756
// SPAWN-RECR: calc equil spawn-recr for Btarget calcs; need to make area-specific
17581757
SSBpR_temp = SSB_equil;
1759-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
1758+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
17601759
yld1(ii) = Equ_SpawnRecr_Result(1);
17611760
}
17621761

‎SS_param.tpl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -166,7 +166,7 @@ PARAMETER_SECTION
166166
number half_sigmaRsq;
167167
number sigmaR;
168168
number SSBpR_virgin;
169-
number SSBpR_virgin_adj;
169+
number SSBpR_virgin_4SRR;
170170
number regime_change;
171171
number rho;
172172
number dirichlet_Parm;

‎SS_popdyn.tpl‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -344,7 +344,7 @@ FUNCTION void get_initial_conditions()
344344
equ_Recr = 1.0;
345345
Do_Equil_Calc(equ_Recr); // call function to do equilibrium calculation. Returns SPR because R = 1.0
346346
SSBpR_virgin = SSB_equil; // spawners per recruit. Needed for Sr_fxn = 10
347-
SSBpR_virgin_adj = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
347+
SSBpR_virgin_4SRR = SSB_equil; // also needed for Sr_fxn 10. Will get revised in benchmark to use averaged biology if requested.
348348
if(SR_fxn == 10) // B-H with a,b
349349
{
350350
// WHAM based on R = A*S/(1+B*S)
@@ -355,8 +355,8 @@ FUNCTION void get_initial_conditions()
355355
356356
alpha = mfexp(SR_parm(3));
357357
beta = mfexp(SR_parm(4));
358-
steepness = alpha * SSBpR_virgin_adj / (4. + alpha * SSBpR_virgin_adj);
359-
Recr_virgin = 1. / beta * (alpha - (1. / SSBpR_virgin_adj));
358+
steepness = alpha * SSBpR_virgin_4SRR / (4. + alpha * SSBpR_virgin_4SRR);
359+
Recr_virgin = 1. / beta * (alpha - (1. / SSBpR_virgin_4SRR));
360360
// warning << " before AB_calcs " << "parm " << SR_parm(1) << " calc " << log(Recr_virgin) << endl;
361361
SR_parm(1) = log(Recr_virgin);
362362
SR_parm(2) = steepness;
@@ -536,7 +536,7 @@ FUNCTION void get_initial_conditions()
536536
SSBpR_temp = SSB_equil / equ_Recr; // spawners per recruit at initial F
537537
// get equilibrium SSB and recruitment from SSBpR_temp, Recr_virgin and virgin steepness
538538
// this is the initial year, so no time-vary effects available, so uses _virgin quantities for spawner-recruitment
539-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
539+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_virgin, Recr_virgin, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
540540
R1_exp = Equ_SpawnRecr_Result(2); // set the expected recruitment equal to this equilibrium
541541
exp_rec(eq_yr, 1) = R1_exp;
542542
@@ -1040,7 +1040,7 @@ FUNCTION void get_time_series()
10401040
bio_yr = y;
10411041
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation with current year's biology and adjusted R0
10421042
SSB_use = SSB_equil;
1043-
SSBpR_virgin_adj = SSB_use / R0_use; // update
1043+
SSBpR_virgin_4SRR = SSB_use / R0_use; // update
10441044
if (fishery_on_off == 1)
10451045
{
10461046
Fishon = 1;
@@ -1500,7 +1500,7 @@ FUNCTION void get_time_series()
15001500
bio_yr = y;
15011501
Do_Equil_Calc(R0_use); // call function to do equilibrium calculation
15021502
SSB_use = SSB_equil;
1503-
SSBpR_virgin_adj = SSB_use / R0_use; // update
1503+
SSBpR_virgin_4SRR = SSB_use / R0_use; // update
15041504
if (fishery_on_off == 1)
15051505
{
15061506
Fishon = 1;

‎SS_recruit.tpl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -268,7 +268,7 @@ FUNCTION void apply_recdev(prevariable& NewRecruits, const prevariable& Recr_vir
268268
/* SS_Label_FUNCTION 44 Equil_Spawn_Recr_Fxn */
269269
// SPAWN-RECR: function Equil_Spawn_Recr_Fxn
270270
FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
271-
const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSBpR_virgin_adj, const prevariable& SSBpR_current)
271+
const prevariable& SSB_virgin_use, const prevariable& Recr_virgin_use, const prevariable& SSBpR_virgin_use, const prevariable& SSBpR_current)
272272
{
273273
RETURN_ARRAYS_INCREMENT();
274274
dvar_vector Equil_Spawn_Recr_Calc(1, 2); // values to return 1 is B_equil, 2 is R_equil
@@ -316,7 +316,7 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm,
316316
317317
// SS3 previously used alternative formulation: R = A*S/(B+S)
318318
// converting SS3 to align with WHAM
319-
alpha = 4.0 * steepness / (SSBpR_virgin_adj * (1. - steepness));
319+
alpha = 4.0 * steepness / (SSBpR_virgin_use * (1. - steepness));
320320
beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin_use);
321321
// " h " << steepness << " derive " << alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin) << " " << endl;
322322
// " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SSBpR_virgin) << endl;

‎SS_write_report.tpl‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1858,7 +1858,7 @@ FUNCTION void write_bigoutput()
18581858
{
18591859
case 3: // Beverton-Holt with steepness
18601860
{
1861-
alpha = 4.0 * steepness / (SSBpR_virgin_adj * (1. - steepness));
1861+
alpha = 4.0 * steepness / (SSBpR_virgin_4SRR * (1. - steepness));
18621862
beta = (5.0 * steepness - 1.0) / ((1. - steepness) * SSB_virgin);
18631863
SS2out << "Ln(R0): " << SR_parm(1) << endl << "R0: " << mfexp(SR_parm(1)) << endl;
18641864
SS2out << "steepness: " << steepness << endl;
@@ -1870,8 +1870,8 @@ FUNCTION void write_bigoutput()
18701870
{
18711871
SS2out << "Ln(alpha): " << SR_parm(3) << " alpha " << mfexp(SR_parm(3)) << endl;
18721872
SS2out << "Ln(beta): " << SR_parm(4) << " beta " << mfexp(SR_parm(4)) << endl;
1873-
SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SSBpR_virgin_adj))) << endl; // virgin R0
1874-
SS2out << "steepness_derived: " << alpha * SSBpR_virgin_adj / (4. + alpha * SSBpR_virgin_adj) << endl; // steepness virgin
1873+
SS2out << "ln(R0)_derived: " << log( 1. / beta * (alpha - (1. / SSBpR_virgin_4SRR))) << endl; // virgin R0
1874+
SS2out << "steepness_derived: " << alpha * SSBpR_virgin_4SRR / (4. + alpha * SSBpR_virgin_4SRR) << endl; // steepness virgin
18751875
break;
18761876
}
18771877
case 8:
@@ -4741,10 +4741,12 @@ FUNCTION void SPR_profile()
47414741
Wt_Age_beg(s) = Wt_Age_t(t, 0); // used for smrybio
47424742
Wt_Age_mid(s) = Wt_Age_t(t, -1);
47434743
if (s == spawn_seas)
4744+
{
47444745
fec = Wt_Age_t(t, -2);
4745-
report5 << " repro_output for spr/ypr: " << fec(1) << endl;
4746+
SS2out << " repro_output for spr/ypr: " << fec(1) << endl;}
47464747
}
47474748

4749+
SS2out << "unfished values for SRR: SSB " << SSB_unf << " R " << Recr_unf << " SSBpR " << SSBpR_virgin_4SRR << endl;
47484750
SS2out << "SPRloop Iter Bycatch Fmult F_std SSBpR YpR_dead YpR_dead*Recr YpR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
47494751
for (f = 1; f <= Nfleet; f++)
47504752
{
@@ -4902,7 +4904,7 @@ FUNCTION void SPR_profile()
49024904
Do_Equil_Calc(equ_Recr);
49034905
// SPAWN-RECR: calc equil spawn-recr in the SPR loop
49044906
SSBpR_temp = SSB_equil;
4905-
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_adj, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
4907+
Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_virgin_4SRR, SSBpR_temp); // returns 2 element vector containing equilibrium biomass and recruitment at this SPR
49064908
Btgt_prof = Equ_SpawnRecr_Result(1);
49074909
Btgt_prof_rec = Equ_SpawnRecr_Result(2);
49084910
if (Btgt_prof < 0.001 || Btgt_prof_rec < 0.001)

0 commit comments

Comments
 (0)