revise growth code wrt season duration

.github/workflows/run-ss3-bootstrap.yml

@@ -45,9 +45,7 @@ jobs:
           sudo add-apt-repository ppa:ubuntu-toolchain-r/test
           sudo apt-get update
           sudo apt-get install --only-upgrade libstdc++6
-          sudo apt-get install -y libcurl4-openssl-dev
-          sudo apt-get install -y libfontconfig1-dev
-          sudo apt-get install -y libharfbuzz-dev libfribidi-dev
+          sudo apt-get install -y libcurl4-openssl-dev libfontconfig1-dev libharfbuzz-dev libfribidi-dev libuv1-dev pkg-config
 
       - name: Setup R
         uses: r-lib/actions/setup-r@v2

.github/workflows/test-r4ss-with-ss3.yml

@@ -45,9 +45,7 @@ jobs:
           sudo add-apt-repository ppa:ubuntu-toolchain-r/test
           sudo apt-get update
           sudo apt-get install --only-upgrade libstdc++6
-          sudo apt-get install -y libcurl4-openssl-dev
-          sudo apt-get install -y libfontconfig1-dev
-          sudo apt-get install -y libharfbuzz-dev libfribidi-dev
+          sudo apt-get install -y libcurl4-openssl-dev libfontconfig1-dev libharfbuzz-dev libfribidi-dev libuv1-dev pkg-config
 
       - name: Setup R
         uses: r-lib/actions/setup-r@v2

SS_param.tpl

@@ -83,7 +83,7 @@ PARAMETER_SECTION
   vector CV_const(1,N_GP*gender)
   matrix mgp_save(styr,YrMax,1,N_MGparm2);
   vector mgp_adj(1,N_MGparm2);
-  matrix Cohort_Growth(styr,YrMax,0,nages)
+  matrix Cohort_Growth(styr,YrMax,0,nages)  // multiplier on the growth increment
   3darray Cohort_Lmin(1,N_GP*gender,styr,YrMax,0,nages)
   vector VBK_seas(0,nseas);
 

SS_popdyn.tpl

@@ -126,7 +126,7 @@ FUNCTION void get_initial_conditions()
       t++;
       for (subseas = 1; subseas <= N_subseas; subseas++) //  do all subseasons in first year
       {
-        get_growth3(y, t, s, subseas); //  in case needed for Lorenzen M
+        get_growth3(styr, t, s, subseas); //  in case needed for Lorenzen M
         Make_AgeLength_Key(s, subseas);
       }
     }
@@ -797,6 +797,7 @@ FUNCTION void get_time_series()
       {
         ALK_subseas_update = 1; // indicate that all ALKs will need re-estimation
         get_growth2(y); //  propagates growth to each season this year and to begin next year
+        t = t_base + 1;
         get_growth3(y, t, 1, 1); //  cleans up the linear growth range for begin of this year
       }
       if (timevary_MG(y, 3) > 0)

SS_prelim.tpl

@@ -1319,6 +1319,15 @@
             write_message (WARN, 0);
           }
         }
+        //  check for vulnerability in plusgroup size
+        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 using -997 or setting a Linf_decay value like initial Z";
+            write_message (WARN, 0);
+          }
+        }
       }
       }
     }

SS_readcontrol_330.tpl

