diff --git a/Change_log_for_SS_3.30.xlsx b/Change_log_for_SS_3.30.xlsx index ed80a6bf..e69bea34 100644 Binary files a/Change_log_for_SS_3.30.xlsx and b/Change_log_for_SS_3.30.xlsx differ diff --git a/SS_benchfore.tpl b/SS_benchfore.tpl index 982910ba..7404eef6 100644 --- a/SS_benchfore.tpl +++ b/SS_benchfore.tpl @@ -1651,6 +1651,28 @@ FUNCTION void Get_Forecast() { Make_FishSelex(); // calcs fishery selex by current season, all fleets, current gmorph } + if(N_pred>0) + { + for(f1=1;f1<=N_pred;f1++) + { + f=predator(f1); + pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect + if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s); + + for (gp=1;gp<=N_GP*gender;gp++) + {g=g_Start(gp); // base platoon + for (settle=1;settle<=N_settle_timings;settle++) + { + g+=N_platoon; + int gpi=GP3(g); // GP*gender*settlement + natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f); + surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s)); + surv2(s,gpi)=square(surv1(s,gpi)); + } + } + } + } + } // end of seasonal biology if(s==nseas) {adv_age=1;} else {adv_age=0;} // advance age or not when doing survivorship diff --git a/SS_expval.tpl b/SS_expval.tpl index a78b117d..1cd3fdbb 100644 --- a/SS_expval.tpl +++ b/SS_expval.tpl @@ -285,7 +285,10 @@ FUNCTION void Get_expected_values(const int y,const int t); } case 2: // F rate { - vbio=Hrate(f,t); + if(fleet_type(f)<4) + {vbio=Hrate(f,t);} + else + {vbio=pred_M2(predator_rev(f),t);} break; } case 30: // spawning biomass #30 @@ -448,31 +451,48 @@ FUNCTION void Get_expected_values(const int y,const int t); case(2): // DISCARD_OUTPUT /* SS_Label_46.2 expected discard amount */ { - j=have_data(ALK_time,f,data_type,0); // number of observations - if(j>0) + if(have_data(ALK_time,f,data_type,0)>0) // number of observations { j=have_data(ALK_time,f,data_type,1); // only getting first observation for now - if(catch_ret_obs(f,t)>0.0 || y>endyr) + if(fleet_type(f)<=2) { - if(disc_units(f)==3) // numbers regardless of catchunits for retained catch + if(catch_ret_obs(f,t)>0.0 || y>endyr) + { + if(disc_units(f)==3) // numbers regardless of catchunits for retained catch + { + exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6); + } + else if(catchunits(f)==1) // biomass units for retained and discarded catch + { + exp_disc(f,j)=catch_fleet(t,f,1)-catch_fleet(t,f,3); // discard in biomass + if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,1) + 0.0000001); + } + else // numbers for retained and discarded catch + { + exp_disc(f,j)=catch_fleet(t,f,4)-catch_fleet(t,f,6); // discard in numbers + if(disc_units(f)==2) exp_disc(f,j) /= (catch_fleet(t,f,4) + 0.0000001); + } + if(exp_disc(f,j)<0.0) { + N_warn++; + warning<0) get_natmort(); + if(MG_active(1)>0) + { + get_natmort(); + natM_M1=natM; + } #ifdef DO_ONCE if(do_once==1) cout<<" natmort OK"<=Bmark_Yr(1)&&y<=Bmark_Yr(2)) - { - for (s=1;s<=nseas;s++) - for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) - { - natM_unf(s,gp)+=natM(s,gp); // need nseas to capture differences due to settlement - surv1_unf(s,gp)+=surv1(s,gp); // need nseas to capture differences due to settlement - surv2_unf(s,gp)+=surv2(s,gp); // need nseas to capture differences due to settlement - } - } if(MG_active(4)>0) get_recr_distribution(); if(y>=Bmark_Yr(7)&&y<=Bmark_Yr(8)) @@ -240,6 +234,41 @@ FUNCTION void get_initial_conditions() // SS_Label_Info_23.3.3 #for each platoon, combine size_at_age distribution with length selectivity and weight-at-length to get combined selectivity vectors Make_FishSelex(); } + +// SS_Label_Info_23.3.4 #add predator M2 to M1 to update seasonal M in styr + if(N_pred>0) + { + for(f1=1;f1<=N_pred;f1++) + { + f=predator(f1); + pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect + if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s); + + for (gp=1;gp<=N_GP*gender;gp++) + {g=g_Start(gp); // base platoon + for (settle=1;settle<=N_settle_timings;settle++) + { + g+=N_platoon; + int gpi=GP3(g); // GP*gender*settlement + natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f); + surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s)); + surv2(s,gpi)=square(surv1(s,gpi)); + if(do_once==1) echoinput<=Bmark_Yr(1)&&y<=Bmark_Yr(2)) + { + for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) + { + natM_unf(s,gp)+=natM(s,gp); // need nseas to capture differences due to settlement + surv1_unf(s,gp)+=surv1(s,gp); // need nseas to capture differences due to settlement + surv2_unf(s,gp)+=surv2(s,gp); // need nseas to capture differences due to settlement + } + } } #ifdef DO_ONCE @@ -630,17 +659,12 @@ FUNCTION void get_time_series() // but Maunder's M in get_natmort() may use the fecundity vector, so would be using the most recently calculated Problem?? if(Hermaphro_Option!=0) get_Hermaphro(); } - if(timevary_MG(y,1)>0) get_natmort(); - if(y>=Bmark_Yr(1)&&y<=Bmark_Yr(2)) + if(timevary_MG(y,1)>0) { - for (s=1;s<=nseas;s++) - for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) - { - natM_unf(s,gp)+=natM(s,gp); - surv1_unf(s,gp)+=surv1(s,gp); - surv2_unf(s,gp)+=surv2(s,gp); - } + get_natmort(); + natM_M1=natM; } + if(timevary_MG(y,4)>0) get_recr_distribution(); if(y>=Bmark_Yr(7)&&y<=Bmark_Yr(8)) { @@ -752,8 +776,42 @@ FUNCTION void get_time_series() for (g=1;g<=gmorph;g++) if(use_morph(g)>0) {Make_FishSelex();} + +// SS_Label_Info_24.x.x #add predator M2 inside the yr,seas loop + if(N_pred>0) + { + for(f1=1;f1<=N_pred;f1++) + { + f=predator(f1); + pred_M2(f1,t)=mgp_adj(predparm_pointer(f1)); // base with no seasonal effect + if(nseas>1) pred_M2(f1,t)*=mgp_adj(predparm_pointer(f1)+s); + for (gp=1;gp<=N_GP*gender;gp++) + {g=g_Start(gp); // base platoon + for (settle=1;settle<=N_settle_timings;settle++) + { + g+=N_platoon; + int gpi=GP3(g); // GP*gender*settlement + natM(s,gpi)+=pred_M2(f1,t)*sel_al_3(s,g,f); + surv1(s,gpi)=mfexp(-natM(s,gpi)*seasdur_half(s)); + surv2(s,gpi)=square(surv1(s,gpi)); + if(do_once==1) echoinput<=Bmark_Yr(1)&&y<=Bmark_Yr(2)) + { + for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) + { + natM_unf(s,gp)+=natM(s,gp); + surv1_unf(s,gp)+=surv1(s,gp); + surv2_unf(s,gp)+=surv2(s,gp); + } } + } // SS_Label_Info_24.2.2 #Compute spawning biomass if this is spawning season so recruits could occur later this season // SPAWN-RECR: calc SPB in time series if spawning is at beginning of the season if(s==spawn_seas && spawn_time_seas<0.0001) // compute spawning biomass if spawning at beginning of season so recruits could occur later this season @@ -892,7 +950,7 @@ FUNCTION void get_time_series() // SS_Label_Info_24.3.3 #Do fishing mortality catage_tot.initialize(); - if(catch_seas_area(t,p,0)==1 && fishery_on_off==1) + if( (catch_seas_area(t,p,0)==1 && fishery_on_off==1)) { if(F_Method>1) // not Pope's { @@ -933,7 +991,7 @@ FUNCTION void get_time_series() for (g=1;g<=gmorph;g++) if(use_morph(g)>0) { - Z_rate(t,p,g)=natM(s,GP3(g)); + Z_rate(t,p,g)=natM(s,GP3(g)); // where natM already includes M2 for (f=1;f<=Nfleet;f++) //loop over fishing fleets to get Z if (catch_seas_area(t,p,f)!=0) { @@ -1041,7 +1099,7 @@ FUNCTION void get_time_series() } } //close gmorph loop } // close fishery - } // end continuous F method + } // end continuous F method } else @@ -1129,8 +1187,28 @@ FUNCTION void get_time_series() // SS_Label_Info_24.3.3.4 #No catch or fishery turned off, so set Z=M for (g=1;g<=gmorph;g++) if(use_morph(g)>0) - {Z_rate(t,p,g)=natM(s,GP3(g));} + {Z_rate(t,p,g)=natM(s,GP3(g)); + Zrate2(p,g)=elem_div( (1.-mfexp(-seasdur(s)*Z_rate(t,p,g))), Z_rate(t,p,g)); + } } + if(N_predparms>0) + { + for(f1=1;f1<=N_pred;f1++) + { + f=predator(f1); + for (g=1;g<=gmorph;g++) + if(use_morph(g)>0) + { + catch_fleet(t,f,1)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_1(s,g,f))*Zrate2(p,g); + catch_fleet(t,f,2)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),deadfish_B(s,g,f))*Zrate2(p,g); + catch_fleet(t,f,3)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_2(s,g,f))*Zrate2(p,g); // retained bio + catch_fleet(t,f,4)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_3(s,g,f))*Zrate2(p,g); + catch_fleet(t,f,5)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),deadfish(s,g,f))*Zrate2(p,g); + catch_fleet(t,f,6)+=pred_M2(f1,t)*elem_prod(natage(t,p,g),sel_al_4(s,g,f))*Zrate2(p,g); // retained numbers + catage(t,f,g)=pred_M2(f1,t)*elem_prod(elem_prod(natage(t,p,g),deadfish(s,g,f)),Zrate2(p,g)); + } //close gmorph loop + } + } } //close area loop if(s==1 && save_for_report==1){ Smry_Table(y,1)=totbio; diff --git a/SS_prelim.tpl b/SS_prelim.tpl index 3753e6b7..ea82799b 100644 --- a/SS_prelim.tpl +++ b/SS_prelim.tpl @@ -951,6 +951,30 @@ PRELIMINARY_CALCS_SECTION // SS_Label_Info_6.8.4 #Call fxn get_natmort() echoinput<<"ready to do natmort "<0) + { + for(f1=1;f1<=N_pred;f1++) + { + pred_M2(f1,styr)=MGparm(predparm_pointer(f1)); // base with no seasonal effect + if(nseas==1) + { + for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) {natM(s,gp)+=pred_M2(f1,styr);} + } + else + { + t=styr-1; // resets for each predator + for (s=1;s<=nseas;s++) + { + t++; + for (gp=1;gp<=N_GP*gender*N_settle_timings;gp++) {natM(s,gp)+=pred_M2(f1,styr)*MGparm(predparm_pointer(f1)+s);} + } + } + } + } + natM = value(natM); surv1 = value(surv1); surv2 = value(surv2); diff --git a/SS_readcontrol_330.tpl b/SS_readcontrol_330.tpl index 45f507de..c7c6f5d4 100644 --- a/SS_readcontrol_330.tpl +++ b/SS_readcontrol_330.tpl @@ -600,6 +600,8 @@ // read setup for natmort parameters: LO, HI, INIT, PRIOR, PR_type, CV, PHASE, use_env, use_dev, dev_minyr, dev_maxyr, dev_phase, Block, Block_type int N_MGparm int N_natMparms + int N_predparms + ivector predparm_pointer(1,Nfleet+1) int N_growparms int N_M_Grow_parms int recr_dist_parms @@ -618,6 +620,7 @@ // read natmort setup LOCAL_CALCS N_natMparms=0; + N_predparms=N_pred+N_pred*nseas; natM_5_opt=0; MGparm_point.initialize(); // 0=1Parm; 1=segmented; 2=Lorenzen; 3=agespecific; 4=agespec with seas interpolate; 5=Maunder_M @@ -1044,6 +1047,7 @@ break; } } + switch (Grow_type) { case 1: @@ -1245,6 +1249,32 @@ { ParCount++; ParmLabel+="FracFemale_GP_"+NumLbl(gp); } + predparm_pointer=-1; + for(int pc=1;pc<=N_pred;pc++) + { + ParCount++; + predparm_pointer(pc)=ParCount; // first parm for this predator + onenum=" "; + sprintf(onenum, "%d", pc); + ParmLabel+="M2_pred"+onenum; + Parm_info+="val"; + Parm_minmax.push_back (3); + if(nseas>1) + { + for(s=1;s<=nseas;s++) + { + ParCount++; + onenum2=" "; + sprintf(onenum2, "%d", s); + ParmLabel+="M2_pred"+onenum+"_s"+onenum2; + Parm_info+="val"; + Parm_minmax.push_back (3); + + } + } + } + predparm_pointer(Nfleet+1)=ParCount; + echoinput<<" predparm pointer"<0) mgp_type(AgeKeyParm,N_MGparm)=6; if(catch_mult_pointer>0) mgp_type(catch_mult_pointer,N_MGparm)=7; for(f=frac_female_pointer; f<=frac_female_pointer+N_GP-1;f++) mgp_type(f)=4; + if(N_pred>0) mgp_type(predparm_pointer(1),predparm_pointer(1)+N_predparms-1)=1; echoinput<<"mgparm_type for each parm: 1=M; 2=growth; 3=wtlen,mat,fec,hermo; 4=recr&femfrac; 5=migr; 6=ageerror; 7=catchmult"<1 && need_catch_mult(f)>0) {N_warn++; cout<<"exit with warning"< 0.0 && fleet_type(f)<=2) // excludes survey fleets + if(fleet_area(f)==p && catch_ret_obs(f,t) > 0.0 && fleet_type(f)<=2) // excludes survey and predator fleets { catch_seas_area(t,p,f)=1; catch_seas_area(t,p,0)=1; @@ -888,7 +899,7 @@ disc_maxval(f)=max(disc_maxval(f),obs_disc(f,j)); if( discdata[i](4)<0.0) discdata[i](3)=-fabs( discdata[i](3)); // convert to new format using negative fleet if( discdata[i](3)<0) {yr_disc_use(f,j)=-1;} else {yr_disc_use(f,j)=1;disc_N_fleet_use(f)++;} - if(catch_ret_obs(f,t)<=0.0) + if(fleet_type(f)<3 && catch_ret_obs(f,t)<=0.0) { N_warn++; warning<0 && fishery_on_off==0) SS2out<<";_reported_M_includes_PredM2"; SS2out << endl; SS2out << "Bio_Pattern Sex Yr "<0) { - for (f=1;f<=Nfleet;f++) - if(fleet_type(f)<=2) - if(disc_units(f)>0) SS2out<0) SS2out<0) for (f=1;f<=Nfleet;f++) - if(fleet_type(f)<=2) + if(fleet_type(f)<=2 || fleet_type(f)==4) { for (i=1;i<=disc_N_fleet(f);i++) { @@ -2283,7 +2290,7 @@ FUNCTION void write_bigoutput() SS2out<0) @@ -2310,7 +2317,7 @@ FUNCTION void write_bigoutput() SS2out<0) + if((fleet_type(f)<=2 && Do_Retain(f)>0) || fleet_type(f)==4) for (g=1;g<=gmorph;g++) { if(use_morph(g)>0) @@ -2409,6 +2416,33 @@ FUNCTION void write_bigoutput() if(use_morph(g)>0) {for (s=1;s<=nseas;s++) SS2out<0) + { + SS2out<=3 && Grow_type<=6) { diff --git a/SS_write_ssnew.tpl b/SS_write_ssnew.tpl index 77743018..6b931beb 100644 --- a/SS_write_ssnew.tpl +++ b/SS_write_ssnew.tpl @@ -50,7 +50,7 @@ FUNCTION void write_nudata() report1 << nages<<" #_Nages=accumulator age, first age is always age 0"<< endl; report1 << pop<<" #_Nareas"<