add HCR_anchor in forecast.ss

SS_benchfore.tpl

@@ -776,8 +776,6 @@ FUNCTION void Get_Benchmarks(const int show_MSY)
       Equ_SpawnRecr_Result = Equil_Spawn_Recr_Fxn(SR_parm_work, SSB_unf, Recr_unf, SSBpR_unf); //  returns 2 element vector containing equilibrium biomass and recruitment at this SPR
       report5 << " Benchmark SSB, R0, SPR0: " << SSB_unf << " " << Recr_unf << " "  << SSBpR_unf << " equil: " << Equ_SpawnRecr_Result << endl;
     }
-    report5 << 0 << " y: " << y << " Repro_output_by_age_for_morph_1 bench " << fec(1) << endl;
-//    report5 << "Repro_output_by_age_for_morph_1: " << fec(1) << endl;
 
   }
   SR_parm_work(N_SRparm2 + 1) = SSB_unf;
@@ -2211,6 +2209,11 @@ FUNCTION void Get_Forecast()
         break;
       }
     }
+    if (Fcast_Loop_Control(5) <= 1)
+    {HCR_anchor = SSB_unf;}
+    else if (Fcast_Loop_Control(5) ==2)
+    {HCR_anchor = SSB_virgin;}
+    report5 << "Control_rule_anchor_approach: " << Fcast_Loop_Control(5) << " HCR_anchor: " << HCR_anchor << endl;
     report5 << "#" << endl;
   }
 
@@ -2227,7 +2230,7 @@ FUNCTION void Get_Forecast()
 
   for (int Fcast_Loop1 = 1; Fcast_Loop1 <= jloop; Fcast_Loop1++) //   for different forecast conditions
   {
-    report5 << Fcast_Loop1 << " y: " << 0 << " Repro_output_by_age_for_morph_1 top_forecast: " << fec(1) << endl;
+//    report5 << Fcast_Loop1 << " y: " << 0 << " Repro_output_by_age_for_morph_1 top_forecast: " << fec(1) << endl;
 
     switch (Fcast_Loop1) //  select which ABC_loops to use
     {
@@ -2315,7 +2318,7 @@ FUNCTION void Get_Forecast()
         }
       }
     }
-    report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 endyr: " << fec(1) << endl;
+//    report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 endyr: " << fec(1) << endl;
     for (y = endyr + 1; y <= YrMax; y++)
     {
       t_base = styr + (y - styr) * nseas - 1;
@@ -2505,7 +2508,7 @@ FUNCTION void Get_Forecast()
                   Wt_Age_mid(s, g) = ALK(ALK_idx, g) * wt_len(s, GP(g)); // use for fisheries with no size selectivity
                 }
             }
-    report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 annual: " << fec(1) << endl;
+//    report5 << Fcast_Loop1 << " y: " << y << " updated_Repro_output_by_age_for_morph_1 annual: " << fec(1) << endl;
             Wt_Age_t(t, 0) = Wt_Age_beg(s);
             for (g = 1; g <= gmorph; g++)
               if (use_morph(g) > 0)