@@ -487,11 +487,13 @@
         {
           ALK_idx = (s - 1) * N_subseas + subseas;
           settle_time = settle_g(g);
-          //   real_age is real age since settlement and is used in growth calculations
-          //   calen_age is real age since the beginning of the year in which spawning occurred
+          //  real_age is real age since settlement and is used in growth calculations
+          //  calen_age is real age since the beginning of the year in which spawning occurred
+          //  except with seas-as-years, sumseas is < 1.0 so calen_age needs to be in terms of season count, not age
+          //  so, divide by sumseas
           for (a = 0; a <= nages; a++)
           {
-            calen_age(g, ALK_idx, a) = r_ages(a) + azero_seas(s) + double(subseas - 1) / double(N_subseas) * seasdur(s);
+            calen_age(g, ALK_idx, a) = r_ages(a) + azero_seas(s) + double(subseas - 1) / double(N_subseas) * seasdur(s) / (sumseas / 12.);
             if (a < Settle_age(settle_time))
             {
               real_age(g, ALK_idx, a) = 0.;
@@ -883,9 +885,17 @@
   number AFIX_plus;
   number Linf_decay;  //  decay factor to calculate mean L at maxage from Linf and the decaying abundance above maxage
                       //  forced equal to 0.20 in 3.24 (which also assumed linear, not VBK, growth)
+                      //  -999 code uses the 3.24 approach
+                      //  -998 code sets no growth within the plus group, even if Linf > Length_at_maxage.  Sets plusgroupsize_update = 0
+                      //  -997 code sets no growth within the plus group, but keeps plusgroupsize_update = 1
+
+  int plusgroupsize_update;  // created internally from value of Linf_decay
+                      //  = 0 keeps mean size at maxage (i.e. plusgroup) constant during time series, even if there is time-varying growth
+                      //  = 1 updates mean size at maxage using weighted average of current mean size and size of incoming cohort, but only in years with timevarying growth
+                      //  note that if Linf is > mean size at maxage, then Z during time series logically should cause reduction in mean size at maxage because there are fewer old fish
+                      //  SO! there needs to be an additional option to update the plus group mean size even if growth is not time-varying
+                      //  but this should be implemented in popdyn, not in biofxn to avoid having to update calculations for all ages
   int do_ageK;
-  ivector first_grow_age(1,gmorph);
-!! first_grow_age.initialize();
 !! k=0;
 !! do_ageK=0;
 !! if(Grow_type<=2 || Grow_type==8) {k=4;}  //  AFIX and AFIX2
@@ -944,7 +954,15 @@
     N_growparms = 2; // for the two CV parameters
     k1 = N_GP * gender; // for reading empirical length_at_age
   }
-
+  if (AFIX2 != 999 && AFIX2 > nages)
+  {
+    warnstream << "AFIX2 (Age post-settlement for L2, aka Amax) must be <= nages, entered value = " << AFIX2;
+    write_message (FATAL, 0);
+    AFIX2 = 999.;
+  }
+
+  plusgroupsize_update = 1;  // update mean size in plus group with weighted mean of current size and growth of incoming cohort.  Only implemented beginning in the year for which there is time-varying growth
+  if (Linf_decay == -998) plusgroupsize_update = 0;  // option which omits the updating
   echoinput << " N_growparms  " << N_growparms << endl;
   AFIX2_forCV = AFIX2;
   if (AFIX2_forCV > nages) AFIX2_forCV = nages;
@@ -961,10 +979,11 @@
   N_M_Grow_parms = N_natMparms + N_growparms;
   lin_grow.initialize();
 
-  echoinput << "g a seas subseas ALK_idx real_age calen_age lin_grow first_grow_age" << endl;
+  echoinput << "g a seas subseas ALK_idx real_age calen_age lin_grow" << endl;
   for (g = 1; g <= gmorph; g++)
     if (use_morph(g) > 0)
     {
+      int setfirstgrowage = -1;
       for (a = 0; a <= nages; a++)
       {
         for (s = 1; s <= nseas; s++)
@@ -981,11 +1000,11 @@
             {
               lin_grow(g, ALK_idx, a) = 1.0; //  at the transition from linear to VBK growth
             }
-            else if (first_grow_age(g) == 0)
+            else if (setfirstgrowage == -1)
             {
               lin_grow(g, ALK_idx, a) = -1.0; //  flag for first age on growth curve beyond AFIX
               if (subseas == N_subseas) {
-                first_grow_age(g) = a;
+                setfirstgrowage = a;
               } //  so that lingrow will be -1 for rest of this season
             }
             else
@@ -994,7 +1013,7 @@
             } //  flag for being in growth curve
 
             if (a < 4) echoinput << g << " " << a << " " << s << " " << subseas << " " << ALK_idx << " " << real_age(g, ALK_idx, a)
-                                 << " " << calen_age(g, ALK_idx, a) << " " << lin_grow(g, ALK_idx, a) << " " << first_grow_age(g) << endl;
+                                 << " " << calen_age(g, ALK_idx, a) << " " << lin_grow(g, ALK_idx, a) << endl;
           }
       }
     }
@@ -1360,9 +1379,9 @@
         }
         case 8:
         {
-          ParmLabel += "L_at_Amin_" + GenderLbl(gg) + GP_Lbl(gp);
+          ParmLabel += "L_at_Age0_" + GenderLbl(gg) + GP_Lbl(gp);
           ParmLabel += "Linf_" + GenderLbl(gg) + GP_Lbl(gp);
-          ParmLabel += "VonBert_K_" + GenderLbl(gg) + GP_Lbl(gp);
+          ParmLabel += "rmax_" + GenderLbl(gg) + GP_Lbl(gp);
           ParmLabel += "Cessation_" + GenderLbl(gg) + GP_Lbl(gp);
           ParCount += 4;
           break;
@@ -1689,6 +1708,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 +1744,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 +1850,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
       }

SS_readdata_330.tpl

@@ -92,6 +92,7 @@
   int eq_yr;
   int bio_yr;
   number sumseas;
+  number sumseas_yr;  // sum of season durations in nits of years.
 
   //  SS_Label_Info_2.1.3 #Set up seasons
   vector seasdur_half(1,nseas);   // half a season
