diff --git a/SS_biofxn.tpl b/SS_biofxn.tpl
index b9781597..bd4efd23 100644
--- a/SS_biofxn.tpl
+++ b/SS_biofxn.tpl
@@ -1,7 +1,7 @@
// SS_Label_file #9. **SS_biofxn.tpl**
// SS_Label_file # * get_MGsetup() // get parameter values for this year
-// SS_Label_file # * get_growth1() // prep growth quantities
-// SS_Label_file # * get_growth2() // growth to beginning of each season of upcoming year
+// SS_Label_file # * get_growth1() // prep growth quantities that are not time-varying, especially CV_growth
+// SS_Label_file # * get_growth2() // growth to beginning of each season or first season of upcoming year
// SS_Label_file # * get_growth3() // growth to particular time point in a season
// SS_Label_file # * get_natmort()
// SS_Label_file # * get_recr_distribution()
@@ -79,30 +79,30 @@ FUNCTION void get_MGsetup(const int yz)
FUNCTION void get_growth1()
{
// SS_Label_Info_15.1 #create seasonal effects for growth K, and for wt_len parameters
- if (MGparm_doseas > 0)
+ // VBK_seas(0) contains sum of season-specific durations, it will equal 1.0 in annual models and seasonal models
+ // except as modified by seasonal K effects. It will be less than 1.0 in a seasons-as-years model
+
+ if (MGparm_doseas > 0) // there are seasonal effects
{
if (MGparm_seas_effects(10) > 0) // for seasonal K
{
VBK_seas(0) = 0.0;
for (s = 1; s <= nseas; s++)
{
- VBK_seas(s) = mfexp(MGparm(MGparm_seas_effects(10) + s));
- VBK_seas(0) += VBK_seas(s) * seasdur(s);
+ VBK_seas(s) = mfexp(MGparm(MGparm_seas_effects(10) + s)); // so equals 1.0 if the seasonal parameter = 0.0
+ VBK_seas(0) += VBK_seas(s) * seasdur(s); // sum seasonal effects weighted by season duration
}
}
else
{
VBK_seas = sum(seasdur); // set vector to null effect
}
+
for (gp = 1; gp <= N_GP; gp++)
for (j = 1; j <= 8; j++)
{
- #ifdef DO_ONCE
- {
- if (do_once == 1)
- echoinput << j << " wt_len seas " << MGparm_seas_effects(j) << endl;
- }
- #endif
+ if (do_once == 1)
+ echoinput << j << " wt_len seas_effects " << MGparm_seas_effects(j) << endl;
if (MGparm_seas_effects(j) > 0)
{
wtlen_seas(0, gp, j) = 0.0;
@@ -121,13 +121,15 @@ FUNCTION void get_growth1()
}
}
}
- else
+ else // no seasonal effects
{
- VBK_seas = sum(seasdur); // set vector to null effect
+ VBK_seas(1, nseas) = 1.;
+
+ VBK_seas(0) = sum(seasdur); // set vector to null effect // note that seasons_as_years will have value < 1.0
for (s = 1; s <= nseas; s++)
wtlen_seas(s) = 1.0; // set vector to null effect
}
-
+ if (do_once == 1) echoinput << "VBK_seasonal_adjustment: " << VBK_seas(1, nseas) << " total: " << VBK_seas(0) << endl;
// SS_Label_Info_15.2 #create variability of size-at-age factors using direct assignment or offset approaches
gp = 0;
for (gg = 1; gg <= gender; gg++)
@@ -206,19 +208,25 @@ FUNCTION void get_growth2(const int y)
// called at beginning of each year, so y is known
// if y=styr, then does equilibrium size-at-age according to start year growth parameters
// for any year, calculates for each season the size at the beginning of the next season, with growth increment calculated according to that year's parameters
+ // VBK(gp,a) is the growth parameter for growth type gp and age. VBK for age 0 is used for all ages. set to negative of parameter
+ // VBK_seas(1, nseas) is optional seasonal adjustment to VBK
+ // VBK_by_seas is VBK adjusted by VBK_seas
+ // VBK_work is current value of VBK for the selected model condition
+
//Growth Cessation Model code added by Mark Maunder October 2018
//The growth cessation model is described in
//Maunder, M.N., Deriso, R.B., Schaefer, K.M., Fuller, D.W., Aires-da-Silva, A.M., Minte‑Vera, C.V., Campana, S.E. 2018. The growth cessation model: a growth model for species showing a near cessation in growth with application to bigeye tuna (Thunnus obesus). Marine Biology (2018) 165:76.
- //Ian Taylor derived the formula for Linf
+ //Ian Taylor derived the formula for L_inf
int k2;
int add_age;
int ALK_idx2; // beginning of first subseas of next season
+ int L_inf_plus = 0;
dvariable plusgroupsize;
dvariable current_size;
- dvariable VBK_temp;
- dvariable VBK_temp2; // with VBKseas(s) multiplied
+ dvariable VBK_work; // working version of VBK
+ dvariable VBK_by_seas; // with VBKseas(s) multiplied
dvariable LminR;
dvariable LmaxR;
dvariable LinfR;
@@ -240,7 +248,7 @@ FUNCTION void get_growth2(const int y)
#ifdef DO_ONCE
{
if (do_once == 1)
- echoinput << "GROWTH, yr= " << y << endl;
+ echoinput << "GROWTH2 by season in year " << y << endl;
}
#endif
for (gg = 1; gg <= gender; gg++)
@@ -248,32 +256,36 @@ FUNCTION void get_growth2(const int y)
{
gp++;
Ip = MGparm_point(gg, GPat) + N_natMparms;
- switch (Grow_type) // create specific growth parameters from the mgp_adj list of current MGparms
+ // SS_Label_Info_16.2.1 #set Lmin, Lmax, VBK, Richards to this year's values for mgp_adj
+ // SS_Label_Info_16.2.2 #Set up age specific k
+ if (Grow_type != 7)
{
- case 7: // empirical length
+ // get basic parameters used by all grow_type
+ if (MGparm_def > 1 && gp > 1) // do offset approach
{
- break;
+ Lmin(gp) = Lmin(1) * mfexp(mgp_adj(Ip));
+ Lmax_temp(gp) = Lmax_temp(1) * mfexp(mgp_adj(Ip + 1));
+ VBK(gp) = VBK(1) * mfexp(mgp_adj(Ip + 2)); // assigns to all ages for which VBK is defined
+ }
+ else
+ {
+ Lmin(gp) = mgp_adj(Ip);
+ Lmax_temp(gp) = mgp_adj(Ip + 1); // size at A2; could be 999 to indicate L_inf
+ VBK(gp) = -mgp_adj(Ip + 2); // because always used as negative; assigns to all ages for which VBK is defined
}
+ L_inf(gp) = Lmax_temp(gp); // because this is the typical assignment
+ VBK_work = VBK(gp, 0); // will be reset to VBK(gp,nages) if using age-specific K; note that this inherits the negative sign
- default: // process parameters for all other grow_type
+ // get specific growth parameters: Richards/Gompertz, growth cessation, and age-specific VBK
+ switch (Grow_type)
{
- // SS_Label_Info_16.2.1 #set Lmin, Lmax, VBK, Richards to this year's values for mgp_adj
- if (MGparm_def > 1 && gp > 1) // do offset approach
+ case 1: // standard von Bertallanfy no action needed
{
- Lmin(gp) = Lmin(1) * mfexp(mgp_adj(Ip));
- Lmax_temp(gp) = Lmax_temp(1) * mfexp(mgp_adj(Ip + 1));
- VBK(gp) = VBK(1) * mfexp(mgp_adj(Ip + 2)); // assigns to all ages for which VBK is defined
- }
- else
- {
- Lmin(gp) = mgp_adj(Ip);
- Lmax_temp(gp) = mgp_adj(Ip + 1); // size at A2; could be 999 to indicate Linf
- VBK(gp) = -mgp_adj(Ip + 2); // because always used as negative; assigns to all ages for which VBK is defined
+ if (AFIX2 != 999) L_inf_plus = 1;
+ break;
}
- VBK_temp = VBK(gp, 0); // will be reset to VBK(gp,nages) if using age-specific K
- // SS_Label_Info_16.2.2 #Set up age specific k
- if (Grow_type == 3) // age specific k
+ case 3: // age specific k ascending sequence
{
j = 1;
for (a = 1; a <= nages; a++)
@@ -289,12 +301,15 @@ FUNCTION void get_growth2(const int y)
VBK(gp, a) = VBK(gp, a - 1);
}
}
- VBK_temp = VBK(gp, nages);
+ VBK_work = VBK(gp, nages);
+ if (AFIX2 != 999) L_inf_plus = 1;
+ break;
}
- else if (Grow_type == 4) // age specific k reverse order, so age_k_points need to be descending
+
+ case 4: // age specific k reverse order, so age_k_points need to be descending
{
j = 1;
- VBK(gp, nages) = VBK_temp;
+ VBK(gp, nages) = VBK_work; // oldest age has the base parameter value
for (a = nages - 1; a >= 0; a--)
{
if (a == Age_K_points(j))
@@ -308,29 +323,32 @@ FUNCTION void get_growth2(const int y)
VBK(gp, a) = VBK(gp, a + 1);
}
}
+ if (AFIX2 != 999) L_inf_plus = 1;
+ break;
}
- else if (Grow_type == 5) // age specific k replacement, so age_k_points need to be descending
+
+ case 5: // age specific k replacement, age_k_points need to be descending
{
+ // all ages start at base parameter value, which is VBK_work
j = 1;
for (a = nages; a >= 0; a--)
{
if (a == Age_K_points(j))
{
- VBK(gp, a) = mgp_adj(Ip + 2 + j) * VBK_temp;
+ VBK(gp, a) = mgp_adj(Ip + 2 + j) * VBK_work;
if (j < Age_K_count)
j++;
}
else
{
- VBK(gp, a) = VBK_temp;
+ VBK(gp, a) = VBK_work;
}
}
- VBK_temp = VBK(gp, nages);
+ if (AFIX2 != 999) L_inf_plus = 1;
+ break;
}
- // get Linf from Lmax
- // get Richards or growth cessation parameter if appropriate
- if (Grow_type == 2) // Richards
+ case 2: // Richards or Gompertz growth
{
if (MGparm_def > 1 && gp > 1)
{
@@ -344,21 +362,19 @@ FUNCTION void get_growth2(const int y)
inv_Richards = 1.0 / Richards(gp);
if (AFIX2 == 999)
{
- L_inf(gp) = Lmax_temp(gp);
+ // L_inf(gp) = Lmax_temp(gp); // already set
LinfR = pow(L_inf(gp), Richards(gp));
}
else
{
LmaxR = pow(Lmax_temp(gp), Richards(gp));
- LinfR = LminR + (LmaxR - LminR) / (1. - mfexp(VBK_temp * VBK_seas(0) * (AFIX_delta)));
+ LinfR = LminR + (LmaxR - LminR) / (1. - mfexp(VBK_work * VBK_seas(0) * (AFIX_delta)));
L_inf(gp) = pow(LinfR, inv_Richards);
}
- #ifdef DO_ONCE
- if (do_once == 1)
- echoinput << " linf " << L_inf(gp) << " VBK: " << VBK_temp << endl;
- #endif
+ break;
}
- else if (Grow_type == 8)
+
+ case 8:
{
if (MGparm_def > 1 && gp > 1)
{
@@ -369,51 +385,68 @@ FUNCTION void get_growth2(const int y)
Richards(gp) = mgp_adj(Ip + 3);
}
L_inf(gp) = Lmax_temp(gp);
- VBK_temp = -VBK(gp, 0) * VBK_seas(0);
// t50 is the calculated inflection age for the decline in K
- t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / VBK_temp) - 1.0) / Richards(gp);
- }
- else
- {
- if (AFIX2 == 999)
- {
- L_inf(gp) = Lmax_temp(gp);
- }
- else
- {
- L_inf(gp) = Lmin(gp) + (Lmax_temp(gp) - Lmin(gp)) / (1. - mfexp(VBK_temp * VBK_seas(0) * (AFIX_delta)));
- #ifdef DO_ONCE
- if (do_once == 1)
- echoinput << VBK_temp << " " << VBK_seas(0) << " " << VBK_temp * VBK_seas(0) << " " << Lmax_temp(gp) << " " << L_inf(gp) << endl;
- #endif
- }
+ // note use of -VBK_work because VBK_work has negative sign
+ t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / (-VBK_work * VBK_seas(0))) - 1.0) / Richards(gp);
+ break;
}
+ }
- // SS_Label_Info_16.2.3 #Set up Lmin and Lmax in Start Year
- if (y == styr)
- {
- Cohort_Lmin(gp) = Lmin(gp); // sets for all years and ages
- }
- else if (timevary_MG(y, 2) > 0) // using time-vary growth
+ if (L_inf_plus == 1)
+ {
+ L_inf(gp) = Lmin(gp) + (Lmax_temp(gp) - Lmin(gp)) / (1. - mfexp(VBK_work * VBK_seas(0) * (AFIX_delta)));
+ }
+
+ // SS_Label_Info_16.2.3 #Set up Cohort_Lmin in Start Year
+ if (y == styr)
+ {
+ Cohort_Lmin(gp) = Lmin(gp); // sets for all years and ages
+ }
+ else if (timevary_MG(y, 2) > 0) // using time-vary growth
+ {
+ k = min(nages, (YrMax - y));
+ for (a = 0; a <= k; a++)
{
- k = min(nages, (YrMax - y));
- for (a = 0; a <= k; a++)
- {
- Cohort_Lmin(gp, y + a, a) = Lmin(gp);
- } // sets for future years so cohort remembers its size at birth; with Lmin(gp) being size at birth this year
- }
- } // end setup of parametric growth parameters
- } // end switch between parametric and non-parametric growth
+ Cohort_Lmin(gp, y + a, a) = Lmin(gp);
+ } // sets for future years so cohort remembers its size at birth; with Lmin(gp) being size at birth this year
+ }
+ } // end parametric growth parameters
+ else
+ {
+ // grow_type == 7} // not implemented placeholder for empirical length entry
+ }
#ifdef DO_ONCE
if (do_once == 1)
{
- echoinput << "sex: " << gg << " GP: " << gp << " Lmin: " << Lmin(gp) << " Linf: " << L_inf(gp) << " VBK_temp: " << VBK_temp << " VBK@age: " << -VBK(gp) << endl;
+ echoinput << "sex: " << gg << " GP: " << gp << " Lmin: " << Lmin(gp) << " L_max: " << Lmax_temp(gp) << " Linf: " << L_inf(gp) << " VBK@age: " << -VBK(gp) << endl;
if (Grow_type == 2)
echoinput << " Richards: " << Richards(gp) << endl;
if (Grow_type == 8)
- echoinput << " Cessation_decay: " << Richards(gp) << endl;
+ echoinput << "t50 for growth cessation: " << t50 << " using decay_rate: " << Richards(gp) << endl;
}
#endif
+ if (y == styr)
+ {
+ // SS_Label_Info_16.2.4.1 #set up the delta in growth variability across ages if needed
+ // get more CV factors that may depend upon growth parameters, but cannot be time-varying
+ if (CV_const(gp) > 0)
+ {
+ if (CV_depvar_a == 0)
+ {
+ CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (Lmax_temp(gp) - Lmin(gp));
+ }
+ else
+ {
+ CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (AFIX2_forCV - AFIX);
+ }
+ }
+ else
+ {
+ CV_delta(gp) = 0.0;
+ CV_G(gp) = CVLmin(gp); // sets all seasons and whole age range
+ }
+ } // end preliminary calcs
+
// SS_Label_Info_16.2.4 #Loop settlement events because growth starts at time of settlement
g = g_Start(gp); // base platoon
for (settle = 1; settle <= N_settle_timings; settle++)
@@ -421,56 +454,28 @@ FUNCTION void get_growth2(const int y)
g += N_platoon; // increment by N_platoon because only middle platoon has growth modeled
if (use_morph(g) > 0)
{
- if (y == styr)
- {
- switch (Grow_type)
- {
- case 7: // non-parametric
- {
- break;
- }
- default:
- {
- // SS_Label_Info_16.2.4.1 #set up the delta in growth variability across ages if needed
- if (CV_const(gp) > 0)
- {
- if (CV_depvar_a == 0)
- {
- CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (Lmax_temp(gp) - Lmin(gp));
- }
- else
- {
- CV_delta(gp) = (CVLmax(gp) - CVLmin(gp)) / (AFIX2_forCV - AFIX);
- }
- }
- else
- {
- CV_delta(gp) = 0.0;
- CV_G(gp) = CVLmin(gp); // sets all seasons and whole age range
- }
- }
- }
-
- // SS_Label_Info_16.2.4.1.1 #if y=styr, get size-at-age in first subseason of first season of this first year
+ // SS_Label_Info_16.2.4.1.1 #if y=styr, get size-at-age in first subseason of first season of this first year
+ // VBK_seas(0) is the sum of all the season durations * seasonal K adjustments. Normally it has a value of 1.0
+ // but for seasons-as-years, it will be the value equal to the selected duration of the single season
+// Ave_Size(styr, 1, g)(0, first_grow_age(g)) = Lmin(gp);
+ if (y == styr)
+ {
+
switch (Grow_type)
{
case 1:
{
- VBK_temp2 = VBK_temp * VBK_seas(0);
for (a = 0; a <= nages; a++)
{
- // Ave_Size(styr,1,g,a) = Lmin(gp) + (Lmin(gp)-L_inf(gp))* (mfexp(VBK_temp2*(real_age(g,1,a)-AFIX))-1.0);
- Ave_Size(styr, 1, g, a) = L_inf(gp) + (Lmin(gp) - L_inf(gp)) * mfexp(VBK_temp2 * (real_age(g, 1, a) - AFIX));
+ Ave_Size(styr, 1, g, a) = L_inf(gp) + (Lmin(gp) - L_inf(gp)) * mfexp(VBK_work * VBK_seas(0) * (real_age(g, 1, a) - AFIX));
} // done ageloop
break;
}
case 2: // Richards
{
- Ave_Size(styr, 1, g)(0, first_grow_age(g)) = Lmin(gp);
- VBK_temp2 = VBK_temp * VBK_seas(0);
- for (a = first_grow_age(g); a <= nages; a++)
+ for (a = 0; a <= nages; a++)
{
- temp = LinfR + (LminR - LinfR) * mfexp(VBK_temp2 * (real_age(g, 1, a) - AFIX));
+ temp = LinfR + (LminR - LinfR) * mfexp(VBK_work * VBK_seas(0) * (real_age(g, 1, a) - AFIX));
Ave_Size(styr, 1, g, a) = pow(temp, inv_Richards);
} // done ageloop
break;
@@ -484,19 +489,18 @@ FUNCTION void get_growth2(const int y)
case 3: // age-specific K, so need age-by-age calculations
{
ALK_idx = 1;
- //VBK_seas(0) accounts for season duration
+
for (a = 0; a <= nages; a++)
{
- k2 = a - 1;
+ k2 = max(0, a - 1);
if (lin_grow(g, ALK_idx, a) >= -1.0) // linear segment, or first time point beyond AFIX;
{
- Ave_Size(styr, 1, g, a) = Lmin(gp) + (Lmin(gp) - L_inf(gp)) * (mfexp(VBK(gp, 0) * VBK_seas(0) * (real_age(g, 1, a) - AFIX)) - 1.0);
+ Ave_Size(styr, 1, g, a) = Lmin(gp) + (Lmin(gp) - L_inf(gp)) * (mfexp(VBK(gp, k2) * VBK_seas(0) * (real_age(g, 1, a) - AFIX)) - 1.0);
}
else
{
Ave_Size(styr, 1, g, a) = Ave_Size(styr, 1, g, k2) + (mfexp(VBK(gp, k2) * VBK_seas(0)) - 1.0) * (Ave_Size(styr, 1, g, k2) - L_inf(gp));
}
- // echoinput< -997.) // decay rate has been read; uses same code for Richards and standard
+ if (Linf_decay > 0.0001) // decay rate has been read. Should be an estimate of initial Z
{
temp1 = 1.0;
temp4 = 1.0;
temp = current_size;
temp2 = mfexp(-Linf_decay); // cannot use natM or Z because growth is calculated first
- #ifdef DO_ONCE
- if (do_once == 1)
- echoinput << " L_inf " << L_inf(gp) << " size@exactly maxage " << current_size << endl;
- #endif
- if (Grow_type < 3)
- {
- VBK_temp2 = VBK(gp, 0);
- }
- else
- {
- VBK_temp2 = VBK(gp, nages);
- }
- VBK_temp2 = (1.0 - mfexp(VBK_temp2 * VBK_seas(0)));
+
+ dvariable VBK_temp = (1.0 - mfexp(VBK_work * VBK_seas(0)));
for (a = nages + 1; a <= 3 * nages; a++)
{
temp4 *= temp2; // decay numbers at age by exp(-0.xxx)
- current_size += (L_inf(gp) - current_size) * VBK_temp2;
+ current_size += (L_inf(gp) - current_size) * VBK_temp;
temp += temp4 * current_size;
temp1 += temp4; // accumulate numbers to create denominator for mean size calculation
}
@@ -567,15 +560,15 @@ FUNCTION void get_growth2(const int y)
}
Ave_Size(styr, 1, g, nages) = temp / temp1; // this is weighted mean size at nages
}
- else
+ else if (Linf_decay == -998.) // decay code set to -998 which omits growth in the plus group
{
// no adjustment
}
#ifdef DO_ONCE
if (do_once == 1)
- echoinput << " adjusted size at maxage " << Ave_Size(styr, 1, g, nages) << " using decay of: " << Linf_decay << endl;
+ echoinput << " adjusted size in plusgroup: " << Ave_Size(styr, 1, g, nages) << " using decay of: " << Linf_decay << endl;
#endif
- } // end initial year calcs
+ } // end initial year calcs
// SS_Label_Info_16.2.4.2 #loop seasons for growth calculation
for (s = 1; s <= nseas; s++)
@@ -595,12 +588,13 @@ FUNCTION void get_growth2(const int y)
else
add_age = 0; // advance age or not
// growth to next season
- VBK_temp2 = (mfexp(VBK_temp * seasdur(s) * VBK_seas(s)) - 1.0);
- // warning< styr && s == nseas) // do plus group mean size update
{
- if (y > styr && Linf_decay != -998.)
+ if (Linf_decay != -998.)
{
-
+ if (do_once == 1 && timevary_MG_firstyr == y)
+ {
+ warnstream << "plus group mean size is updated in years with time-vary growth beginning in: " << y << "; can turn this off with Linf_decay = -998";
+ write_message (NOTE, 0);
+ }
// 3.24 code
- #ifdef DO_ONCE
if (do_once == 1)
- echoinput << " plus group calc: "
- << " N _entering: " << natage(t, 1, g, nages - 1) << " N_inplus: " << natage(t, 1, g, nages) << " size in: " << Ave_Size(t + 1, 1, g, nages) << " old size: " << plusgroupsize << " ";
- #endif
+ echoinput << niter << " " << y << " plus group calc: "
+ << " N _entering: " << natage(t, 1, g, nages - 1) << " N_inplus: " << natage(t, 1, g, nages) << " size in: " << Ave_Size(t + 1, 1, g, nages) << " old size: " << plusgroupsize << " ";
temp = ((natage(t, 1, g, nages - 1) + 0.01) * Ave_Size(t + 1, 1, g, nages) + (natage(t, 1, g, nages) + 0.01) * plusgroupsize) / (natage(t, 1, g, nages - 1) + natage(t, 1, g, nages) + 0.02);
Ave_Size(t + 1, 1, g, nages) = temp;
- #ifdef DO_ONCE
+
+ #ifdef DO_ONCE
if (do_once == 1 && g == 1)
echoinput << " final_val " << Ave_Size(t + 1, 1, g, nages) << endl;
#endif
@@ -845,7 +825,6 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse
dvariable LminR;
dvariable inv_Richards;
dvariable t50;
- dvariable VBK_temp2;
ALK_idx = (s - 1) * N_subseas + subseas; // note that this changes a global value
for (g = g_Start(1) + N_platoon; g <= gmorph; g += N_platoon) // looping the middle platoons for each sex*gp
@@ -919,10 +898,9 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse
} // done Richards
case 8: // Cessation
{
- // VBK_temp2=-VBK_temp*seasdur(s); // negative to restore positive
// t50 is the calculated inflection age for the decline in K
- dvariable VBK_temp = -VBK(gp, 0);
- t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / (VBK_temp)) - 1.0) / Richards(gp);
+ dvariable VBK_temp = -VBK(gp, 0) * VBK_seas(0);
+ t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / VBK_temp) - 1.0) / Richards(gp);
for (a = 0; a <= nages; a++)
{
// calculate a full year's growth increment, then multiple by seasdur(s)
@@ -968,6 +946,10 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse
break;
}
} // done switch
+ #ifdef DO_ONCE
+ if (do_once == 1)
+ echoinput << " seas.subseas: " << s << "." << subseas <<" g:" << g <<" size: " << Ave_Size(t, subseas, g)(0, min(6, nages)) << " plusgroup: " << Ave_Size(t + 1, 1, g, nages) << endl;
+ #endif
} // end need this platoon
} // done platoon
} // end calc size-at-age at a particular subseason
diff --git a/SS_prelim.tpl b/SS_prelim.tpl
index f02b9636..d270d012 100644
--- a/SS_prelim.tpl
+++ b/SS_prelim.tpl
@@ -1319,6 +1319,16 @@
write_message (WARN, 0);
}
}
+ // check for vulnerability in plusgroup size
+ warning << Ave_Size(styr, 1, g, nages - 1) << " " << 0.99 * Ave_Size(styr, 1, g, nages) << endl;
+ if (Ave_Size(styr, 1, g, nages - 1) < 0.99 * Ave_Size(styr, 1, g, nages))
+ {
+ if ( Linf_decay == -998 && MG_active_firstyr(2) > 0 )
+ {
+ warnstream << "There is greater than 1% growth in plusgroup, timevary growth starts in : " << MG_active_firstyr(2) << " and Linf_decay = -998 disables updating plusgroup size; suggest setting a Linf_decay value like initial Z";
+ write_message (WARN, 0);
+ }
+ }
}
}
}
diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl
index 77df33bd..e5ad3679 100644
--- a/SS_readcontrol_330.tpl
+++ b/SS_readcontrol_330.tpl
@@ -1689,6 +1689,7 @@
// stores years to calc non-constant MG parms (1=natmort; 2=growth; 3=wtlen & fec; 4=recr_dist&femfrac; 5=movement; 6=ageerrorkey; 7=catchmult)
ivector timevary_pass(styr-3,YrMax+1) // extracted column
ivector MG_active(0,7) // 0=all, 1=M, 2=growth 3=wtlen, 4=recr_dist&femfrac, 5=migration, 6=ageerror, 7=catchmult
+ ivector MG_active_firstyr(0,7) // first year in which a MGparm is timevarying
vector env_data_pass(1,2) // holds min-max year with env data
int do_densitydependent;
@@ -1724,6 +1725,7 @@
echoinput << "Now read env, block/trend, and dev adjustments to MGparms " << endl;
timevary_MG.initialize(); // stores years to calc non-constant MG parms (1=natmort; 2=growth; 3=wtlen & fec; 4=recr_dist; 5=movement)
MG_active.initialize();
+ MG_active_firstyr.initialize();
CGD_onoff = 0;
timevary_parm_start_MG = 0;
@@ -1829,6 +1831,7 @@
if (timevary_MG(y, f) > 0)
{
MG_active(f) = 1;
+ if (MG_active_firstyr(2) == 0 ) MG_active_firstyr(f) = y;
timevary_MG(y, 0) = 1; // tracks active status for all MG types
if(timevary_MG_firstyr == YrMax) timevary_MG_firstyr = y; // save for reporting in MSY and spawn_recruit output
}