Skip to content

Commit 3716a97

Browse files
Rick-Methot-NOAARick-Methot-NOAA
committed
add cohort devs and prevent shrink to Cessation
1 parent 6788c0e commit 3716a97

1 file changed

Lines changed: 54 additions & 24 deletions

File tree

‎SS_biofxn.tpl‎

Lines changed: 54 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -471,6 +471,7 @@ FUNCTION void get_growth2(const int y)
471471
}
472472
case 2: // Richards
473473
{
474+
warning<<LminR << " " << LinfR << " " << Richards(gp) << " " << inv_Richards << endl;
474475
for (a = 0; a <= nages; a++)
475476
{
476477
temp = LinfR + (LminR - LinfR) * mfexp(VBK_work * VBK_seas(0) * (real_age(g, 1, a) - AFIX));
@@ -651,9 +652,12 @@ FUNCTION void get_growth2(const int y)
651652
} // where add_age =1 if s=nseas, else 0 (k2 assignment could be in a matrix so not recalculated
652653
if (lin_grow(g, ALK_idx, a) == -2.0)
653654
{
654-
temp = pow(Ave_Size(t, 1, g, a), Richards(gp));
655+
temp = pow(Ave_Size(t, 1, g, a), Richards(gp)); // transform current size
655656
t2 = temp - LinfR; // remaining growth potential in transformed units
656-
temp1 = Ave_Size(t, 1, g, a) + pow((VBK_by_seas * t2), inv_Richards) * Cohort_Growth(y, a);
657+
temp1 = pow(temp + (VBK_by_seas * t2), inv_Richards); // new size without cohort growth dev
658+
// if(a<4 && g==1) echoinput << y<<" "<<a <<" lingrow: "<<lin_grow(g, ALK_idx, a)<<" start "<<Ave_Size(t, 1, g, a)<<" t2 "<<t2<<" result "<<temp1;
659+
temp1 = Ave_Size(t, 1, g, a) + (temp1 - Ave_Size(t, 1, g, a)) * Cohort_Growth(y, a);
660+
// if(a<4 && g==1) echoinput << " after CGD " << temp1 << endl;
657661
if (a == nages && s == nseas)
658662
{
659663
plusgroupsize = temp1;
@@ -664,8 +668,10 @@ FUNCTION void get_growth2(const int y)
664668
}
665669
else if (lin_grow(g, ALK_idx, a) == -1.0) // first time point beyond AFIX; lin_grow will stay at -1 for all remaining subseas of this season
666670
{
667-
temp = (LminR - LinfR) * (mfexp(VBK_work * sumseas_yr * VBK_seas(s) * (real_age(g, ALK_idx2, k2) - AFIX)) - 1.0); // growth potential in transform units
668-
Ave_Size(t + 1, 1, g, k2) = Cohort_Lmin(gp, y, a) + pow(temp, inv_Richards) * Cohort_Growth(y, a); // growth in real units
671+
t2 = (LminR - LinfR) * (mfexp(VBK_work * sumseas_yr * VBK_seas(s) * (real_age(g, ALK_idx2, k2) - AFIX)) - 1.0); // growth potential in transform units
672+
temp1 = pow((LminR + t2), inv_Richards); // new size in real units, w/o cohort_growth
673+
Ave_Size(t + 1, 1, g, k2) = Lmin(gp) + (temp1 - Lmin(gp)) * Cohort_Growth(y, a);
674+
// if(a<4 && g==1) echoinput << y<<" "<<a <<" lingrow: "<<lin_grow(g, ALK_idx, a)<<" first "<<Ave_Size(t, 1, g, a)<<" t2 "<<t2<<" result "<<Ave_Size(t + 1, 1, g, k2)<<endl;
669675
}
670676
else // in linear phase for subseas
671677
{
@@ -678,8 +684,8 @@ FUNCTION void get_growth2(const int y)
678684
{
679685
join1 = 1.0 / (1.0 + mfexp(100. * t2 ));
680686
Ave_Size(t + 1, 1, g, k2) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t + 1, 1, g, k2) * (1 - join1);
681-
if (do_once == 1 && g == 1)
682-
echoinput << y << " " << a << " prevent shrink-G2 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t + 1, 1, g, k2) << endl;
687+
// if (do_once == 1 && g == 1)
688+
// echoinput << y << " " << a << " prevent shrink-G2 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t + 1, 1, g, k2) << endl;
683689
}
684690
}
685691
} // done ageloop
@@ -699,19 +705,24 @@ FUNCTION void get_growth2(const int y)
699705
if (timevary_MG(y, 2) > 0)
700706
{
701707
t2 = Ave_Size(t + 1, 1, g, k2) - Ave_Size(t, 1, g, a); // growth increment
708+
t2 *= Cohort_Growth(y, a); // adjust increment with cohort growth dev
702709
if (t2 < 0.0) // trap to prevent decrease in size-at-age
703710
{
704711
join1 = 1.0 / (1.0 + mfexp(100. * t2 ));
705-
Ave_Size(t + 1, 1, g, k2) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t + 1, 1, g, k2) * (1 - join1);
706-
if (do_once == 1 && g == 1)
707-
echoinput << y << " " << a << " prevent shrink-G2 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t + 1, 1, g, k2) << endl;
712+
t2 *= (1. - join1); // trap to prevent negative growth increment
713+
// Ave_Size(t + 1, 1, g, k2) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t + 1, 1, g, k2) * (1 - join1);
714+
// if (do_once == 1 && g == 1)
715+
// echoinput << y << " " << a << " prevent shrink-G2 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t + 1, 1, g, k2) << endl;
708716
}
717+
Ave_Size(t + 1, 1, g, k2) = Ave_Size(t, 1, g, a) + t2;
709718
}
719+
Ave_Size(t + 1, 1, g, 0) = Cohort_Lmin(gp, y, a);
710720
if (a == nages && s == nseas)
711721
{
712722
plusgroupsize = Ave_Size(t, 1, g, nages) +
713723
(-VBK_work - (-VBK_work / Richards(gp)) * (log(exp(Richards(gp) * (real_age(g, 1, a) + 1 - t50)) + 1) - log(exp(Richards(gp) * (real_age(g, 1, a) - t50)) + 1))) * seasdur(s);
714724
}
725+
// if(a<4 && g==1) echoinput <<" G2:-T50 "<<t50<<" VBK: "<<-VBK_work<<" start_size "<<Ave_Size(t, 1, g, a)<<" incre " << t2<<" result "<<Ave_Size(t + 1, 1, g, k2)<<endl;
715726
} // done ageloop
716727
break;
717728
}
@@ -887,60 +898,79 @@ FUNCTION void get_growth3(const int y, const int t, const int s, const int subse
887898
inv_Richards = 1.0 / Richards(gp);
888899
// uses VBK(nages) because age-specific K not allowed
889900
// and Cohort_Lmin has already had the power function applied
901+
echoinput << t << " G3-before seas.sub: " << s << "." << subseas <<" g:" << g <<" size: " << Ave_Size(t, subseas, g)(0, min(6, nages)) << " plusgroup: " << Ave_Size(t, 1, g, nages) << endl;
890902
for (a = 0; a <= nages; a++)
891903
{
904+
// if(a<4 && g==1) echoinput << y<<" G3: "<<a <<" lingrow: "<<lin_grow(g, ALK_idx, a) <<" ss " << subseas << " dur "<<subseasdur(s, subseas)<<endl;
905+
892906
if (lin_grow(g, ALK_idx, a) == -2.0) // so doing growth curve
893907
{
908+
// temp = pow(Ave_Size(t, 1, g, a), Richards(gp));
909+
// t2 = temp - LinfR; // remaining growth potential
910+
// temp += (mfexp(VBK(gp, 0) * subseasdur(s, subseas) * VBK_seas(s)) - 1.0) * t2 * Cohort_Growth(y, a);
911+
// Ave_Size(t, subseas, g, a) = pow(temp, inv_Richards);
912+
894913
temp = pow(Ave_Size(t, 1, g, a), Richards(gp));
895914
t2 = temp - LinfR; // remaining growth potential in transformed units
896-
temp1 = (mfexp(VBK(gp, 0) * VBK_seas(s) * subseasdur(s, subseas)) - 1.0) * t2;
897-
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) + pow(temp1, inv_Richards) * Cohort_Growth(y, a);
915+
temp1 = (mfexp(VBK(gp, 0) * VBK_seas(s) * subseasdur(s, subseas)) - 1.0) * t2; // growth increment in transform units w/o cohort growth dev
916+
Ave_Size(t, subseas, g, a) = pow((temp + temp1), inv_Richards); // new size
917+
// if(a<4 && g==1) echoinput <<" grow "<<Ave_Size(t, 1, g, a)<<" result "<<Ave_Size(t, subseas, g, a)<<" t2 "<<t2 <<" temp1 " << temp1<<" VB*time " <<(VBK(gp, 0) * VBK_seas(s) * subseasdur(s, subseas))<<endl;
898918
}
899919
else if (lin_grow(g, ALK_idx, a) >= 0.0) // in linear phase for subseas
900920
{
901921
Ave_Size(t, subseas, g, a) = len_bins(1) + lin_grow(g, ALK_idx, a) * (Cohort_Lmin(gp, y, a) - len_bins(1));
922+
// if(a<4 && g==1) echoinput <<" lin "<<Ave_Size(t, 1, g, a)<<" result "<<Ave_Size(t, subseas, g, a)<<" t2 "<<t2 <<" temp1 " << temp1<<" VB*time " <<(VBK(gp, 0) * VBK_seas(s) * subseasdur(s, subseas))<<endl;
902923
}
903924
else if (lin_grow(g, ALK_idx, a) == -1.0) // first time point beyond AFIX; lin_grow will stay at -1 for all remaining subseas of this season
904925
{
905-
// temp=Cohort_Lmin(gp,y,a) + (Cohort_Lmin(gp,y,a)-LinfR)*
906-
temp = (LminR - LinfR) * (mfexp(VBK(gp, 0) * sumseas_yr * VBK_seas(s) * (real_age(g, ALK_idx, a) - AFIX)) - 1.0);
907-
Ave_Size(t, subseas, g, a) = Cohort_Lmin(gp, y, a) + pow(temp, inv_Richards) * Cohort_Growth(y, a); // growth in real units
926+
temp = LminR + (LminR - LinfR) * (mfexp(VBK(gp, 0)* VBK_seas(s) * sumseas_yr * (real_age(g, ALK_idx, a) - AFIX)) - 1.0);
927+
// if(a<4 && g==1) echoinput <<" temp " << temp;
928+
temp = LinfR + (LminR - LinfR) * mfexp(VBK(gp, 0)* VBK_seas(s) * sumseas_yr * (real_age(g, ALK_idx, a) - AFIX));
929+
// if(a<4 && g==1) echoinput <<" alt " << temp << endl;
930+
Ave_Size(t, subseas, g, a) = pow(temp, inv_Richards); // growth in real units
931+
// if(a<4 && g==1) echoinput <<" first "<<Ave_Size(t, 1, g, a)<<" result "<<Ave_Size(t, subseas, g, a)<<" temp " << temp<<" LminR " <<LminR <<" " << LinfR<<endl;
908932
}
933+
/*
909934
if (timevary_MG(y, 2) > 0)
910935
{
911-
t2 = Ave_Size(t, subseas, g, a) - Ave_Size(t, 1, g, a); // growth increment
936+
t2 = (Ave_Size(t, subseas, g, a) - Ave_Size(t, 1, g, a)) * Cohort_Growth(y, a); // growth increment with cohort growth dev
937+
if (do_once == 1 && g == 1 && a<4)
938+
echoinput << " ; delta " <<Ave_Size(t, subseas, g, a) <<" "<<Ave_Size(t, 1, g, a) << " t2: "<< t2<<endl;
912939
if (t2 < 0.0) // trap to prevent decrease in size-at-age
913940
{
914941
join1 = 1.0 / (1.0 + mfexp(100. * t2 ));
915-
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t, subseas, g, a) * (1 - join1);
942+
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t, subseas, g, a) * (1.0 - join1);
916943
if (do_once == 1 && g == 1)
917-
echoinput << y << " " << a << " prevent shrink-G3 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t, subseas, g, a) << endl;
944+
echoinput << y << " " << a << " prevent shrink-G3 " << t2 << " " << join1 << " start " << Ave_Size(t, 1, g, a) << " after " << Ave_Size(t, subseas, g, a) << endl;
918945
}
919946
}
947+
*/
920948
} // done ageloop
949+
echoinput << t << " G3-after seas.sub: " << s << "." << subseas <<" g:" << g <<" size: " << Ave_Size(t, subseas, g)(0, min(6, nages)) << " plusgroup: " << Ave_Size(t, 1, g, nages) << endl;
921950
break;
922951
} // done Richards
923952
case 8: // Cessation
924953
{
954+
dvariable VBK_temp = -VBK(gp, 0) * VBK_seas(0) * sumseas_yr;
925955
// t50 is the calculated inflection age for the decline in K
926-
dvariable VBK_temp = -VBK(gp, 0) * VBK_seas(0);
927956
t50 = log(exp((L_inf(gp) - Lmin(gp)) * Richards(gp) / VBK_temp) - 1.0) / Richards(gp);
928957
for (a = 0; a <= nages; a++)
929958
{
930959
// calculate a full year's growth increment, then multiple by seasdur(s)
931-
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) +
932-
(VBK_temp - (VBK_temp / Richards(gp)) * (log(exp(Richards(gp) * (real_age(g, ALK_idx, a) + 1 - t50)) + 1) - log(exp(Richards(gp) * (real_age(g, ALK_idx, a) - t50)) + 1))) * subseasdur(s, subseas);
960+
t2 = (VBK_temp - (VBK_temp / Richards(gp)) * (log(exp(Richards(gp) * (real_age(g, ALK_idx, a) + 1 - t50)) + 1) - log(exp(Richards(gp) * (real_age(g, ALK_idx, a) - t50)) + 1))) * subseasdur(s, subseas);
933961
if (timevary_MG(y, 2) > 0)
934962
{
935-
t2 = Ave_Size(t, subseas, g, a) - Ave_Size(t, 1, g, a); // growth increment
963+
t2 *= Cohort_Growth(y, a); // growth increment
936964
if (t2 < 0.0) // trap to prevent decrease in size-at-age
937965
{
938966
join1 = 1.0 / (1.0 + mfexp(100. * t2 ));
939-
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) * join1 + Ave_Size(t, subseas, g, a) * (1 - join1);
940-
if (do_once == 1 && g == 1)
941-
echoinput << y << " " << a << " prevent shrink-G3 " << t2 << " " << join1 << " " << Ave_Size(t, 1, g, a) << " " << Ave_Size(t, subseas, g, a) << endl;
967+
t2 *= (1. - join1); // trap to prevent negative growth increment
968+
// if (do_once == 1 && g == 1)
969+
// echoinput << y << " " << a << " prevent shrink-G3 " << t2 << " " << join1 << endl;
942970
}
943971
}
972+
Ave_Size(t, subseas, g, a) = Ave_Size(t, 1, g, a) + t2;
973+
// if(a<4 && g==1) echoinput <<" G3:-T50 "<<t50<<" VBK: "<<VBK_temp<<" start_size "<<Ave_Size(t, 1, g, a)<<" incre " << t2<<" result "<<Ave_Size(t, subseas, g, a)<<endl;
944974
} // done ageloop
945975
break;
946976
}

0 commit comments

Comments
 (0)