@@ -107,12 +108,14 @@
     seasdur /= sumseas;
     seas_as_year = 0;
     sumseas = 12.0; // to be sure it is exactly 12.
+    sumseas_yr = 1.0;
   }
   else
   {
     seasdur /= 12.;
     seas_as_year = 1;
-    // sumseas will now be used as the duration of the pseudo-year, rather than assuming year has 12 months;
+    sumseas_yr = sumseas / 12.;
+    // sumseas_yr will now be used as the duration of the pseudo-year, rather than assuming year has 12 months;
     if (nseas > 1)
     {
       warnstream << "Error.  Can only have 1 season when during seasons as psuedo-years.";
@@ -146,6 +149,7 @@
   echoinput << subseasdur_delta << " processed subseason duration (frac. of year) " << endl;
   echoinput << " processed subseason cumulative annual time within season " << endl
             << subseasdur << endl;
+  echoinput << sumseas << " sum of season durations; used in seas-as-years" << endl;
   if (seas_as_year == 1)
   {
     warnstream << "Season durations sum to <11.9, so SS3 assumes you are doing years as pseudo-seasons." << endl

SS_write_report.tpl

@@ -3932,10 +3932,6 @@ FUNCTION void write_bigoutput()
       {
         for (y = styr - 3; y <= YrMax; y++)
         {
-          yz = y;
-          if (yz > endyr + 2)
-            yz = endyr + 2;
-          //        if(y==styr-3 || y==styr ||  timevary_MG(yz,2)>0)
           {
             for (s = 1; s <= nseas; s++)
             {
@@ -3948,7 +3944,41 @@ FUNCTION void write_bigoutput()
           }
         }
       }
-    s = 1;
+    if(timevary_MG_firstyr < YrMax)  // timevary growth occurs, so display mean size at age by cohort
+    {
+      int max_t = styr + (YrMax - styr) * nseas + nseas - 1;
+      SS2out << "#" << endl
+             << "mean_size_by_cohort " << endl;
+      SS2out << "Morph YearClass Seas SubSeas" << age_vector << endl;
+      for (g = 1; g <= gmorph; g++)
+      if (use_morph(g) > 0)
+      {
+        for (int yc = styr - nages; yc <= YrMax; yc++)  //  loop year classes
+        {
+//          t = styr + (yc - styr - 1) * nseas;  //  t at age 0 for YC  (assumes birthseason = 1)
+          for (s = 1; s <= nseas; s++)
+          {
+            for (i = 1; i <= N_subseas; i++)
+            {
+              int t_out = styr + (yc - styr - 1) * nseas + s - 1;
+              SS2out << g << " " << yc << " " << s << " " << i;
+              for (a = 0; a <= nages; a++)
+              {
+                t_out += nseas;
+                if (t_out <= max_t && t_out >= styr)
+                {SS2out << " " << Ave_Size(t_out, i, g, a);}
+//                {SS2out << " " << t_out;}  // for checking on the index algorithm
+                else
+                {SS2out << " NA";};
+              }
+              SS2out << endl;
+            }
+          }
+        }
+      }
+    }
+
+      s = 1;
     for (i = 1; i <= gender; i++)
     {
       SS2out << "#" << endl

SS_write_ssnew.tpl

@@ -1908,12 +1908,18 @@ FUNCTION void write_nucontrol()
   {
     report4 << AFIX << " #_Age(post-settlement) for L1 (aka Amin); first growth parameter is size at this age; linear growth below this" << endl
             << AFIX2 << " #_Age(post-settlement) for L2 (aka Amax); 999 to treat as Linf" << endl
-            << Linf_decay << " #_exponential decay of numbers for calc of size in plus group in the initial year (value should approx initial Z; -999 replicates 3.24; -998 to not calc growth above maxage)" << endl;
+            << Linf_decay << " #_exponential decay for growth within plus group and control for time-varying plus group size " << endl
+            << "#_only important when growth does not get near Linf by maxage" << endl
+            << "#_value should approx initial Z; or use a code:  -999 replicates 3.24 (with Z=0.2 and numbers weighted updating in years with time-varying growth)" << endl
+            << "#_-998 ignores growth within plus group in initial year and disables time-varying changes in plus group mean size" << endl
+            << "#_-997 ignores growth within plus group in initial year and enables updating time-varying plus group" << endl
+            << "#" << endl;
+
     report4 << "0  #_placeholder for future growth feature" << endl;
     if (Grow_type >= 3 && Grow_type <= 5)
     {
       report4 << Age_K_count << " # number of K multipliers to read" << endl
-              << Age_K_points << " # ages for K multiplier" << endl;
+              << Age_K_points << " # list ages for K multiplier; in descending order for options 4 and 5" << endl;
     }
   }
   else