@@ -2659,9 +2662,9 @@ FUNCTION void Get_Forecast()
           }
           else if (ABC_Loop == 2 && s == 1) // Calc the buffer in season 1, will use last year's spawnbio if multiseas and spawnseas !=1
           {
-            temp = SSB_unf;
-            join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * temp)));
-            join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * temp)));
+
+            join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * HCR_anchor)));
+            join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * HCR_anchor)));
 
             switch (HarvestPolicy)
             {
@@ -2674,8 +2677,8 @@ FUNCTION void Get_Forecast()
                 // ramp scales catch as f(B) and buffer (H4010_scale) applied to F
                 {
                   ABC_buffer(y) = H4010_scale_vec(y) *
-                          ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
-                                      + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+                          ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+                                      + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
                                       ) *
                           (join2) // scale combo
                       +
@@ -2686,8 +2689,8 @@ FUNCTION void Get_Forecast()
                 // ramp scales F as f(B) and buffer (H4010_scale) applied to F
                 {
                   ABC_buffer(y) = H4010_scale_vec(y) *
-                          ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
-                                      + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+                          ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+                                      + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
                                       ) *
                           (join2) // scale combo
                       +
@@ -2698,8 +2701,8 @@ FUNCTION void Get_Forecast()
                 // ramp scales catch as f(B) and buffer (H4010_scale) applied to catch
                 {
                   ABC_buffer(y) = 1.0 *
-                          ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
-                                      + (0.0001 + (1.0 - 0.0001) * (H4010_top * temp / SSB_current) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+                          ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+                                      + (0.0001 + (1.0 - 0.0001) * (H4010_top * HCR_anchor / SSB_current) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
                                       ) *
                           (join2) // scale combo
                       +
@@ -2710,8 +2713,8 @@ FUNCTION void Get_Forecast()
                 // ramp scales F as f(B) and buffer (H4010_scale) applied to catch
                 {
                   ABC_buffer(y) = 1.0 *
-                          ((0.0001 * SSB_current / (H4010_bot * temp)) * (join1) // low
-                                      + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * temp) / (H4010_top * temp - H4010_bot * temp)) * (1.0 - join1) // curve
+                          ((0.0001 * SSB_current / (H4010_bot * HCR_anchor)) * (join1) // low
+                                      + (0.0001 + (1.0 - 0.0001) * (SSB_current - H4010_bot * HCR_anchor) / (H4010_top * HCR_anchor - H4010_bot * HCR_anchor)) * (1.0 - join1) // curve
                                       ) *
                           (join2) // scale combo
                       +
@@ -3514,7 +3517,7 @@ FUNCTION void Get_Forecast()
             f = fish_fleet_area(0, ff);
             if (fleet_type(f) == 1)
             {
-              if (ABC_Loop == 2 && HarvestPolicy >= 3)
+              if (ABC_Loop == 2 && HarvestPolicy >= 3)  //  alternative ABC_buffer approach
               {
                 catch_fleet(t, f) *= H4010_scale_vec(y);
               }

SS_param.tpl

@@ -235,6 +235,7 @@ PARAMETER_SECTION
   init_bounded_vector Fcast_recruitments(recdev_end+1,s,recdev_LO,recdev_HI,Fcast_recr_PH2)
   init_bounded_vector Fcast_impl_error(endyr+1,j,-1,1,k)
   vector ABC_buffer(endyr+1,YrMax);
+  number HCR_anchor  //  basis (denominator) bor inflection in control rule
 
 //  SPAWN-RECR:   define some spawning biomass and recruitment entities
   number SSB_virgin

SS_readdata_330.tpl

@@ -4344,8 +4344,8 @@
       "even when the base is set to the mean of earlier recruitments" << endl;
     }
 
-    echoinput << Fcast_Loop_Control(5) << " #echo: loop control 5 not used" << endl;
-
+    echoinput << Fcast_Loop_Control(5) << " #control rule anchor: 1=unfished_benchmark_SSB(old_approach), 2=virgin_SSB " << endl;
+    
     echoinput << "#next enter year in which Fcast loop 3 caps and allocations begin to be applied" << endl;
     *(ad_comm::global_datafile) >> Fcast_Cap_FirstYear;
     echoinput << Fcast_Cap_FirstYear << " # echoed value" << endl;

SS_write_report.tpl

@@ -4687,7 +4687,7 @@ FUNCTION void SPR_profile()
     Wt_Age_mid(s) = Wt_Age_t(t, -1);
     if (s == spawn_seas)
       fec = Wt_Age_t(t, -2);
-    report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
+//    report5 << 0 << " y: " << y << " updated_Repro_output spr/ypr: " << fec(1) << endl;
   }
 
   SS2out << "SPRloop Iter Bycatch Fmult F_report SPR YPR_dead YPR_dead*Recr YPR_ret*Recr Revenue Cost Profit SSB Recruits SSB/Bzero Tot_Catch ";
@@ -5020,9 +5020,9 @@ FUNCTION void Global_MSY()
       SS2out << "Actual ";
       show_MSY = 2; //  invokes just brief output in benchmark
       did_MSY = 0;
-    report5 << 0 << " y: " << y << " updated_Repro_output global_1: " << fec(1) << endl;
+//    report5 << 0 << " y: " << y << " updated_Repro_output global_1: " << fec(1) << endl;
       Get_Benchmarks(show_MSY);
-    report5 << 0 << " y: " << y << " updated_Repro_output global_2: " << fec(1) << endl;
+//    report5 << 0 << " y: " << y << " updated_Repro_output global_2: " << fec(1) << endl;
       did_MSY = 0;
     }
   }

SS_write_ssnew.tpl

@@ -1631,6 +1631,7 @@ FUNCTION void write_nucontrol()
   NuFore << HarvestPolicy << " # Control rule method (0: none; 1: ramp does catch=f(SSB), buffer on F; 2: ramp does F=f(SSB), buffer on F; 3: ramp does catch=f(SSB), buffer on catch; 4: ramp does F=f(SSB), buffer on catch) " << endl;
   NuFore << "# values for top, bottom and buffer exist, but not used when Policy=0" << endl;
   NuFore << H4010_top_rd << " # Control rule inflection for constant F (as frac of Bzero, e.g. 0.40); must be > control rule cutoff, or set to -1 to use Bmsy/SSB_unf " << endl;
+  NuFore << "# Also see HCR_anchor below" << endl;
   NuFore << H4010_bot << " # Control rule cutoff for no F (as frac of Bzero, e.g. 0.10) " << endl;
   NuFore << H4010_scale_rd << " # Buffer:  enter Control rule target as fraction of Flimit (e.g. 0.75), negative value invokes list of [year, scalar] with filling from year to YrMax " << endl;
   if (H4010_scale_rd < 0)
@@ -1654,7 +1655,7 @@ FUNCTION void write_nucontrol()
   {
     NuFore << Fcast_Loop_Control(4) << " # multiplier on base recruitment " << endl;
   }
-  NuFore << Fcast_Loop_Control(5) << " # not used" << endl << "#" << endl;
+  NuFore << Fcast_Loop_Control(5) << " # HCR_anchor: 0 or 1 uses unfished benchmark SSB (old hardwired approach), 2 = virgin SSB" << endl << "#" << endl;
 
   NuFore << Fcast_Cap_FirstYear << "  # FirstYear for caps and allocations (should be after years with fixed inputs) " << endl;