diff --git a/SS_prelim.tpl b/SS_prelim.tpl index c37b8fc0..736a698c 100644 --- a/SS_prelim.tpl +++ b/SS_prelim.tpl @@ -859,6 +859,11 @@ echoinput << " MG_parms after check " << MGparm << endl; MGparm_use = value(MGparm); + if (SR_fxn == 3) + { + if (SRparm_PH(2) < 0 && SRparm(2) >= 0.999999) use_steepness = 0; + } + echoinput << endl << " now check SRparm bounds and priors and do jitter if requested " << endl; for (i = 1; i <= N_SRparm3; i++) diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl index a8f871a3..d7f3f463 100644 --- a/SS_readcontrol_330.tpl +++ b/SS_readcontrol_330.tpl @@ -2191,6 +2191,7 @@ ivector SRparm_PRtype(1,N_SRparm3) vector SRparm_CV(1,N_SRparm3) ivector SRparm_PH(1,N_SRparm3) + int use_steepness // 0/1 to do full steepness calculations with Beverton-Holt; set to 0 if h = 1 and not estimated LOCAL_CALCS // clang-format on @@ -5956,6 +5957,14 @@ active_parm(active_count) = ParCount; } } + + // flag to bypass steepness calculations + use_steepness = 1; // default + // following test will occur in preliminary calcs such that it will work even if values are read from ss3.par + // if (SR_fxn == 3) + // { + // if (SRparm_PH(2) < 0 || SRparm(2) >= 0.999999) use_steepness = 0; + // } if (recdev_cycle > 0) { diff --git a/SS_recruit.tpl b/SS_recruit.tpl index b718fed4..9a1ab894 100644 --- a/SS_recruit.tpl +++ b/SS_recruit.tpl @@ -50,8 +50,15 @@ FUNCTION dvariable Spawn_Recr(const dvar_vector& SRparm_work, const prevariable& case 3: // Beverton-Holt { steepness = SRparm_work(2); - NewRecruits = (4. * steepness * Recr_virgin_use * SSB_curr_adj) / + if (use_steepness == 1) + { + NewRecruits = (4. * steepness * Recr_virgin_use * SSB_curr_adj) / (SSB_virgin_use * (1. - steepness) + (5. * steepness - 1.) * SSB_curr_adj); + } + else + { + NewRecruits = Recr_virgin_use; + } break; } @@ -318,16 +325,22 @@ FUNCTION dvar_vector Equil_Spawn_Recr_Fxn(const dvar_vector& SRparm, // SS3 previously used alternative formulation: R = A*S/(B+S) // converting SS3 to align with WHAM + if (use_steepness == 1) + { alpha = 4.0 * steepness / (SSBpR_virgin_use * (1. - steepness)); beta = (5.0 * steepness - 1.0) / ((1 - steepness) * SSB_virgin_use); - // " h " << steepness << " derive " << alpha * SSBpR_virgin / (4. + alpha * SSBpR_virgin) << " " << endl; - // " R0 " << Recr_virgin_use << " derive " << 1. / beta * (alpha - 1./SSBpR_virgin) << endl; // report5 <<" SSB_unf "<