revise growth code wrt season duration

SS_popdyn.tpl

@@ -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,6 +885,16 @@
   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();
@@ -944,7 +956,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 +981,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 first_grow_age setfirstgrow" << 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 +1002,12 @@
             {
               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 +1016,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) << " " << first_grow_age(g) << " " << setfirstgrowage << endl;
           }
       }
     }
@@ -1689,6 +1711,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 +1747,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 +1853,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

@@ -146,6 +146,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,16 @@ 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 "
+            << " (value should approx initial Z; -999 replicates 3.24 (Z=0.2); -998 or -997 to not allow growth within plus group)" << endl
+            << " #_ -998 also disables time-varying changes in plus group; -997 allows updating time-varying plus group" << endl
+            << " #_ best option is to provide a value for initial Z" << 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