///////////////////////////////////////////////////////////////////
// JJM.tpl (SPRFMO Joint jack mackerel model)
// Naming Conventions:
//
//  GENERAL:
//    styr, endyr begining year and ending year of model (catch data available)
//    nages       number of age groups considered
//    nyrs_        number of observations available to specific data set
//
//  DATA SPECIFIC:

//    catch_bio   Observed catch biomass
//    fsh        fishery data
//
//  Define surveys
//    nsrv        number of surveys
//  Index values
//    nyrs_srv      Number of years of survey index value (annual)
//    yrs_srv        Years of survey index value (annual)
//    obs_srv        Observed survey index value (annual)
//    obs_se_srv    Observed survey standard errors (annual)
//  Age-comp values
//    nyrs_srv_age  Number of years survey age data available
//    yrs_srv_age   Years of survey age value (annual)
//    oac_srv       Observed age comp from surveys
//    n_sample_srv_age    Observed age comp sample sizes from surveys
//
//    eac_srv       Expected age comp from surveys
//
//    sel_srv       selectivity for egg production index
//
//    pred_srv ...
//
//    oac_fsh      Observed age comp from survey
//    obs_srv_size  Observed size comp from survey
//
//    pred_fsh_age    Predicted age comp from fishery data
//    eac_fsh            Expected age comp for fishery data (only years where data available)
//    eac_ ...
//
//    pred_tmp_srv   Predicted index value for trawl survey
//
//    sel_fsh    selectivity for fishery                
//  
//     sel_ch indicates time-varying selectivity change
//  
//    Add bit for historical F
//
//////////////////////////////////////////////////////////////////////////////
 // To ADD/FIX:
 //   parameterization of steepness to work the same (wrt prior) for ricker and bholt
 //   splines for selectivity
//////////////////////////////////////////////////////////////////////////////

DATA_SECTION
  !!version_info+="SPRFMO_Joint_Jack_mackerel_assessment_model_version_0.1;August_2010";
  int iseed 
  int nnodes_tmp;
  !!CLASS ofstream mceval("mceval.dat")
  !!long int lseed=iseed;
  !!CLASS random_number_generator rng(iseed);
  
  int oper_mod
  int mcmcmode
  int mcflag

  !! oper_mod = 0;
  !! mcmcmode = 0;
  !! mcflag   = 1;
 LOCAL_CALCS
  write_input_log<<version_info<<endl;
  tmpstring=adprogram_name + adstring(".dat");
  // if (argc > 1)
  // {
    int on=0;
    if ( (on=option_match(argc,argv,"-ind"))>-1)
    {
      if (on>argc-2 | argv[on+1][0] == '-') 
      { 
        cerr << "Invalid input data command line option"
           " -- ignored" << endl;  
      }
      else
      {
        cntrlfile_name = adstring(argv[on+1]);
      }
    }
    else
    {
        cntrlfile_name =   tmpstring;
    }
    if ( (on=option_match(argc,argv,"-om"))>-1)
    {
      oper_mod  = 1;
      cout<<"Got to operating model option "<<oper_mod<<endl;
    }
    if ( (on=option_match(argc,argv,"-mcmc"))>-1)
    {
      mcmcmode = 1;
    }
  // }
  global_datafile= new cifstream(cntrlfile_name);
  if (!global_datafile)
  {
  }
  else
  {
    if (!(*global_datafile))
    {
      delete global_datafile;
      global_datafile=NULL;
    }
  }
 END_CALCS
 // Read in "name" of this model...
  !! *(ad_comm::global_datafile) >>  datafile_name; // First line is datafile (not used by this executable)
  !! *(ad_comm::global_datafile) >>  model_name; 
  !! ad_comm::change_datafile_name(datafile_name);
  init_int styr
  init_int endyr
  init_int rec_age
  init_int oldest_age
  int nages
  !!  nages = oldest_age - rec_age + 1;
  int styr_rec
  int styr_sp
  int endyr_sp
  int nyrs
  !!  nyrs = endyr - styr + 1;
  !! styr_rec = (styr - nages) + 1;     // First year of recruitment
  !! styr_sp  = styr_rec - rec_age ;  // First year of spawning biomass  
  !! endyr_sp = endyr   - rec_age - 1;// endyr year of (main) spawning biomass
  vector yy(styr,endyr);
  !! yy.fill_seqadd(styr,1) ;
  vector aa(1,nages);
  !! aa.fill_seqadd(rec_age,1) ;

  int junk;
// Fishery specifics
  init_int nfsh                                   //Number of fisheries
  imatrix pfshname(1,nfsh,1,2)
  init_adstring fshnameread;
 LOCAL_CALCS
  for(k=1;k<=nfsh;k++) 
  {
    pfshname(k,1)=1; 
    pfshname(k,2)=1;
  }    // set whole array to equal 1 in case not enough names are read
  adstring_array CRLF;   // blank to terminate lines
  CRLF+="";
  k=1;
  for(i=1;i<=strlen(fshnameread);i++)
  if(adstring(fshnameread(i))==adstring("%")) {
    pfshname(k,2)=i-1; 
    k++;  
    pfshname(k,1)=i+1;
  }
  pfshname(nfsh,2)=strlen(fshnameread);
  for(k=1;k<=nfsh;k++)
  {
    fshname += fshnameread(pfshname(k,1),pfshname(k,2))+CRLF(1);
  }
  log_input(datafile_name);
  log_input(model_name);
  log_input(styr);
  log_input(endyr);
  log_input(rec_age);
  log_input(oldest_age);
  log_input(nfsh);
  log_input(fshname);
 END_CALCS
  init_matrix catch_bio(1,nfsh,styr,endyr)
  init_matrix catch_bio_sd(1,nfsh,styr,endyr)   // Specify catch-estimation precision
  !! log_input(catch_bio);
  !! log_input(catch_bio_sd);
  matrix catch_bio_lsd(1,nfsh,styr,endyr)      //Catch biomass standard errors (for lognormal)
  matrix catch_bio_lva(1,nfsh,styr,endyr)      //Catch biomass standard errors (for lognormal)
  !! catch_bio_lsd = sqrt(log(square(catch_bio_sd) + 1.));
  !! catch_bio_lva = log(square(catch_bio_sd) + 1.);
  matrix catch_bioT(styr,endyr,1,nfsh)
  !! catch_bioT = trans(catch_bio);
  // !! cout << catch_bio<<endl;exit(1);

  init_ivector nyrs_fsh_age(1,nfsh)
  !! log_input(nyrs_fsh_age);

  init_imatrix yrs_fsh_age(1,nfsh,1,nyrs_fsh_age)
  !! log_input(yrs_fsh_age);

  init_matrix n_sample_fsh_age(1,nfsh,1,nyrs_fsh_age)    //Years of survey index value (annual)
  !! log_input(n_sample_fsh_age);

  init_3darray oac_fsh(1,nfsh,1,nyrs_fsh_age,1,nages)
    3darray sim_catage(1,nfsh,1,nyrs_fsh_age,1,nages);
  !! log_input(oac_fsh);

  init_3darray wt_fsh(1,nfsh,styr,endyr,1,nages)  //values of Survey proportions at age
  !! log_input(wt_fsh);

//  Define surveys
  init_int nsrv                                   //number of surveys
  !! log_input(nsrv);
  int nfsh_and_srv
  !! nfsh_and_srv = nfsh+nsrv;
  imatrix psrvname(1,nsrv,1,2)
  init_adstring srvnameread;
 LOCAL_CALCS
  for(k=1;k<=nsrv;k++) 
  {
    psrvname(k,1)=1; 
    psrvname(k,2)=1;
  }    // set whole array to equal 1 in case not enough names are read
  k=1;
  for(i=1;i<=strlen(srvnameread);i++)
  if(adstring(srvnameread(i))==adstring("%")) {
    psrvname(k,2)=i-1; 
    k++;  
    psrvname(k,1)=i+1;
  }
  psrvname(nsrv,2)=strlen(srvnameread);
  for(k=1;k<=nsrv;k++)
  {
    srvname += srvnameread(psrvname(k,1),psrvname(k,2))+CRLF(1);
  }
  log_input(srvname);
 END_CALCS

//  Index values
  init_ivector nyrs_srv(1,nsrv)                   //Number of years of survey index value (annual)
  !! log_input(nyrs_srv);
  init_imatrix yrs_srv(1,nsrv,1,nyrs_srv)         //Years of survey index value (annual)
  !! log_input(yrs_srv);
  init_vector mo_srv(1,nsrv)                      //Month surveys occur 
  !! log_input(mo_srv);
  vector srv_month_frac(1,nsrv)
  !! srv_month_frac = (mo_srv-1.)/12.;
  init_matrix obs_srv(1,nsrv,1,nyrs_srv)          //values of survey index value (annual)
  !! log_input(obs_srv);
  init_matrix obs_se_srv(1,nsrv,1,nyrs_srv)       //values of survey serrs
  !! log_input(obs_se_srv);
  matrix     obs_lse_srv(1,nsrv,1,nyrs_srv) //Survey standard errors (for lognormal)
  matrix     obs_lva_srv(1,nsrv,1,nyrs_srv) //Survey standard errors (for lognormal)
  matrix        corr_dev(1,nsrv,1,nyrs_srv) //Survey standard errors (for lognormal)
  matrix        corr_eff(1,nfsh,styr,endyr) //Survey standard errors (for lognormal)
  matrix         act_eff(1,nfsh,styr,endyr) //Survey standard errors (for lognormal)
  vector              ac(1,nsrv);
  !! obs_lse_srv = elem_div(obs_se_srv,obs_srv);
  !! obs_lse_srv = sqrt(log(square(obs_lse_srv) + 1.));
  !! obs_lva_srv = square(obs_lse_srv);

  init_ivector nyrs_srv_age(1,nsrv)               //Number of years of survey index value (annual)
  !! log_input(nyrs_srv_age);
  init_matrix yrs_srv_age(1,nsrv,1,nyrs_srv_age)  //Years of survey index value (annual)
  !! log_input(yrs_srv_age);
  init_matrix n_sample_srv_age(1,nsrv,1,nyrs_srv_age)         //Years of survey index value (annual)
  !! log_input(n_sample_srv_age);
  init_3darray oac_srv(1,nsrv,1,nyrs_srv_age,1,nages)  //values of Survey proportions at age
  !! log_input(oac_srv);
  init_3darray  wt_srv(1,nsrv,styr,endyr,1,nages)      //values of Survey proportions at age
  !! log_input(wt_srv);

  vector age_vector(1,nages);
    !! for (int j=1;j<=nages;j++)
      !!  age_vector(j) = double(j+rec_age-1);
  init_vector wt_pop(1,nages)
  !! log_input(wt_pop);
  init_vector maturity(1,nages)
  !! log_input(maturity);
  // !! if (max(maturity)>.9) maturity /=2.;
  vector wt_mature(1,nages);
  !! wt_mature = elem_prod(wt_pop,maturity) ;

  //Spawning month-----
  init_number spawnmo
  number spmo_frac
  !! spmo_frac = (spawnmo-1)/12.;

  init_matrix age_err(1,nages,1,nages)
  !! log_input(age_err);

  int k // Index for fishery or survey
  int i // Index for year
  int j // Index for age
 LOCAL_CALCS
  // Rename data file to the control data section... 
  ad_comm::change_datafile_name(cntrlfile_name);
  *(ad_comm::global_datafile) >>  datafile_name; 
  *(ad_comm::global_datafile) >>  model_name; 
  log_input(cntrlfile_name);
 END_CALCS
  // Matrix of selectivity mappings--row 1 is type (1=fishery, 2=index) and row 2 is index within that type
  //  e.g., the following for 2 fisheries and 4 indices means that index 3 uses fishery 1 selectivities,
  //         the other fisheries and indices use their own parameterization
  //  1 1 2 2 1 2 
  //  1 2 1 2 1 4
  init_imatrix sel_map(1,2,1,nfsh_and_srv) 
  // maps fisheries and surveys into sequential sel_map for sharing purposes
  !! write_input_log<< "# Map shared selectivity: "<< endl;log_input(sel_map);
  !! log_input(datafile_name);
  !! log_input(model_name);
  !! projfile_name = cntrlfile_name(1,length(cntrlfile_name)-4) + ".prj";
  init_int    SrType        // 2 Bholt, 1 Ricker
  init_int use_age_err      // nonzero value means use...
  init_int use_effort       // Effort time series (not implemented fully)
  init_number steepnessprior
  init_number cvsteepnessprior
  init_int    phase_srec

  init_number sigmarprior
  number log_sigmarprior
  init_number cvsigmarprior
  init_int    phase_sigmar
  init_int    styr_rec_est
  init_int    endyr_rec_est
  int nrecs_est;

  init_number natmortprior
  init_number cvnatmortprior
  init_int    phase_M

  init_vector qprior(1,nsrv)      
  vector log_qprior(1,nsrv)      
  init_vector cvqprior(1,nsrv)     
  init_ivector phase_q(1,nsrv)

  init_vector q_power_prior(1,nsrv)      
  vector log_q_power_prior(1,nsrv)      
  init_vector cvq_power_prior(1,nsrv)     
  init_ivector phase_q_power(1,nsrv)

  init_ivector    q_age_min(1,nsrv)     // Age that q relates to...
  init_ivector    q_age_max(1,nsrv)     // Age that q relates to...

  init_number cv_catchbiomass
  number catchbiomass_pen
  !!catchbiomass_pen= 1./(2*cv_catchbiomass*cv_catchbiomass);
  init_int nproj_yrs

  int styr_fut
  int endyr_fut            // LAst year for projections
  int phase_Rzero
  int phase_nosr
  number Steepness_UB
  !! phase_Rzero =  4;
  !! phase_nosr  = -3;
  !! styr_fut    = endyr+1;
  !! endyr_fut   = endyr + nproj_yrs; 

  // Selectivity controls
  // read in options for each fishery
  // Loop over fisheries and surveys to read in data (conditional on sel_options)
  ivector   fsh_sel_opt(1,nfsh)
  ivector phase_sel_fsh(1,nfsh)
  vector   curv_pen_fsh(1,nfsh)
  matrix   sel_slp_in_fsh(1,nfsh,1,nyrs)
  matrix   logsel_slp_in_fsh(1,nfsh,1,nyrs)
  matrix   sel_inf_in_fsh(1,nfsh,1,nyrs)
  vector   logsel_slp_in_fshv(1,nfsh)
  vector   sel_inf_in_fshv(1,nfsh)
  vector   logsel_dslp_in_fshv(1,nfsh)
  vector   sel_dinf_in_fshv(1,nfsh)
  matrix   sel_dslp_in_fsh(1,nfsh,1,nyrs)
  matrix   logsel_dslp_in_fsh(1,nfsh,1,nyrs)
  matrix   sel_dinf_in_fsh(1,nfsh,1,nyrs)
  !! logsel_slp_in_fshv.initialize();
  !! sel_inf_in_fshv.initialize();
  !! logsel_dslp_in_fshv.initialize();
  !! sel_inf_in_fshv.initialize();
  !! sel_dinf_in_fshv.initialize();

  !! logsel_slp_in_srvv.initialize();
  !! sel_inf_in_srvv.initialize();
  !! logsel_dslp_in_srvv.initialize();
  !! sel_inf_in_srvv.initialize();
  !! sel_dinf_in_srvv.initialize();

  vector seldec_pen_fsh(1,nfsh) ;
  vector nnodes_fsh(1,nfsh) ;
  int seldecage ;
  !! seldecage = int(nages/2);
  matrix sel_change_in_fsh(1,nfsh,styr,endyr);
  ivector nselages_in_fsh(1,nfsh)

  ivector n_sel_ch_fsh(1,nfsh);
  ivector n_sel_ch_srv(1,nsrv);
  imatrix yrs_sel_ch_tmp(1,20,1,endyr-styr+1);
  imatrix yrs_sel_ch_tmp_srv(1,nsrv,1,endyr-styr+1);
  !! yrs_sel_ch_tmp.initialize();
  !! yrs_sel_ch_tmp_srv.initialize();

  ivector   srv_sel_opt(1,nsrv)
  ivector phase_sel_srv(1,nsrv)

  vector   curv_pen_srv(1,nsrv)

  matrix   logsel_slp_in_srv(1,nsrv,1,nyrs)
  matrix   sel_inf_in_srv(1,nsrv,1,nyrs)
  matrix   sel_dslp_in_srv(1,nsrv,1,nyrs)
  matrix   logsel_dslp_in_srv(1,nsrv,1,nyrs)
  matrix   sel_dinf_in_srv(1,nsrv,1,nyrs)
  matrix   sel_slp_in_srv(1,nsrv,1,nyrs)

  vector   logsel_slp_in_srvv(1,nsrv)
  vector   sel_inf_in_srvv(1,nsrv)
  vector   logsel_dslp_in_srvv(1,nsrv)
  vector   sel_dinf_in_srvv(1,nsrv)

  vector seldec_pen_srv(1,nsrv) ;
  matrix sel_change_in_srv(1,nsrv,styr,endyr);
  ivector nselages_in_srv(1,nsrv)

  // Phase of estimation
  ivector phase_selcoff_fsh(1,nfsh)
  ivector phase_logist_fsh(1,nfsh)
  ivector phase_dlogist_fsh(1,nfsh)
  ivector phase_sel_spl_fsh(1,nfsh)

  ivector phase_selcoff_srv(1,nsrv)
  ivector phase_logist_srv(1,nsrv)
  ivector phase_dlogist_srv(1,nsrv)
  vector  sel_fsh_tmp(1,nages); 
  vector  sel_srv_tmp(1,nages); 
  3darray log_selcoffs_fsh_in(1,nfsh,1,nyrs,1,nages)
  3darray log_selcoffs_srv_in(1,nsrv,1,nyrs,1,nages)
  3darray  log_sel_spl_fsh_in(1,nfsh,1,nyrs,1,nages) // use nages for input to start
  // 3darray log_selcoffs_srv_in(1,nsrv,1,nyrs,1,nages)

 LOCAL_CALCS
  phase_selcoff_srv.initialize();
  phase_logist_srv.initialize();
  phase_dlogist_srv.initialize();
  sel_fsh_tmp.initialize() ;
  sel_srv_tmp.initialize() ;
  log_selcoffs_fsh_in.initialize();
  log_selcoffs_srv_in.initialize();

  // nselages_in_fsh.initialize()   ;  
  // nselages_in_srv.initialize()   ;  
  nselages_in_fsh = nages-1;
  nselages_in_srv = nages-1;
  sel_change_in_fsh.initialize()   ;  
  sel_change_in_srv.initialize()   ;  
  sel_slp_in_fsh.initialize()   ;  // ji
  sel_slp_in_srv.initialize()   ;  // ji
  sel_inf_in_fsh.initialize()   ;  // ji
  sel_inf_in_srv.initialize()   ;  // ji
  logsel_slp_in_fsh.initialize();  // ji
  logsel_slp_in_fshv.initialize();  // ji
  logsel_dslp_in_fsh.initialize(); // ji
  logsel_slp_in_srv.initialize();  // ji
  logsel_slp_in_srvv.initialize();  // ji
  logsel_dslp_in_srv.initialize(); // ji
  for (k=1;k<=nfsh;k++)
  {
    *(ad_comm::global_datafile) >> fsh_sel_opt(k)  ;  
    write_input_log<<fshname(k)<<" fishery Sel Option: "<<fsh_sel_opt(k)<<endl;
    switch (fsh_sel_opt(k))
    {
      case 1 : // Selectivity coefficients 
      {
        *(ad_comm::global_datafile) >> nselages_in_fsh(k)   ;  
        *(ad_comm::global_datafile) >> phase_sel_fsh(k);  
        *(ad_comm::global_datafile) >> curv_pen_fsh(k) ;
        *(ad_comm::global_datafile) >> seldec_pen_fsh(k) ;
        seldec_pen_fsh(k) *= seldec_pen_fsh(k) ;
        for (int i=styr;i<=endyr;i++)
          *(ad_comm::global_datafile) >> sel_change_in_fsh(k,i) ;

        sel_change_in_fsh(k,styr)=1.; 
        n_sel_ch_fsh(k)=0.;
        int j=1; 
        yrs_sel_ch_tmp(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_fsh(k,i)>0) {
            j++; 
            yrs_sel_ch_tmp(k,j) = i; } }

        n_sel_ch_fsh(k) = j; 

        // This to read in pre-specified selectivity values...
        sel_fsh_tmp.initialize();
        log_selcoffs_fsh_in.initialize();
        for (j=1;j<=nages;j++) 
            *(ad_comm::global_datafile) >> sel_fsh_tmp(j);  
        for (int jj=2;jj<=n_sel_ch_fsh(k);jj++) 
        {
          // Set the selectivity for the oldest group
          for (j=nselages_in_fsh(k)+1;j<=nages;j++) 
          {
            sel_fsh_tmp(j)  = sel_fsh_tmp(nselages_in_fsh(k));  
          }
          // Set tmp to actual initial vectors...
          log_selcoffs_fsh_in(k,jj)(1,nselages_in_fsh(k)) = log(sel_fsh_tmp(1,nselages_in_fsh(k))+1e-7/mean(sel_fsh_tmp(1,nselages_in_fsh(k))+1e-7) );
          write_input_log<<"Sel_in_fsh "<< mfexp(log_selcoffs_fsh_in(k,jj))<<endl;
        }
        // exit(1);
        phase_selcoff_fsh(k) = phase_sel_fsh(k);
        phase_logist_fsh(k)  = -1;
        phase_dlogist_fsh(k) = -1;
        phase_sel_spl_fsh(k) = -1;
      }
        break;
      case 2 : // Single logistic
      {
        *(ad_comm::global_datafile) >> phase_sel_fsh(k);  
        for (int i=styr;i<=endyr;i++) { 
          *(ad_comm::global_datafile) >> sel_change_in_fsh(k,i) ; }

        sel_change_in_fsh(k,styr)=1.; 
        n_sel_ch_fsh(k)=0;
        int j=1; 
        yrs_sel_ch_tmp(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_fsh(k,i)>0) {
            j++; yrs_sel_ch_tmp(k,j) = i; } }
        n_sel_ch_fsh(k) = j; 

        // This to read in pre-specified selectivity values...
        *(ad_comm::global_datafile) >> sel_slp_in_fsh(k,1) ;
        *(ad_comm::global_datafile) >> sel_inf_in_fsh(k,1) ;
        for (int jj=2;jj<=n_sel_ch_fsh(k);jj++) 
        {
          sel_inf_in_fsh(k,jj)    =     sel_inf_in_fsh(k,1) ;
          logsel_slp_in_fsh(k,jj) = log(sel_slp_in_fsh(k,1)) ;
        }

        phase_selcoff_fsh(k) = -1;
        phase_logist_fsh(k) = phase_sel_fsh(k);
        phase_dlogist_fsh(k) = -1;

        logsel_slp_in_fshv(k) = logsel_slp_in_fsh(k,1);
           sel_inf_in_fshv(k) =    sel_inf_in_fsh(k,1);
        break;
      }
      case 3 : // Double logistic 
      {
        write_input_log << "Double logistic abandoned..."<<endl;exit(1);
        *(ad_comm::global_datafile) >> phase_sel_fsh(k);  
        for (int i=styr;i<=endyr;i++) { 
          *(ad_comm::global_datafile) >> sel_change_in_fsh(k,i) ; }
        sel_change_in_fsh(k,styr)=1.; 
        n_sel_ch_fsh(k)=0;
        int j=1; 
        yrs_sel_ch_tmp(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_fsh(k,i)>0) {
            j++; yrs_sel_ch_tmp(k,j) = i; } }
        n_sel_ch_fsh(k) = j; 

        // This to read in pre-specified selectivity values...
        for (int jj=1;jj<=n_sel_ch_fsh(k);jj++) 
        {
          *(ad_comm::global_datafile) >> sel_slp_in_fsh(k,jj) ;
          *(ad_comm::global_datafile) >> sel_inf_in_fsh(k,jj) ;
          *(ad_comm::global_datafile) >> sel_dslp_in_fsh(k,jj) ;
          *(ad_comm::global_datafile) >> sel_dinf_in_fsh(k,jj) ;
          logsel_slp_in_fsh(k,jj) = log(sel_slp_in_fsh(k,jj)) ;
          logsel_dslp_in_fsh(k,jj) = log(sel_dslp_in_fsh(k,jj)) ;
        }


        phase_selcoff_fsh(k) = -1;
        phase_logist_fsh(k)  = phase_sel_fsh(k);
        phase_dlogist_fsh(k) = phase_sel_fsh(k)+1;


        logsel_slp_in_fshv(k) = logsel_slp_in_fsh(k,1);
           sel_inf_in_fshv(k) =    sel_inf_in_fsh(k,1);
        logsel_dslp_in_fshv(k) = logsel_dslp_in_fsh(k,1);
           sel_dinf_in_fshv(k) =    sel_dinf_in_fsh(k,1);
        break;
      }
      case 4 : // Splines         
      {
        write_input_log << "Spline selectivity option for fishery "<<fshname(k)<<endl;
        *(ad_comm::global_datafile) >> nselages_in_fsh(k)   ; // Number of selected ages 
        *(ad_comm::global_datafile) >> phase_sel_fsh(k);       
        *(ad_comm::global_datafile) >> nnodes_fsh(k)   ;      // Number of nodes, to run over selected ages only
        *(ad_comm::global_datafile) >> seldec_pen_fsh(k) ;
        seldec_pen_fsh(k) *= seldec_pen_fsh(k) ;
        // Read in year-year changes in selectivity
        for (int i=styr;i<=endyr;i++)
          *(ad_comm::global_datafile) >> sel_change_in_fsh(k,i) ;

        sel_change_in_fsh(k,styr)=1.; 
        n_sel_ch_fsh(k)=0;
        int j=1; 
        yrs_sel_ch_tmp(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_fsh(k,i)>0) {
            j++; 
            yrs_sel_ch_tmp(k,j) = i; } }
        n_sel_ch_fsh(k) = j; 

        log_input(nselages_in_fsh);
        log_input(nnodes_fsh);
        log_input(seldec_pen_fsh);
        log_input(n_sel_ch_fsh);
        log_input(sel_change_in_fsh);
        log_input(yrs_sel_ch_tmp(1,n_sel_ch_fsh(k)));

        // This to read in pre-specified selectivity values...
        sel_fsh_tmp.initialize();
        log_sel_spl_fsh_in.initialize();
        for (j=1;j<=nnodes_fsh(k);j++) 
            *(ad_comm::global_datafile) >> sel_fsh_tmp(j);  

        for (int jj=2;jj<=n_sel_ch_fsh(k);jj++) 
        {
          log_sel_spl_fsh_in(k,jj)(1,nnodes_fsh(k)) = log(sel_fsh_tmp(1,nnodes_fsh(k)));
          write_input_log<<"Initial spline nodal values "<< log_sel_spl_fsh_in(k,jj)(1,nnodes_fsh(k))<<endl;
        }
        // exit(1);
        phase_sel_spl_fsh(k) = phase_sel_fsh(k);
        log_input(phase_sel_spl_fsh);
        phase_selcoff_fsh(k) = -1;
        phase_logist_fsh(k)  = -1;
        phase_dlogist_fsh(k) = -1;
      }
      break;
      write_input_log << fshname(k)<<" fish sel opt "<<endl<<fsh_sel_opt(k)<<" "<<endl<<"Sel_change"<<endl<<sel_change_in_fsh(k)<<endl;
    }
  }
  // Surveys here..............
  for(k=1;k<=nsrv;k++)
  {
    *(ad_comm::global_datafile) >> srv_sel_opt(k)  ;  
    write_input_log << endl<<"Survey "<<srvname(k)<<" Selectivity option: "<<srv_sel_opt(k)  <<endl;  
    switch (srv_sel_opt(k))
    {
      case 1 : // Selectivity coefficients  surveys
      {
        *(ad_comm::global_datafile) >> nselages_in_srv(k)   ;  
        *(ad_comm::global_datafile) >> phase_sel_srv(k);  
        *(ad_comm::global_datafile) >> curv_pen_srv(k) ;
        *(ad_comm::global_datafile) >> seldec_pen_srv(k) ;
        seldec_pen_srv(k) *= seldec_pen_srv(k) ;
        write_input_log << "Survey: "          <<srvname(k)<<endl; // cin>>junk;
        write_input_log << "Survey selages: "  <<nselages_in_srv(k)<<endl; // cin>>junk;
        write_input_log << "Survey sel_phase: "<< phase_sel_srv(k)<<endl; // cin>>junk;
        write_input_log << "srv seldec_penalt: "<<seldec_pen_srv(k) << endl;
        for (int i=styr;i<=endyr;i++)
          *(ad_comm::global_datafile) >> sel_change_in_srv(k,i) ;

        sel_change_in_srv(k,styr)=1.; 
        n_sel_ch_srv(k)=0; 
        int j=1; 
        yrs_sel_ch_tmp_srv(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_srv(k,i)>0) {
            j++; yrs_sel_ch_tmp_srv(k,j) = i; } }
        n_sel_ch_srv(k) = j;
        write_input_log << "srv sel_change_in: "<<sel_change_in_srv(k) << endl;
        write_input_log <<"Number of select changes: survey "<<k<<": "<<n_sel_ch_srv(k)<<endl;
        write_input_log<<"Year sel change 1 "<<yrs_sel_ch_tmp_srv(1,n_sel_ch_srv(k))<<endl;// exit(1);

        // This to read in pre-specified selectivity values...
        sel_srv_tmp.initialize();
        log_selcoffs_srv_in.initialize();
        for (j=1;j<=nages;j++) 
            *(ad_comm::global_datafile) >> sel_srv_tmp(j);  
        for (int jj=2;jj<=n_sel_ch_srv(k);jj++) 
        {
          for (j=nselages_in_srv(k)+1;j<=nages;j++) 
          {
            sel_srv_tmp(j)  = sel_srv_tmp(nselages_in_srv(k));  
          }
          // Set tmp to actual initial vectors...
          log_selcoffs_srv_in(k,jj)(1,nselages_in_srv(k)) = log(sel_srv_tmp(1,nselages_in_srv(k))+1e-7/mean(sel_fsh_tmp(1,nselages_in_fsh(k))+1e-7) );
          write_input_log<<"Sel_in_srv "<< mfexp(log_selcoffs_srv_in(k,jj))<<endl;
        }
        phase_selcoff_srv(k) = phase_sel_srv(k);
        phase_logist_srv(k)  = -1;
        phase_dlogist_srv(k)  = -1;
      }

      break;

      case 2 : // Single logistic
      {
        *(ad_comm::global_datafile) >> phase_sel_srv(k);  
        // cout<<"PHase sel srv "<<phase_sel_srv(k)<<endl;

        for (int i=styr;i<=endyr;i++) { 
          *(ad_comm::global_datafile) >> sel_change_in_srv(k,i) ; }

        sel_change_in_srv(k,styr)=1.; n_sel_ch_srv(k)=0; int j=1; 
        write_input_log << "srv sel_change_in: "<<sel_change_in_srv(k) << endl;
        yrs_sel_ch_tmp_srv(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) { if(sel_change_in_srv(k,i)>0) { j++; yrs_sel_ch_tmp_srv(k,j) = i; } }
        n_sel_ch_srv(k) = j; write_input_log <<"Number of select changes: survey "<<k<<": "<<n_sel_ch_srv(k)<<endl;
        write_input_log<<"Year sel change 1 "<<yrs_sel_ch_tmp_srv(1,n_sel_ch_srv(k))<<endl;// exit(1);

        // This to read in pre-specified selectivity values...
        *(ad_comm::global_datafile) >> sel_slp_in_srv(k,1) ;
        *(ad_comm::global_datafile) >> sel_inf_in_srv(k,1) ;
        for (int jj=2;jj<=n_sel_ch_srv(k);jj++) {
          sel_inf_in_srv(k,jj)    =     sel_inf_in_srv(k,1) ;
          logsel_slp_in_srv(k,jj) = log(sel_slp_in_srv(k,1)) ;
        }

        write_input_log<<k<<" "<< phase_sel_srv(k)<<endl;
        write_input_log<<k<<" "<< sel_slp_in_srv(k)(1,n_sel_ch_srv(k)) <<endl;// exit(1);
        phase_selcoff_srv(k) = -1;
        phase_logist_srv(k) = phase_sel_srv(k);
        phase_dlogist_srv(k)  = -1;

        logsel_slp_in_srvv(k) = logsel_slp_in_srv(k,1);
           sel_inf_in_srvv(k) =    sel_inf_in_srv(k,1);
      }
      break;
      case 3 : // Double logistic 
      {
        write_input_log << "Double logistic abandoned..."<<endl;exit(1);
        *(ad_comm::global_datafile) >> phase_sel_srv(k);  
        for (int i=styr;i<=endyr;i++) { 
          *(ad_comm::global_datafile) >> sel_change_in_srv(k,i) ; }
        sel_change_in_srv(k,styr)=1.; n_sel_ch_srv(k)=0; int j=1; 
        write_input_log << "srv sel_change_in: "<<sel_change_in_srv(k) << endl;
        yrs_sel_ch_tmp_srv(k,j) = styr;
       // Number of selectivity changes is equal to the number of vectors (yr 1 is baseline)
        for (int i=styr+1;i<=endyr;i++) {
          if(sel_change_in_srv(k,i)>0) {
            j++; yrs_sel_ch_tmp_srv(k,j) = i; } }
        n_sel_ch_srv(k) = j; write_input_log <<"Number of select changes: survey "<<k<<": "<<n_sel_ch_srv(k)<<endl;
        write_input_log<<"Year sel change 1 "<<yrs_sel_ch_tmp_srv(1,n_sel_ch_srv(k))<<endl;// exit(1);

        // This to read in pre-specified selectivity values...
        for (int jj=1;jj<=n_sel_ch_srv(k);jj++) {
          *(ad_comm::global_datafile) >> sel_slp_in_srv(k,jj) ;
          *(ad_comm::global_datafile) >> sel_inf_in_srv(k,jj) ;
          *(ad_comm::global_datafile) >> sel_dslp_in_srv(k,jj) ;
          *(ad_comm::global_datafile) >> sel_dinf_in_srv(k,jj) ;
          logsel_slp_in_srv(k,jj) = log(sel_slp_in_srv(k,jj)) ;
          logsel_dslp_in_srv(k,jj) = log(sel_dslp_in_srv(k,jj)) ;
        }


        phase_selcoff_srv(k) = -1;
        phase_logist_srv(k)  = phase_sel_srv(k);
        phase_dlogist_srv(k) = phase_sel_srv(k)+1;

        logsel_slp_in_srvv(k) = logsel_slp_in_srv(k,1);
           sel_inf_in_srvv(k) =    sel_inf_in_srv(k,1);
        logsel_dslp_in_srvv(k) = logsel_dslp_in_srv(k,1);
           sel_dinf_in_srvv(k) =    sel_dinf_in_srv(k,1);
      }
        break;
      case 4 : // spline for surveys
      {
      }
      break;
    }
    write_input_log << srvname(k)<<" srv sel opt "<<srv_sel_opt(k)<<" "<<sel_change_in_srv(k)<<endl;
  }
  write_input_log<<"Phase survey Sel_Coffs: "<<phase_selcoff_srv<<endl; 
 END_CALCS
  init_number test;
  !! write_input_log<<" Test: "<<test<<endl;
 !! if (test!=123456789) {cerr<<"Control file not read in correctly... "<<endl;exit(1);}


  ivector nopt_fsh(1,2) // number of options...
  !! nopt_fsh.initialize();
  !! for (k=1;k<=nfsh;k++) if(fsh_sel_opt(k)==1) nopt_fsh(1)++;else nopt_fsh(2)++;

  // Fishery selectivity description:
  // type 1
  
  // Number of ages

  !! write_input_log << "# Fshry Selages: " << nselages_in_fsh  <<endl;
  !! write_input_log << "# Srvy  Selages: " << nselages_in_srv <<endl;



  !! write_input_log << "# Phase for age-spec fishery "<<phase_selcoff_fsh<<endl;
  !! write_input_log << "# Phase for logistic fishery "<<phase_logist_fsh<<endl;
  !! write_input_log << "# Phase for dble logistic fishery "<<phase_dlogist_fsh<<endl;

  !! write_input_log << "# Phase for age-spec survey  "<<phase_selcoff_srv<<endl;
  !! write_input_log << "# Phase for logistic survey  "<<phase_logist_srv<<endl;
  !! write_input_log << "# Phase for dble logistic srvy "<<phase_dlogist_srv<<endl;

  !! for (k=1;k<=nfsh;k++) if (phase_selcoff_fsh(k)>0) curv_pen_fsh(k) = 1./ (square(curv_pen_fsh(k))*2);
  !! write_input_log<<"# Curv_pen_fsh: "<<endl<<curv_pen_fsh<<endl;
  !! for (k=1;k<=nsrv;k++) if (phase_selcoff_srv(k)>0) curv_pen_srv(k) = 1./ (square(curv_pen_srv(k))*2);
  !! write_input_log<<"# Curv_pen_srv: "<<endl<<curv_pen_fsh<<endl;

  int  phase_fmort;
  int  phase_proj;
  imatrix yrs_sel_ch_fsh(1,nfsh,1,n_sel_ch_fsh);
  imatrix   nselages_fsh(1,nfsh,1,n_sel_ch_fsh);
  matrix xnodes_fsh(1,nfsh,1,nnodes_fsh)
  matrix xages_fsh(1,nfsh,1,nages)

  imatrix yrs_sel_ch_srv(1,nsrv,1,n_sel_ch_srv);
  imatrix   nselages_srv(1,nsrv,1,n_sel_ch_srv);
 LOCAL_CALCS
  for (k=1; k<=nfsh;k++)
  {
    // xages_fsh increments from 0-1 by number of ages, say
    xages_fsh.initialize();
    log_input(xages_fsh);
    xages_fsh(k).fill_seqadd(0.,1.0/(nages-1));
    log_input(xages_fsh);
    //  xnodes increments from 0-1 by number of nodes
    xnodes_fsh.initialize();
    xnodes_fsh(k).fill_seqadd(0.,1.0/(nnodes_fsh(k)-1));
    log_input(xnodes_fsh);
    // xages_fsh(k).fill_seqadd(0,1.0/(nselages_in_fsh(k)-1)); //prefer to use nselages but need 3d version to work
    yrs_sel_ch_fsh(k) = yrs_sel_ch_tmp(k)(1,n_sel_ch_fsh(k));
  }
  write_input_log<<"Yrs fsh_sel change: "<<yrs_sel_ch_fsh<<endl;
  for (k=1; k<=nsrv;k++)
    yrs_sel_ch_srv(k) = yrs_sel_ch_tmp_srv(k)(1,n_sel_ch_srv(k));
  write_input_log<<"Yrs srv_sel change: "<<yrs_sel_ch_srv<<endl;
    log_sigmarprior = log(sigmarprior);
    log_input(steepnessprior);
    log_input(sigmarprior);
    nrecs_est = endyr_rec_est-styr_rec_est+1;
    nrecs_est = endyr_rec_est-styr_rec_est+1;
    write_input_log<<"#  SSB estimated in styr endyr: " <<styr_sp    <<" "<<endyr_sp      <<" "<<endl;
    write_input_log<<"#  Rec estimated in styr endyr: " <<styr_rec    <<" "<<endyr        <<" "<<endl;
    write_input_log<<"#  SR Curve fit  in styr endyr: " <<styr_rec_est<<" "<<endyr_rec_est<<" "<<endl;
    write_input_log<<"#             Model styr endyr: " <<styr        <<" "<<endyr        <<" "<<endl;
    log_qprior = log(qprior);
    log_input(qprior);
    log_q_power_prior = log(q_power_prior);
    write_input_log<<"# q_power_prior " <<endl<<q_power_prior<<" "<<endl;
    write_input_log<<"# cv_catchbiomass " <<endl<<cv_catchbiomass<<" "<<endl;
    write_input_log<<"# CatchbiomassPen " <<endl<<catchbiomass_pen<<" "<<endl;
    write_input_log<<"# Number of projection years " <<endl<<nproj_yrs<<" "<<endl;// cin>>junk;

 END_CALCS
  number R_guess;

  vector offset_srv(1,nsrv)
  vector offset_fsh(1,nfsh)
  int do_fmort;
  !! do_fmort=0;
  int Popes;
 LOCAL_CALCS
  Popes=0; // option to do Pope's approximation (not presently flagged outside of code)
  if (Popes) 
    phase_fmort = -2;
  else
    phase_fmort = 1;

  phase_proj  =  5;

  Steepness_UB = .9999; // upper bound of steepness
  offset_srv.initialize();
  offset_fsh.initialize();

  double sumtmp;
  for (k=1;k<=nfsh;k++)
    for (i=1;i<=nyrs_fsh_age(k);i++)
    {
      oac_fsh(k,i) /= sum(oac_fsh(k,i)); // Normalize to sum to one
      offset_fsh(k) -= n_sample_fsh_age(k,i)*(oac_fsh(k,i) + 0.001) * log(oac_fsh(k,i) + 0.001 ) ;
    }

  for (k=1;k<=nsrv;k++)
    for (i=1;i<=nyrs_srv_age(k);i++)
    {
      oac_srv(k,i) /= sum(oac_srv(k,i)); // Normalize to sum to one
      offset_srv(k) -= n_sample_srv_age(k,i)*(oac_srv(k,i) + 0.001) * log(oac_srv(k,i) + 0.001 ) ;
    }
  write_input_log << "# Multinomial Offset surveys: "<<endl;
  write_input_log<<offset_srv<<endl;
  write_input_log << "# Multinomial Offset fisheries: "<<endl;
  write_input_log<<offset_fsh<<endl;

  if (ad_comm::argc > 1) // Command line argument to profile Fishing mortality rates...
  {
    int on=0;
    if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-uFmort"))>-1)
      do_fmort=1;
  }

  // Compute an initial Rzero value based on exploitation 
   double btmp=0.;
   double ctmp=0.;
   dvector ntmp(1,nages);
   ntmp(1) = 1.;
   for (int a=2;a<=nages;a++)
     ntmp(a) = ntmp(a-1)*exp(-natmortprior-.05);
   btmp = wt_pop * ntmp;
   write_input_log << "Mean Catch"<<endl;
   ctmp = mean(catch_bio);
   write_input_log << ctmp <<endl;
   R_guess = log((ctmp/.02 )/btmp) ;
   write_input_log << "R_guess "<<endl;
   write_input_log << R_guess <<endl;
 END_CALCS

PARAMETER_SECTION
 // Biological Parameters
  init_bounded_number M(.02,.8,phase_M)
  matrix  natage(styr,endyr+1,1,nages)
  // vector Sp_Biom(styr_sp,endyr)
  vector pred_rec(styr_rec,endyr)
  vector mod_rec(styr_rec,endyr) // As estimated by model
  matrix  Z(styr,endyr,1,nages)
  matrix  S(styr,endyr,1,nages)
  number  surv
  number  natmort
 
 // Stock rectuitment params
  init_number mean_log_rec(1); 
  init_bounded_number steepness(0.21,Steepness_UB,phase_srec)
  init_number log_Rzero(phase_Rzero)  
  // OjO
  // init_bounded_vector initage_dev(2,nages,-15,15,4)
  init_bounded_vector rec_dev(styr_rec,endyr,-15,15,2)
  // init_vector rec_dev(styr_rec,endyr,2)
  init_number log_sigmar(phase_sigmar);
  number m_sigmarsq  
  number m_sigmar
  number sigmarsq  
  number sigmar
  number alpha   
  number beta   
  number Bzero   
  number Rzero   
  number phizero
  number avg_rec_dev   

 // Fishing mortality parameters
  // init_vector         log_avg_fmort(1,nfsh,phase_fmort)
  // init_bounded_matrix fmort_dev(1,nfsh,styr,endyr,-15,15.,phase_fmort)
  init_bounded_matrix fmort(1,nfsh,styr,endyr,0.00,5.,phase_fmort)
  vector Fmort(styr,endyr);  // Annual total Fmort
  number hrate
  number catch_tmp
  number Fnew 

  !! for (k=1;k<=nfsh;k++) nselages_fsh(k)=nselages_in_fsh(k); // Sets all elements of a vector to one scalar value...
  !! for (k=1;k<=nsrv;k++) nselages_srv(k)=nselages_in_srv(k); // Sets all elements of a vector to one scalar value...

 //  init_3darray log_selcoffs_fsh(1,nfsh,1,n_sel_ch_fsh,1,nselages_fsh,phase_selcoff_fsh)
  init_matrix_vector log_selcoffs_fsh(1,nfsh,1,n_sel_ch_fsh,1,nselages_fsh,phase_selcoff_fsh) // 3rd dimension out...
  !! if (fsh_sel_opt(1)==4) nnodes_tmp=nnodes_fsh(1);  // NOTE THIS won't work in general
  init_matrix_vector  log_sel_spl_fsh(1,nfsh,1,n_sel_ch_fsh,1,nnodes_tmp,phase_sel_spl_fsh)

  !! log_input(nfsh);
  !! log_input(n_sel_ch_fsh);
  !! log_input(nselages_fsh);
  !! log_input(phase_selcoff_fsh);
  init_vector_vector logsel_slope_fsh(1,nfsh,1,n_sel_ch_fsh,phase_logist_fsh)
  matrix                sel_slope_fsh(1,nfsh,1,n_sel_ch_fsh)
  init_vector_vector     sel50_fsh(1,nfsh,1,n_sel_ch_fsh,phase_logist_fsh)
  init_vector_vector logsel_dslope_fsh(1,nfsh,1,n_sel_ch_fsh,phase_dlogist_fsh)
  matrix                sel_dslope_fsh(1,nfsh,1,n_sel_ch_fsh)
  !! int lb_d50=nages/2;
  init_bounded_vector_vector     seld50_fsh(1,nfsh,1,n_sel_ch_fsh,lb_d50,nages,phase_dlogist_fsh)

  // !!exit(1);
  3darray log_sel_fsh(1,nfsh,styr,endyr,1,nages)
  3darray sel_fsh(1,nfsh,styr,endyr,1,nages)
  matrix avgsel_fsh(1,nfsh,1,n_sel_ch_fsh);

  matrix  Ftot(styr,endyr,1,nages)
  3darray F(1,nfsh,styr,endyr,1,nages)
  3darray eac_fsh(1,nfsh,1,nyrs_fsh_age,1,nages)
  matrix  pred_catch(1,nfsh,styr,endyr)
  3darray catage(1,nfsh,styr,endyr,1,nages)
  matrix catage_tot(styr,endyr,1,nages)
  matrix expl_biom(1,nfsh,styr,endyr)

 // Parameters for computing SPR rates 
  vector F50(1,nfsh)
  vector F40(1,nfsh)
  vector F35(1,nfsh)

 // Stuff for SPR and yield projections
  number sigmar_fut
  vector f_tmp(1,nfsh)
  number SB0
  number SBF50
  number SBF40
  number SBF35
  vector Fratio(1,nfsh)
  !! Fratio = 1;
  !! Fratio /= sum(Fratio);

  matrix Nspr(1,4,1,nages)
 
  matrix nage_future(styr_fut,endyr_fut,1,nages)

  init_vector rec_dev_future(styr_fut,endyr_fut,phase_proj);


 
  vector Sp_Biom_future(styr_fut-rec_age,endyr_fut);

  3darray F_future(1,nfsh,styr_fut,endyr_fut,1,nages);
  matrix Z_future(styr_fut,endyr_fut,1,nages);
  matrix S_future(styr_fut,endyr_fut,1,nages);
  matrix catage_future(styr_fut,endyr_fut,1,nages);
  number avg_rec_dev_future
  vector avg_F_future(1,5)

 // Survey Observation parameters
  init_number_vector log_q_srv(1,nsrv,phase_q) 
  init_number_vector log_q_power_srv(1,nsrv,phase_q_power) 
  !! write_input_log <<nsrv<<endl<<n_sel_ch_srv<<endl<<nselages_srv<<endl<<phase_selcoff_srv<<endl;
  init_matrix_vector log_selcoffs_srv(1,nsrv,1,n_sel_ch_srv,1,nselages_srv,phase_selcoff_srv)
  init_vector_vector logsel_slope_srv(1,nsrv,1,n_sel_ch_srv,phase_logist_srv) // Need to make positive or reparameterize
  init_vector_vector logsel_dslope_srv(1,nsrv,1,n_sel_ch_srv,phase_dlogist_srv) // Need to make positive or reparameterize
  matrix                sel_slope_srv(1,nsrv,1,n_sel_ch_srv)
  matrix                sel_dslope_srv(1,nsrv,1,n_sel_ch_srv)
  init_vector_vector     sel50_srv(1,nsrv,1,n_sel_ch_srv,phase_logist_srv)
  init_bounded_vector_vector     seld50_srv(1,nfsh,1,n_sel_ch_srv,lb_d50,nages,phase_dlogist_srv)

  3darray log_sel_srv(1,nsrv,styr,endyr,1,nages)
  3darray sel_srv(1,nsrv,styr,endyr,1,nages)
  matrix avgsel_srv(1,nsrv,1,n_sel_ch_srv);

  matrix pred_srv(1,nsrv,styr,endyr)
  3darray eac_srv(1,nsrv,1,nyrs_srv_age,1,nages)

 // Likelihood value names         
  number sigma
  vector rec_like(1,4)
  vector catch_like(1,nfsh)
  vector age_like_fsh(1,nfsh)
  vector age_like_srv(1,nsrv)
  matrix sel_like_fsh(1,nfsh,1,4)       
  matrix sel_like_srv(1,nsrv,1,4)       
  vector surv_like(1,nsrv)
  vector fpen(1,6)    
  vector post_priors(1,4)
  vector post_priors_srvq(1,nsrv)
  objective_function_value obj_fun
  vector obj_comps(1,12)

  sdreport_number F50_est
  sdreport_number F40_est
  sdreport_number F35_est
  vector q_srv(1,nsrv)
  vector q_power_srv(1,nsrv)
  // sdreport_vector q_srv(1,nsrv)
  sdreport_vector totbiom(styr,endyr+1)
  sdreport_vector totbiom_NoFish(styr,endyr)
  sdreport_vector Sp_Biom(styr_sp,endyr+1)
  sdreport_vector Sp_Biom_NoFishR(styr,endyr)
  sdreport_number ABCBiom;
  sdreport_vector recruits(styr,endyr+1)
  // vector recruits(styr,endyr+1)
  sdreport_number depletion
  sdreport_number depletion_dyn
  sdreport_number MSY;
  sdreport_number MSYL;
  sdreport_number Fmsy;
  sdreport_number lnFmsy;
  sdreport_number Fcur_Fmsy;
  sdreport_number Rmsy;
  sdreport_number Bmsy;
  sdreport_number Bcur_Bmsy;
  // NOTE TO DAVE: Need to have a phase switch for sdreport variables(
  matrix catch_future(1,4,styr_fut,endyr_fut); // Note, don't project for F=0 (it will bomb)
  sdreport_matrix future_biomass(1,5,styr_fut,endyr_fut)
  !! write_input_log <<"logRzero "<<log_Rzero<<endl;
  !! write_input_log <<"logmeanrec "<<mean_log_rec<<endl;
  !! write_input_log<< "exp(log_sigmarprior "<<exp(log_sigmarprior)<<endl;
 // Initialize coefficients (if needed)
 LOCAL_CALCS
  for (k=1;k<=nfsh;k++) 
  {
    write_input_log<<"Fish sel phase: "<<phase_selcoff_fsh(k)<<" "<<fshname(k)<<endl;
    switch (fsh_sel_opt(k))
    {
      case 1 : // Selectivity coefficients 
      {
        if(phase_selcoff_fsh(k)<0)
        {
          write_input_log<<"Initial fixing fishery sel to"<<endl<<n_sel_ch_fsh(k)<<endl;
          for (int jj=1;jj<=n_sel_ch_fsh(k);jj++) 
          {
            log_selcoffs_fsh(k,jj)(1,nselages_in_fsh(k)) = log_selcoffs_fsh_in(k,jj)(1,nselages_in_fsh(k));
            write_input_log <<"Init coef:"<<endl<<exp(log_selcoffs_fsh(k,jj)(1,nselages_in_fsh(k))) <<endl;
          }
        }
      }
        break;
      case 2 : // Single logistic
      {
        if(phase_logist_fsh(k)<0)
        {
          write_input_log<<"Fixing fishery sel to"<<endl<<n_sel_ch_fsh(k)<<endl;
          for (int jj=1;jj<=n_sel_ch_fsh(k);jj++) 
          {
            logsel_slope_fsh(k,jj) = logsel_slp_in_fsh(k,jj)  ;
            sel50_fsh(k,jj)        =    sel_inf_in_fsh(k,jj)  ;
          }
        }
      }
      case 3 : // Double logistic 
      {
        if(phase_dlogist_fsh(k)<0)
        {
          write_input_log<<"Fixing fishery sel to"<<endl<<n_sel_ch_fsh(k)<<endl;
          for (int jj=1;jj<=n_sel_ch_fsh(k);jj++) 
          {
            logsel_slope_fsh(k,jj) = logsel_slp_in_fsh(k,jj)  ;
            sel50_fsh(k,jj)        =    sel_inf_in_fsh(k,jj)  ;
          }
        }
      }
      case 4 : // Selectivity spline initialize 
      /* {
        if(phase_sel_spl_fsh(k)<0)
        {
          write_input_log<<"Initial fishery spline to"<<endl<<n_sel_ch_fsh(k)<<endl;
          for (int jj=1;jj<=n_sel_ch_fsh(k);jj++) 
          {
            log_sel_spl_fsh(k,jj)(1,nnodes_tmp) = log_sel_spl_fsh_in(k,jj)(1,nnodes_tmp);
            // write_input_log <<"Init coef:"<<endl<<exp(log_sel_spl_fsh(k,jj)(1,nselages_in_fsh(k))) <<endl;
          }
          log_input(log_sel_spl_fsh);
        }
       }*/
     break;
    }
  }
  for (k=1;k<=nsrv;k++) 
  {
    write_input_log<<"Srvy sel phase: "<<phase_selcoff_srv(k)<<endl;
    if(phase_selcoff_srv(k)<0)
    {
      write_input_log<<"Fixing "<<srvname(k)<<" survey sel to"<<endl<<n_sel_ch_srv(k)<<endl;
    // exit(1);
      for (int jj=1;jj<=n_sel_ch_srv(k);jj++) 
      {
        log_selcoffs_srv(k,jj)(1,nselages_in_srv(k)) = log_selcoffs_srv_in(k,jj)(1,nselages_in_srv(k));
        // write_input_log <<"Init coef:"<<endl<<exp(log_selcoffs_srv(k,jj)(1,nselages_in_srv(k))) <<endl;
      }
    }
  }
  write_input_log <<"Leaving data secton"<<endl;
 END_CALCS

INITIALIZATION_SECTION
  M natmortprior; 
  steepness steepnessprior
  log_sigmar log_sigmarprior;
  log_Rzero R_guess;
  mean_log_rec R_guess;
  // log_avg_fmort -2.065
  log_q_srv log_qprior; 
  log_q_power_srv log_q_power_prior; 

  sel50_fsh sel_inf_in_fshv 

  logsel_dslope_fsh logsel_dslp_in_fshv ;
  seld50_fsh sel_dinf_in_fshv 

  logsel_slope_srv logsel_slp_in_srvv ;
  sel50_srv sel_inf_in_srvv ;

  logsel_dslope_srv logsel_dslp_in_srvv ;
  seld50_srv sel_dinf_in_srvv ;

 //+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==
PROCEDURE_SECTION
  natmort = M;
  fpen.initialize();
  for (k=1;k<=nsrv;k++) 
  {
    q_srv(k) = mfexp(log_q_srv(k) );
    q_power_srv(k) = mfexp(log_q_power_srv(k) );
  }

  // Main model calcs---------------------
  Get_Selectivity();
  Get_Mortality();
  Get_Bzero();
  Get_Numbers_at_Age();
  Get_Survey_Predictions();
  Get_Fishery_Predictions();

  // Output calcs-------------------------
  if (sd_phase())
  {
    compute_spr_rates();
    if (mcmcmode)
    {
      Calc_Dependent_Vars();
      mcflag   = 0;
      mcmcmode = 0;
    }
    else
    {
      if (mcflag)
        Calc_Dependent_Vars();
    }
  }
  // Objective function calcs------------
  evaluate_the_objective_function();
  // Other calcs-------------------------
  if (mceval_phase())
  {
    if (oper_mod)
      Oper_Model();
    else
    {
      compute_spr_rates();
      write_mceval();
    }
  }
  if (do_fmort)
    Profile_F();
 //+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==

FUNCTION write_mceval
  if (mcmcmode != 3)
    write_mceval_hdr();
  mcmcmode = 3;
  mceval<< model_name         << " "  ;
  mceval<< obj_fun            << " "  ;
  mceval<< sel_fsh(1,2007,12) << " "  ;
  mceval<<endl;

 /*
  get_msy();
  Future_projections();
  Calc_Dependent_Vars();
  mceval<<
  q_srv     << " "<< 
  M         << " "<< 
  steepness << " "<< 
  depletion << " "<< 
  MSY       << " "<< 
  MSYL      << " "<< 
  Fmsy      << " "<< 
  Fcur_Fmsy << " "<< 
  Bcur_Bmsy << " "<< 
  Bmsy      << " "<< 
  ABCBiom   << " "<< 
  F35       << " "<<
  F40       << " "<<
  F50       << " "<<
  future_biomass(1,endyr_fut) << " "<< 
  future_biomass(2,endyr_fut) << " "<< 
  future_biomass(3,endyr_fut) << " "<< 
  future_biomass(4,endyr_fut) << " "<< 
  future_biomass(5,endyr_fut) << " "<< 
  catch_future(1,styr_fut)    << " "<<  
  catch_future(2,styr_fut)    << " "<<  
  catch_future(3,styr_fut)    << " "<<  
  catch_future(4,styr_fut)    << " "<<  endl;

 */
FUNCTION Get_Selectivity
  // Calculate the logistic selectivity (Only if being used...)   
  for (k=1;k<=nfsh;k++)
  {
    switch (fsh_sel_opt(k))
    {
      case 1 : // Selectivity coefficients 
      //---Calculate the fishery selectivity from the sel_coffs (Only if being used...)   
      {
        int isel_ch_tmp = 1 ;
        dvar_vector sel_coffs_tmp(1,nselages_fsh(k,isel_ch_tmp));
        for (i=styr;i<=endyr;i++)
        {
          if (i==yrs_sel_ch_fsh(k,isel_ch_tmp)) 
          {
            sel_coffs_tmp.initialize();
            sel_coffs_tmp = log_selcoffs_fsh(k,isel_ch_tmp);
            avgsel_fsh(k,isel_ch_tmp)              = log(mean(mfexp(sel_coffs_tmp)));
            // Increment if there is still space to do so...
            if (isel_ch_tmp<n_sel_ch_fsh(k))
              isel_ch_tmp++;
          }
         // Need to flag for changing selectivity....XXX
          log_sel_fsh(k,i)(1,nselages_fsh(k,isel_ch_tmp))        = sel_coffs_tmp;
          log_sel_fsh(k,i)(nselages_fsh(k,isel_ch_tmp),nages)    = log_sel_fsh(k,i,nselages_fsh(k,isel_ch_tmp));
          log_sel_fsh(k,i)                                  -= log(mean(mfexp(log_sel_fsh(k,i) )));
        }
      }
      break;
      case 2 : // Single logistic
      {
        sel_slope_fsh(k) = mfexp(logsel_slope_fsh(k));
        int isel_ch_tmp = 1 ;
        dvariable sel_slope_tmp = sel_slope_fsh(k,isel_ch_tmp);
        dvariable sel50_tmp     = sel50_fsh(k,isel_ch_tmp);
        for (i=styr;i<=endyr;i++)
        {
          if (i==yrs_sel_ch_fsh(k,isel_ch_tmp)) 
          {
            sel_slope_tmp = sel_slope_fsh(k,isel_ch_tmp);
            sel50_tmp     =     sel50_fsh(k,isel_ch_tmp);
            if (isel_ch_tmp<n_sel_ch_fsh(k))
              isel_ch_tmp++;
          }
          log_sel_fsh(k,i)(1,nselages_fsh(k,isel_ch_tmp))     = -1.*log( 1.0 + mfexp(-1.*sel_slope_tmp * 
                                                ( age_vector(1,nselages_fsh(k,isel_ch_tmp)) - sel50_tmp) ));
          log_sel_fsh(k,i)(nselages_fsh(k,isel_ch_tmp),nages) = log_sel_fsh(k,i,nselages_fsh(k,isel_ch_tmp));
        }
    }
    break;
    case 3 : // Double logistic
    {
      sel_slope_fsh(k)  = mfexp(logsel_slope_fsh(k));
      sel_dslope_fsh(k) = mfexp(logsel_dslope_fsh(k));
      int isel_ch_tmp = 1 ;
      dvariable sel_slope_tmp = sel_slope_fsh(k,isel_ch_tmp);
      dvariable sel50_tmp     = sel50_fsh(k,isel_ch_tmp);
      dvariable sel_dslope_tmp = sel_dslope_fsh(k,isel_ch_tmp);
      dvariable seld50_tmp     = seld50_fsh(k,isel_ch_tmp);
      for (i=styr;i<=endyr;i++)
      {
        if (i==yrs_sel_ch_fsh(k,isel_ch_tmp)) 
        {
          sel_slope_tmp  = sel_slope_fsh(k,isel_ch_tmp);
          sel50_tmp      =     sel50_fsh(k,isel_ch_tmp);
          sel_dslope_tmp = sel_dslope_fsh(k,isel_ch_tmp);
          seld50_tmp     =     seld50_fsh(k,isel_ch_tmp);
          if (isel_ch_tmp<n_sel_ch_fsh(k))
            isel_ch_tmp++;
        }
        log_sel_fsh(k,i)(1,nselages_fsh(k,isel_ch_tmp))     =
                     -log( 1.0 + mfexp(-1.*sel_slope_tmp * 
                     ( age_vector(1,nselages_fsh(k,isel_ch_tmp)) - sel50_tmp) ))+
                     log( 1. - 1/(1.0 + mfexp(-sel_dslope_tmp * 
                     ( age_vector(1,nselages_fsh(k,isel_ch_tmp)) - seld50_tmp))) );

        log_sel_fsh(k,i)(nselages_fsh(k,isel_ch_tmp),nages) = 
                     log_sel_fsh(k,i,nselages_fsh(k,isel_ch_tmp));

        log_sel_fsh(k,i) -= max(log_sel_fsh(k,i));  
      }
    }
    break;
    //---Calculate the fishery selectivity from the sel_spl from nodes...
    case 4 : // Splines
    {
       int isel_ch_tmp = 1 ;
       dvar_matrix tempsel(1,n_sel_ch_fsh(k),1,nselages_fsh(k,1));
       dvector agetmp(1,nselages_fsh(k,1));
       agetmp.fill_seqadd(0.,1.0/(nselages_fsh(k,1)-1));
       // This needs to be dimensioned by the number of changes and the number of coefficients
       vcubic_spline_function_array a(1,n_sel_ch_fsh(k),xnodes_fsh(k),log_sel_spl_fsh(k));
       tempsel = a(agetmp);
       j=1;
       // cout <<tempsel<<endl;exit(1);
       log_sel_fsh(k,styr)(1,nselages_fsh(k,1)) = tempsel(j);;
       avgsel_fsh(k,j)  = mean(log_sel_spl_fsh(k,j)); //log(mean(mfexp(log_sel_fsh(k,styr))));
       log_sel_fsh(k,styr)(nselages_fsh(k,1),nages) = log_sel_fsh(k,styr,nselages_fsh(k,1))  ; 
       log_sel_fsh(k,styr)-= log(mean(mfexp(log_sel_fsh(k,styr) )));
       for (int i=styr+1;i<=endyr;i++) {
         if(sel_change_in_fsh(k,i)>0) {
           avgsel_fsh(k,j)   = mean(log_sel_spl_fsh(k,j)); // log(mean(mfexp(tempsel(j) )));
           j++; 
           log_sel_fsh(k,i)(1,nselages_fsh(k,1)) = tempsel(j);;
           log_sel_fsh(k,i)(nselages_fsh(k,1),nages) = log_sel_fsh(k,i,nselages_fsh(k,1))  ; 
           log_sel_fsh(k,i) -= log(mean(mfexp(log_sel_fsh(k,i) )));
         } 
         else
           log_sel_fsh(k,i) = log_sel_fsh(k,i-1) ; 

       }
     }
     break;
    } // End of switch for fishery selectivity type
  } // End of fishery loop
  // Survey specific---
  for (k=1;k<=nsrv;k++)
  {
    switch (srv_sel_opt(k))
    {
      case 1 : // Selectivity coefficients
      //---Calculate the fishery selectivity from the sel_coffs (Only if being used...)   
      {
        int isel_ch_tmp = 1 ;
  
        dvar_vector sel_coffs_tmp(1,nselages_srv(k,isel_ch_tmp));
        for (i=styr;i<=endyr;i++)
        {
          if (i==yrs_sel_ch_srv(k,isel_ch_tmp)) 
          {
            sel_coffs_tmp.initialize();
            sel_coffs_tmp = log_selcoffs_srv(k,isel_ch_tmp);
            avgsel_srv(k,isel_ch_tmp)              = log(mean(mfexp(sel_coffs_tmp(q_age_min(k),q_age_max(k)))));
            if (isel_ch_tmp<n_sel_ch_srv(k))
              isel_ch_tmp++;
          }
          log_sel_srv(k,i)(1,nselages_srv(k,isel_ch_tmp))        = sel_coffs_tmp;
          log_sel_srv(k,i)(nselages_srv(k,isel_ch_tmp),nages)    = log_sel_srv(k,i,nselages_srv(k,isel_ch_tmp));
          log_sel_srv(k,i)                                  -= log(mean(mfexp(log_sel_srv(k,i)(q_age_min(k),q_age_max(k))))); 
        }
      }
  
        break;
      case 2 : // Asymptotic logistic
        {
          sel_slope_srv(k) = mfexp(logsel_slope_srv(k));
          int isel_ch_tmp = 1 ;
          dvariable sel_slope_tmp = sel_slope_srv(k,isel_ch_tmp);
          dvariable sel50_tmp     = sel50_srv(k,isel_ch_tmp);
          for (i=styr;i<=endyr;i++)
          {
            if (i==yrs_sel_ch_srv(k,isel_ch_tmp)) 
            {
              sel_slope_tmp = sel_slope_srv(k,isel_ch_tmp);
              sel50_tmp     =     sel50_srv(k,isel_ch_tmp);
              if (isel_ch_tmp<n_sel_ch_srv(k))
                isel_ch_tmp++;
            }
            log_sel_srv(k,i) = - log( 1.0 + mfexp(-sel_slope_tmp * ( age_vector - sel50_tmp) ));
            log_sel_srv(k,i)                                  -= log(mean(mfexp(log_sel_srv(k,i)(q_age_min(k),q_age_max(k))))); 
          }
        }
        break;
      case 3 : // Double logistic
        {
          sel_slope_srv(k)  = mfexp(logsel_slope_srv(k));
          sel_dslope_srv(k) = mfexp(logsel_dslope_srv(k));
          int isel_ch_tmp = 1 ;
          dvariable sel_slope_tmp = sel_slope_srv(k,isel_ch_tmp);
          dvariable sel50_tmp     = sel50_srv(k,isel_ch_tmp);
          dvariable sel_dslope_tmp = sel_dslope_srv(k,isel_ch_tmp);
          dvariable seld50_tmp     = seld50_srv(k,isel_ch_tmp);
          for (i=styr;i<=endyr;i++)
          {
            if (i==yrs_sel_ch_srv(k,isel_ch_tmp)) 
            {
              sel_slope_tmp  = sel_slope_srv(k,isel_ch_tmp);
              sel50_tmp      =     sel50_srv(k,isel_ch_tmp);
              sel_dslope_tmp = sel_dslope_srv(k,isel_ch_tmp);
              seld50_tmp     =     seld50_srv(k,isel_ch_tmp);
              if (isel_ch_tmp<n_sel_ch_srv(k))
                isel_ch_tmp++;
            }
            log_sel_srv(k,i)(1,nselages_srv(k,isel_ch_tmp))     =
                         -log( 1.0 + mfexp(-1.*sel_slope_tmp * 
                         ( age_vector(1,nselages_srv(k,isel_ch_tmp)) - sel50_tmp) ))+
                         log( 1. - 1/(1.0 + mfexp(-sel_dslope_tmp * 
                         ( age_vector(1,nselages_srv(k,isel_ch_tmp)) - seld50_tmp))) );

            log_sel_srv(k,i)(nselages_srv(k,isel_ch_tmp),nages) = 
                         log_sel_srv(k,i,nselages_srv(k,isel_ch_tmp));

            log_sel_srv(k,i) -= max(log_sel_srv(k,i));  
            log_sel_srv(k,i)                                  -= log(mean(mfexp(log_sel_srv(k,i)(q_age_min(k),q_age_max(k))))); 
          }
        }
      break;
    }// end of swtiches for survey selectivity
  } // End of survey loop

  // Map selectivities across fisheries and surveys as needed.
  for (k=1;k<=nfsh;k++)
    if (sel_map(2,k)!=k)  // If 2nd row shows a different fishery then use that fishery
      log_sel_fsh(k) = log_sel_fsh(sel_map(2,k));

  for (k=1+nfsh;k<=nfsh_and_srv;k++)
    if (sel_map(1,k)!=2) 
      log_sel_srv(k-nfsh) = log_sel_fsh(sel_map(2,k));
    else if (sel_map(2,k)!=(k-nfsh)) 
      log_sel_srv(k-nfsh) = log_sel_srv(sel_map(2,k));

  sel_fsh = mfexp(log_sel_fsh);
  sel_srv = mfexp(log_sel_srv);

FUNCTION Get_Mortality
  surv    = mfexp(-M);
  natmort = M;
  Z       = M;
  if (!Popes)
  {
    Fmort.initialize();
    for (k=1;k<=nfsh;k++)
    {
      Fmort +=  fmort(k);
      for (i=styr;i<=endyr;i++)
      {
        F(k,i)   =  fmort(k,i) * sel_fsh(k,i) ;
        Z(i)    += F(k,i);
      }
    }
    S  = mfexp(-1.*Z);
  }
  

FUNCTION Get_Numbers_at_Age
  // natage(styr,1) = mfexp(mean_log_rec + rec_dev(styr)); 
  // Recruitment in subsequent years
  for (i=styr+1;i<=endyr;i++)
    natage(i,1)=mfexp(mean_log_rec+rec_dev(i));

  mod_rec(styr)  = natage(styr,1);

  dvariable  t1=mfexp(-natmort*0.5);
  dvariable  t2=mfexp(-natmort);
  for (i=styr;i<=endyr;i++)
  {
    if (Popes)
    {
      Catch_at_Age(i);
      // Pope's approximation //   Next year N     =   This year x NatSurvivl - catch
      natage(i+1)(2,nages) = ++(natage(i)(1,nages-1)*t2 - catage_tot(i)(1,nages-1)*t1);
      Ftot(i)(1,nages-1) = log(natage(i)(1,nages-1)) - --log(natage(i+1)(2,nages)) - natmort;
      natage(i+1,nages)   += natage(i,nages)*t2 - catage_tot(i,nages)*t1;
      // Approximation to "F" continuous form for computing within-year sp biomass
      Ftot(i,nages)      = log(natage(i,nages-1)+natage(i,nages)) -log(natage(i+1,nages)) -natmort;
      // write_input_log <<i<<" "<<Ftot(i)(nages-4,nages)<<endl; // cout <<i<<" "<<natage(i)<<endl; // cout <<i<<" "<<natage(i+1)<<endl;
      dvariable ctmp=sum(catage_tot(i));
      for (k=1;k<=nfsh;k++)
      {
        F(k,i)  = Ftot(i) * sum(catage(k,i))/ctmp;
      }
      Z(i)    = Ftot(i)+natmort;
      S(i)    = mfexp(-Z(i));
    }
    else // Baranov
    {
      // get_Fs( i ); //ojo, add switch here for different catch equation XX
      if (i!=endyr)
      {
        natage(i+1)(2,nages) = ++elem_prod(natage(i)(1,nages-1),S(i)(1,nages-1));
        natage(i+1,nages)+=natage(i,nages)*S(i,nages);
      }
    }
    Catch_at_Age(i);
    Sp_Biom(i)  = elem_prod(natage(i),pow(S(i),spmo_frac)) * wt_mature; 
    if (i<endyr) mod_rec(i+1)  = natage(i+1,1);
  }

FUNCTION Get_Survey_Predictions
  // Survey computations------------------
  dvariable sum_tmp;
  sum_tmp.initialize();
  int yy;
  for (k=1;k<=nsrv;k++)
  {
    for (i=styr;i<=endyr;i++)
    {        
      pred_srv(k,i) = q_srv(k) * pow(elem_prod(natage(i),pow(S(i),srv_month_frac(k))) * 
                                     elem_prod(sel_srv(k,i) , wt_srv(k,i)),q_power_srv(k));
    }
    for (i=1;i<=nyrs_srv_age(k);i++)
    {
      yy = int(yrs_srv_age(k,i)); 
      dvar_vector tmp_n   = elem_prod(pow(S(yy),srv_month_frac(k)),elem_prod(sel_srv(k,yy),natage(yy)));  
      sum_tmp             = sum(tmp_n);
      if (use_age_err)
        eac_srv(k,i)      = age_err * tmp_n/sum_tmp;
      else
        eac_srv(k,i)      = tmp_n/sum_tmp;

    }
  }
FUNCTION Get_Fishery_Predictions
  for (k=1;k<=nfsh;k++)
  {
    for (i=1; i<=nyrs_fsh_age(k); i++)
      if (use_age_err)
        eac_fsh(k,i) = age_err * catage(k,yrs_fsh_age(k,i))/sum(catage(k,yrs_fsh_age(k,i)));
      else
        eac_fsh(k,i) = catage(k,yrs_fsh_age(k,i))/sum(catage(k,yrs_fsh_age(k,i)));
  }

FUNCTION Calc_Dependent_Vars
  // cout<<"In DepVar stage 1"<<endl;
  get_msy();
  if (phase_proj>0) 
    Future_projections();

  dvar_matrix N_nofsh(styr,endyr,1,nages);
  N_nofsh.initialize();
  N_nofsh(styr) = natage(styr);
  for (i=styr;i<=endyr;i++)
  {                 
    recruits(i)  = natage(i,1);
    if (i>styr)
    {
      N_nofsh(i,1) = recruits(i);
      N_nofsh(i)(2,nages) = ++N_nofsh(i-1)(1,nages-1)*exp(-natmort);
      N_nofsh(i,nages)   += N_nofsh(i-1,nages)*exp(-natmort);
    }
    totbiom_NoFish(i) = N_nofsh(i)*wt_pop;
    totbiom(i)        = natage(i)*wt_pop;
    Sp_Biom_NoFishR(i) = Sp_Biom(i) / (N_nofsh(i)*pow(exp(-natmort),spmo_frac) * wt_mature); 
    depletion         = totbiom(endyr)/totbiom(styr);
    depletion_dyn     = totbiom(endyr)/totbiom_NoFish(endyr);
  }
  dvar_vector Ntmp(1,nages);
  Ntmp(2,nages) = ++elem_prod(natage(endyr)(1,nages-1),S(endyr)(1,nages-1));
  Ntmp(nages)  += natage(endyr,nages)*S(endyr,nages);
  Ntmp(1)       = SRecruit(Sp_Biom(endyr+1-rec_age));
  ABCBiom       = Ntmp*wt_pop;
  recruits(endyr+1) = Ntmp(1);
  totbiom(endyr+1)  = ABCBiom;

FUNCTION void Catch_at_Age(const int& i)
  dvariable vbio=0.;
  dvariable pentmp;
  dvar_vector Nmid(1,nages);
  dvar_vector Ctmp(1,nages);
  catage_tot(i).initialize();
  if (Popes)
  {
    Nmid = natage(i)*mfexp(-natmort/2 ); 
  }
  for (k=1;k<=nfsh;k++)
  {
    if (Popes)
    {
      pentmp=0.;
      Ctmp = elem_prod(Nmid,sel_fsh(k,i));
      vbio = Ctmp*wt_fsh(k,i);
      //Kludge to go here...
      // dvariable SK = posfun( (.98*vbio - catch_bio(k,i))/vbio , 0.1 , pentmp );
      dvariable SK = posfun( (vbio - catch_bio(k,i))/vbio , 0.1 , pentmp );
      catch_tmp    = vbio - SK*vbio; 
      hrate        = catch_tmp / vbio;
      fpen(4) += pentmp;
      Ctmp *= hrate;                          
      if (hrate>1) {cout << catch_tmp<<" "<<vbio<<endl;exit(1);}
      catage_tot(i) += Ctmp;                      
      catage(k,i)    = Ctmp;                      
      if (last_phase())
        pred_catch(k,i) = Ctmp*wt_fsh(k,i);
    }
    else
    {
      catage(k,i) = elem_prod(elem_div(F(k,i),Z(i)),elem_prod(1.-S(i),natage(i)));
      pred_catch(k,i) = catage(k,i)*wt_fsh(k,i);
    }
  }
  //+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==
FUNCTION evaluate_the_objective_function
  // if (active(fmort_dev))   
  if (active(fmort))   
  {
    Cat_Like();
    Fmort_Pen();
  }
  Rec_Like();
  if (active(rec_dev))
    Age_Like();
  Srv_Like();
  Sel_Like();
  Compute_priors();
  if (active(log_Rzero)) // OjO
    obj_fun += .5 * square(log_Rzero-mean_log_rec); // A slight penalty to keep Rzero in reality...

  obj_comps.initialize();
  obj_comps(1) = sum(catch_like);
  obj_comps(2) = sum(age_like_fsh);
  obj_comps(3) = sum(sel_like_fsh);
  obj_comps(4) = sum(surv_like);
  obj_comps(5) = sum(age_like_srv);
  obj_comps(6) = sum(sel_like_srv);
  obj_comps(7) = sum(rec_like);
  obj_comps(8) = sum(fpen);
  obj_comps(9) = sum(post_priors_srvq);
  obj_comps(10)= sum(post_priors);
  obj_fun     += sum(obj_comps);

FUNCTION Cat_Like
  // Eases into the catch-biomass likelihoods.  If too far off to start, full constraint to fit can be too aggressive
  catch_like.initialize();
  dvariable catch_pen;
  switch (current_phase())
  {
    case 1:
      catch_pen = .1;
      break;
    case 2:
      catch_pen = .5;
      break;
    case 3:
      catch_pen = .8;
      break;
    case 4:
      catch_pen = 1.0;
      break;
    case 5:
      catch_pen = 1;
      break;
    default:
      catch_pen = 1;
      break;
  }
  if (current_phase()>3)
  {
    for (k=1;k<=nfsh;k++)
      for (i=styr;i<=endyr;i++)
         catch_like(k) += .5*square(log(catch_bio(k,i)+.0001) - log(pred_catch(k,i)+.0001) )/catch_bio_lva(k,i);
  }
  else
  {
    for (k=1;k<=nfsh;k++)
      catch_like(k) += catchbiomass_pen * norm2(log(catch_bio(k)   
                      +.000001) - log(pred_catch(k) +.000001));
  }

  catch_like *= catch_pen;

FUNCTION Rec_Like
  rec_like.initialize();
  if (active(rec_dev))
  {
    sigmar     =  mfexp(log_sigmar);
    sigmarsq   =  square(sigmar);
    dvariable SSQRec;
    SSQRec.initialize();
    if (current_phase()>2)
    {
      if (last_phase())
      {
        pred_rec = SRecruit(Sp_Biom(styr_rec-rec_age,endyr-rec_age).shift(styr_rec)(styr_rec,endyr));
      }
      else 
      {
         pred_rec = 1.+SRecruit(Sp_Biom(styr_rec-rec_age,endyr-rec_age).shift(styr_rec)(styr_rec,endyr));
      }
      dvar_vector chi(styr_rec_est,endyr_rec_est);
      chi = log(mod_rec(styr_rec_est,endyr_rec_est)) - log(pred_rec(styr_rec_est,endyr_rec_est));
      SSQRec   = norm2( chi ) ;
      m_sigmar   =  sqrt( SSQRec  / nrecs_est);
      m_sigmarsq =  m_sigmar * m_sigmar   ;

      if (current_phase()>4||last_phase())
        rec_like(1) = (SSQRec+ sigmarsq/2.)/(2*sigmarsq) + nrecs_est*log_sigmar; 
      else
        rec_like(1) = .1*(SSQRec+ sigmarsq/2.)/(2*sigmarsq) + nrecs_est*log_sigmar; 
    }

    if (last_phase())
    {
      // Variance term for the parts not estimated by sr curve
      rec_like(4) += .5*norm2( rec_dev(styr_rec,styr_rec_est) )/sigmarsq + (styr_rec_est-styr_rec)*log(sigmar) ; 

      if ( endyr > endyr_rec_est)
        rec_like(4) += .5*norm2( rec_dev(endyr_rec_est,endyr  ) )/sigmarsq + (endyr-endyr_rec_est)*log(sigmar) ; 
    }
    else // JNI comment next line
       rec_like(2) += norm2( rec_dev(styr_rec_est,endyr) ) ;

    rec_like(2) += norm2( rec_dev(styr_rec_est,endyr) ) ;

    if (active(rec_dev_future))
    {
      // Future recruitment variability (based on past)
      sigmar_fut   = sigmar ;
      rec_like(3) += norm2(rec_dev_future)/(2*square(sigmar_fut))+ size_count(rec_dev_future)*log(sigmar_fut);
    }
  }

FUNCTION Compute_priors
  post_priors.initialize();
  post_priors_srvq.initialize();
  for (k=1;k<=nsrv;k++)
  {
    if (active(log_q_srv(k)))
      post_priors_srvq(k) += square(q_srv(k)-qprior(k))/(2*cvqprior(k)*cvqprior(k)); 
    if (active(log_q_power_srv(k)))
      post_priors_srvq(k) += square(q_power_srv(k)-q_power_prior(k))/(2*cvq_power_prior(k)*cvq_power_prior(k)); 
  }

  if (active(M))
    post_priors(1) += square(M-natmortprior)/(2*cvnatmortprior*cvnatmortprior); 

  if (active(steepness))
    post_priors(2) += square(steepness-steepnessprior)/(2*cvsteepnessprior*cvsteepnessprior); 

  if (active(log_sigmar))
    post_priors(3) += square(sigmar-sigmarprior)/(2*cvsigmarprior*cvsigmarprior); 

FUNCTION Fmort_Pen
  // Phases less than 3, penalize High F's---------------------------------
  if (current_phase()<3)
    fpen(1) += 1.* norm2(F - .2);
  else 
    fpen(1) += 0.0001*norm2(F - .2); 

  // for (k=1;k<=nfsh;k++)  fpen(2) += 20.*square(mean(fmort_dev(k)) ); // this is just a normalizing constraint (fmort_devs sum to zero) }
    
FUNCTION Sel_Like 
  sel_like_fsh.initialize();
  sel_like_srv.initialize();
  for (k=1;k<=nfsh;k++)
  {
    if (active(log_selcoffs_fsh(k)))
    {
      for (i=1;i<=n_sel_ch_fsh(k);i++)
      {
        int iyr = yrs_sel_ch_fsh(k,i) ;
        sel_like_fsh(k,1) += curv_pen_fsh(k)*norm2(first_difference(
                                                   first_difference(log_sel_fsh(k,iyr ))));
        // This part is the penalty on the change itself--------------
        if (i>1)
        {
          dvariable var_tmp = square(sel_change_in_fsh(k,iyr ));
          sel_like_fsh(k,2)    += .5*norm2( log_sel_fsh(k,iyr-1) - log_sel_fsh(k,iyr) ) / var_tmp ;
        }
        int nagestmp = nselages_fsh(k,1);
        for (j=seldecage;j<=nagestmp;j++)
        {
          dvariable difftmp = log_sel_fsh(k,iyr,j-1)-log_sel_fsh(k,iyr,j) ;
          if (difftmp > 0.)
            sel_like_fsh(k,3)    += .5*square( difftmp ) / seldec_pen_fsh(k);
        }
        obj_fun            += 20 * square(avgsel_fsh(k,i)); // To normalize selectivities
      }
    }
   if (active(log_sel_spl_fsh(k))	)
    {
      for (i=1;i<=n_sel_ch_fsh(k);i++)
      {
        int iyr = yrs_sel_ch_fsh(k,i) ;
        // This part is the penalty on the change itself--------------
        if (i>1)
        {
          dvariable var_tmp = square(sel_change_in_fsh(k,iyr ));
          sel_like_fsh(k,2)    += .5*norm2( log_sel_fsh(k,iyr-1) - log_sel_fsh(k,iyr) ) / var_tmp ;
        }
        int nagestmp = nselages_fsh(k,1);
        /*for (j=seldecage;j<=nagestmp;j++)
        {
          dvariable difftmp = log_sel_fsh(k,iyr,j-1)-log_sel_fsh(k,iyr,j) ;
          if (difftmp > 0.)
            sel_like_fsh(k,3)    += .5*square( difftmp ) / seldec_pen_fsh(k);
        }*/
        obj_fun            += 20. * square(avgsel_fsh(k,i)); // To normalize selectivities
        // obj_fun            += 20. * square(log_sel_fsh(k,i,nages)); // To normalize selectivities
      }
    }
  }

  for (k=1;k<=nsrv;k++)
  {
    if (active(log_selcoffs_srv(k)))
    {
      for (i=1;i<=n_sel_ch_srv(k);i++)
      {
        int iyr = yrs_sel_ch_srv(k,i) ;
        sel_like_srv(k,1) += curv_pen_srv(k)*norm2(first_difference(
                                                 first_difference(log_sel_srv(k,iyr))));
        // This part is the penalty on the change itself--------------
        if (i>1)
        {
          dvariable var_tmp = square(sel_change_in_srv(k,iyr ));
          sel_like_srv(k,2)    += .5*norm2( log_sel_srv(k,iyr-1) - log_sel_srv(k,iyr) ) 
                                   / var_tmp ;
        }
        int nagestmp = nselages_srv(k,1);
        for (j=seldecage;j<=nagestmp;j++)
        {
          dvariable difftmp = log_sel_srv(k,iyr,j-1)-log_sel_srv(k,iyr,j) ;
          if (difftmp > 0.)
            sel_like_srv(k,3)    += .5*square( difftmp ) / seldec_pen_srv(k);
        }
        obj_fun            += 20. * square(avgsel_srv(k,i));  // To normalize selectivities
      }
    }
  }

FUNCTION Srv_Like
  // Fit to indices (log-Normal) -------------------------------------------
  surv_like.initialize();
  int yy;
  for (k=1;k<=nsrv;k++)
    for (i=1;i<=nyrs_srv(k);i++)
    {
      yy = int(yrs_srv(k,i));
      surv_like(k) += square(log(obs_srv(k,i)) - log(pred_srv(k,yy)) ) / 
                                   (2.*obs_lse_srv(k,i)*obs_lse_srv(k,i));
    }
  /* normal distribution option to add someday...
    for (i=1;i<=nyrs_srv(k);i++)
      surv_like(k) += square(obs_srv(k,i) - pred_srv(k,yrs_srv(k,i)) ) / 
                                   (2.*obs_se_srv(k,i)*obs_se_srv(k,i));
  */

FUNCTION Age_Like
  age_like_fsh.initialize();
  for (k=1;k<=nfsh;k++)
    for (int i=1;i<=nyrs_fsh_age(k);i++)
      age_like_fsh(k) -= n_sample_fsh_age(k,i)*(oac_fsh(k,i) + 0.001) * log(eac_fsh(k,i) + 0.001 ) ;
  age_like_fsh -= offset_fsh;

  age_like_srv.initialize();
  for (k=1;k<=nsrv;k++)
    for (int i=1;i<=nyrs_srv_age(k);i++)
      age_like_srv(k) -= n_sample_srv_age(k,i)*(oac_srv(k,i) + 0.001) * log(eac_srv(k,i) + 0.001 ) ;
  age_like_srv -= offset_srv;

FUNCTION Oper_Model
 // Initialize things used here only
  Calc_Dependent_Vars();
  Write_Datafile();

  int nsims;
  ifstream sim_in("nsims.dat");
  sim_in >> nsims; sim_in.close();

  dvector ran_srv_vect(1,nsrv);
  ofstream SaveOM("Om_Out.dat",ios::app);
  double C_tmp;
  dvariable Fnow;
  dvariable meanrec;
  meanrec=mean(recruits);
  dvector new_srv(1,nsrv);
  new_srv.initialize();
  system("cls"); cout<<"Number of replicates: "<<endl;
  // Initialize recruitment in first year
  nage_future(styr_fut,1) = meanrec;
  nage_future(styr_fut)(2,nages)              = ++elem_prod(natage(endyr)(1,nages-1),S(endyr)(1,nages-1));
  nage_future(styr_fut,nages)                += natage(endyr,nages)*S(endyr,nages);
  for (int isim=1;isim<=nsims;isim++)
  {
    cout<<isim<<" ";
    // Copy file to get mean for Mgt Strategies
    system("init_stuff.bat");
    for (i=styr_fut;i<=endyr_fut;i++)
    {
      // Some unit normals...
      ran_srv_vect.fill_randn(rng);
      // Create new survey observations
      for (k = 1 ; k<= nsrv ; k++)
        new_srv(k) = mfexp(ran_srv_vect(k)*.2)*value(nage_future(i)*q_srv(k)*sel_srv(k,endyr)); // use value function since converts to a double
      // Append new survey observation to datafile
      ofstream srv_out("NewSrv.dat",ios::app);
      srv_out <<i<<" "<< new_srv<<endl; srv_out.close();
      system("ComputeTAC.bat"); // commandline function to get TAC (catchnext.dat)

      // Now read in TAC (actual catch)
      ifstream CatchNext("CatchNext.dat");
      CatchNext >> C_tmp; CatchNext.close();
      Fnow = SolveF2(endyr,nage_future(i), C_tmp);

      F_future(1,i) = sel_fsh(1,endyr) * Fnow;
      Z_future(i)   = F_future(1,i) + natmort;
      S_future(i)   = mfexp(-Z_future(i));
      // nage_future(i,1)  = SRecruit( Sp_Biom_future(i-rec_age) ) * mfexp(rec_dev_future(i)) ;     

      nage_future(i,1)  = meanrec;
      Sp_Biom_future(i) = wt_mature * elem_prod(nage_future(i),pow(S_future(i),spmo_frac)) ;
      // Now graduate for the next year....
      if (i<endyr_fut)
      {
        nage_future(i+1)(2,nages) = ++elem_prod(nage_future(i)(1,nages-1),S_future(i)(1,nages-1));
        nage_future(i+1,nages)   += nage_future(i,nages)*S_future(i,nages);
      }
      catage_future(i) = 0.; 
      for (k = 1 ; k<= nfsh ; k++)
        catage_future(i) += elem_prod(nage_future(i) , elem_prod(F_future(k,i) , elem_div( ( 1.- S_future(i) ) , Z_future(i))));
  
      SaveOM << model_name       <<
        " "  << isim             <<
        " "  << i                <<
        " "  << Fnow             <<
        " "  << nage_future(i)(5,nages)*wt_fsh(1,endyr)(5,nages)<<
        " "  << Sp_Biom_future(i)<<
        " "  << catage_future(i)*wt_fsh(1,endyr)<<
        " "  << nage_future(i)<<
        " "  << F_future(1,i)<<
      endl;
    }
  }
  SaveOM.close();
  if (!mceval_phase())
    exit(1);

FUNCTION void get_future_Fs(const int& i,const int& iscenario)
    // get F's
    switch (iscenario)
    {
      case 1:
        f_tmp = F35;
        cout <<"F35 "<<i<<" "<<f_tmp<<endl;
        break;
      case 2:
        // f_tmp = SolveF2(endyr,nage_future(i), 1.0  * sum(catch_bioT(endyr)) );
        for (int k=1;k<=nfsh;k++) f_tmp(k) = mean(F(k,endyr));
        break;
      case 3:
        // f_tmp = SolveF2(endyr,nage_future(i), 0.5 * sum(catch_bioT(endyr)) );
        for (int k=1;k<=nfsh;k++) f_tmp(k) = .5*mean(F(k,endyr));
        break;
      case 4:
        // f_tmp = SolveF2(endyr,nage_future(i), 0.25 * sum(catch_bioT(endyr)) );
        for (int k=1;k<=nfsh;k++) f_tmp(k) = .25*mean(F(k,endyr));
        break;
      case 5:
        f_tmp = 0.0;
        break;
    }
    Z_future(i) = natmort;
    for (k=1;k<=nfsh;k++)
    {
      F_future(k,i) = sel_fsh(k,endyr) * f_tmp(k);
      Z_future(i)  += F_future(k,i);
    }
    S_future(i) = mfexp(-Z_future(i));

FUNCTION Future_projections
  // Need to check on treatment of Fratio--whether it should be included or not
  future_biomass.initialize();
  catch_future.initialize();
  for (int iscen=1;iscen<=5;iscen++)
  {
   // Future Sp_Biom set equal to estimated Sp_Biom w/ right lag
    // Sp_Biom_future(styr_fut-rec_age,styr_fut-1) = Sp_Biom(endyr-rec_age+1,endyr);
    for (i=styr_fut-rec_age;i<=styr_fut-1;i++)
      Sp_Biom_future(i) = wt_mature * elem_prod(natage(i),pow(S(i),spmo_frac)) ;

    // cout<<Sp_Biom(endyr-10,endyr)<<endl<<Sp_Biom_future<<endl;exit(1);
    nage_future(styr_fut)(2,nages) = ++elem_prod(natage(endyr)(1,nages-1),S(endyr)(1,nages-1));
    nage_future(styr_fut,nages)   += natage(endyr,nages)*S(endyr,nages);
    Sp_Biom_future(styr_fut)       = wt_mature * elem_prod(nage_future(i),pow(S_future(i),spmo_frac)) ;
    // Future Recruitment (and Sp_Biom)
    for (i=styr_fut;i<endyr_fut;i++)
    {
      nage_future(i,1)  = SRecruit( Sp_Biom_future(i-rec_age) ) * mfexp(rec_dev_future(i)) ;     
      get_future_Fs(i,iscen);
      // Now graduate for the next year....
      nage_future(i+1)(2,nages) = ++elem_prod(nage_future(i)(1,nages-1),S_future(i)(1,nages-1));
      nage_future(i+1,nages)   += nage_future(i,nages)*S_future(i,nages);
      Sp_Biom_future(i) = wt_mature * elem_prod(nage_future(i),pow(S_future(i),spmo_frac)) ;
    }
    nage_future(endyr_fut,1)  = SRecruit( Sp_Biom_future(endyr_fut-rec_age) ) * mfexp(rec_dev_future(endyr_fut)) ;     
    get_future_Fs(endyr_fut,iscen);
    Sp_Biom_future(endyr_fut)  = wt_mature * elem_prod(nage_future(endyr_fut),pow(S_future(endyr_fut),spmo_frac)) ;
    // Now get catch at future ages
    for (i=styr_fut; i<=endyr_fut; i++)
    {
      catage_future(i) = 0.; 
      for (k = 1 ; k<= nfsh ; k++)
      {
        catage_future(i) += elem_prod(nage_future(i) , elem_prod(F_future(k,i) , 
                              elem_div( ( 1.- S_future(i) ) , Z_future(i))));
        if (iscen!=5) 
          catch_future(iscen,i)   += catage_future(i)*wt_fsh(k,endyr);
      }
      future_biomass(iscen,i) = Sp_Biom_future(i);
      cout<<iscen<<" "<<i<<" "<< future_biomass(iscen,i) <<" "<<nage_future(i)(1,4)<<endl;
    }
  }   //End of loop over F's
  Sp_Biom(endyr+1) = Sp_Biom_future(endyr+1);

FUNCTION get_msy
 /*Function calculates used in calculating MSY and MSYL for a designated component of the
  population, given values for stock recruitment and selectivity...  
  Fmsy is the trial value of MSY example of the use of "funnel" to reduce the amount of storage for derivative calculations */

  dvariable sumF=0.;
  for (k=1;k<=nfsh;k++)
    sumF += sum(F(k,endyr));
  for (k=1;k<=nfsh;k++)
    Fratio(k) = sum(F(k,endyr)) / sumF;

  dvariable Stmp;
  dvariable Rtmp;
  double df=1.e-05;
  dvariable F1;
  F1.initialize();
  F1 = (1.1*natmortprior);
  dvariable F2;
  dvariable F3;
  dvariable yld1;
  dvariable yld2;
  dvariable yld3;
  dvariable dyld;
  dvariable dyldp;
  int breakout=0;
  // Newton Raphson stuff to go here
  for (int ii=1;ii<=8;ii++)
  {
    if (mceval_phase()&&(F1>5||F1<0.01)) 
    {
      ii=8;
      if (F1>5) F1=5.0; 
      else      F1=0.001; 
      breakout    = 1;
    }
    F2     = F1 + df*.5;
    F3     = F2 - df;
    // yld1   = yield(Fratio,F1, Stmp,Rtmp); // yld2   = yield(Fratio,F2,Stmp,Rtmp); // yld3   = yield(Fratio,F3,Stmp,Rtmp);
    yld1   = yield(Fratio,F1);
    yld2   = yield(Fratio,F2);
    yld3   = yield(Fratio,F3);
    dyld   = (yld2 - yld3)/df;                          // First derivative (to find the root of this)
    dyldp  = (yld2 + yld3 - 2.*yld1)/(.25*df*df);   // Second derivative (for Newton Raphson)
    if (breakout==0)
    {
      F1    -= dyld/dyldp;
    }
    else
    {
      if (F1>5) 
        cout<<"Fmsy v. high "<< endl;// yld1<<" "<< yld2<<" "<< yld3<<" "<< F1<<" "<< F2<<" "<< F3<<" "<< endl;
      else      
        cout<<"Fmsy v. low "<< endl;// yld1<<" "<< yld2<<" "<< yld3<<" "<< F1<<" "<< F2<<" "<< F3<<" "<< endl;
    }
  }
  {
    dvar_vector ttt(1,5);
    ttt      = yld(Fratio,F1);
    Fmsy     = F1;
    MSY      = ttt(2);
    Bmsy     = ttt(1);
    MSYL     = ttt(1)/Bzero;
    lnFmsy   = log(MSY/ttt(5)); // Exploitation fraction relative to total biomass
    Bcur_Bmsy= Sp_Biom(endyr)/Bmsy;

    dvariable FFtmp;
    FFtmp.initialize();
    for (k=1;k<=nfsh;k++)
      FFtmp += mean(F(k,endyr));
    Fcur_Fmsy= FFtmp/Fmsy;
    Rmsy     = Rtmp;
  }

 //+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+ 
FUNCTION dvar_vector yld(const dvar_vector& Fratio, const dvariable& Ftmp)
  RETURN_ARRAYS_INCREMENT();
  /*dvariable utmp=1.-mfexp(-(Ftmp)); dvariable Ntmp; dvariable Btmp; dvariable yield; dvariable survtmp=exp(-1.*natmort); dvar_vector seltmp=sel_fsh(endyr); Ntmp = 1.; Btmp = Ntmp*wt(1)*seltmp(1); Stmp = .5*Ntmp*wt(1)*maturity(1); yield= 0.; for ( j=1 ; j < nages ; j++ ) { Ntmp  *= (1.-utmp*seltmp(j))*survtmp; Btmp  += Ntmp*wt(j+1)*seltmp(j+1); Stmp  += .5 * Ntmp *wt(j+1)*maturity(j+1); } //Max Age - 1 yr yield   += utmp * Btmp; Ntmp    /= (1-survtmp*(1.-utmp*seltmp(nages))); Btmp    += Ntmp*wt(nages)*seltmp(nages); Stmp    += 0.5 *wt(nages)* Ntmp *maturity(nages); yield   += utmp * Btmp; //cout<<yield<<" "<<Stmp<<" "<<Btmp<<" ";*/
  dvar_vector msy_stuff(1,5);
  dvariable phi;
  dvar_vector Ntmp(1,nages);
  dvar_vector Ctmp(1,nages);
  msy_stuff.initialize();

  dvar_matrix seltmp(1,nfsh,1,nages);
  for (k=1;k<=nfsh;k++)
   seltmp(k) = sel_fsh(k,endyr); // NOTE uses last-year of fishery selectivity for projections.

  dvar_matrix Fatmp(1,nfsh,1,nages);
  dvar_vector Ztmp(1,nages);

  Ztmp = natmort;
  for (k=1;k<=nfsh;k++)
  { 
    Fatmp(k) = Fratio(k) * Ftmp * seltmp(k);
    Ztmp    += Fatmp(k);
  } 
  dvar_vector survtmp = mfexp(-Ztmp);

  Ntmp(1) = 1.;
  for ( j=1 ; j < nages; j++ )
    Ntmp(j+1)  =   Ntmp(j) * survtmp(j); // Begin numbers in the next year/age class
  Ntmp(nages)  /= (1.- survtmp(nages)); 

  for (k=1;k<=nfsh;k++)
  {
    Ctmp.initialize();
    for ( j=1 ; j <= nages; j++ )
      Ctmp(j)      = Ntmp(j) * Fatmp(k,j) * (1. - survtmp(j)) / Ztmp(j);

    msy_stuff(2)  += wt_fsh(k,endyr) * Ctmp;
  }
  phi    = elem_prod( Ntmp , pow(survtmp,spmo_frac ) ) * wt_mature;
  // Req    = Requil(phi) * exp(sigmarsq/2);
  msy_stuff(5)  = Ntmp * wt_pop;      
  msy_stuff(4)  = phi/phizero ;       // SPR
  msy_stuff(3)  = Requil(phi) ;       // Eq Recruitment
  msy_stuff(5) *= msy_stuff(3);       // BmsyTot
  msy_stuff(2) *= msy_stuff(3);       // MSY
  msy_stuff(1)  = phi*(msy_stuff(3)); // Bmsy
  RETURN_ARRAYS_DECREMENT();
  return msy_stuff;

FUNCTION dvariable yield(const dvar_vector& Fratio, const dvariable& Ftmp)
  RETURN_ARRAYS_INCREMENT();
  /*dvariable utmp=1.-mfexp(-(Ftmp)); dvariable Ntmp; dvariable Btmp; dvariable yield; dvariable survtmp=exp(-1.*natmort); dvar_vector seltmp=sel_fsh(endyr); Ntmp = 1.; Btmp = Ntmp*wt(1)*seltmp(1); Stmp = .5*Ntmp*wt(1)*maturity(1); yield= 0.; for ( j=1 ; j < nages ; j++ ) { Ntmp  *= (1.-utmp*seltmp(j))*survtmp; Btmp  += Ntmp*wt(j+1)*seltmp(j+1); Stmp  += .5 * Ntmp *wt(j+1)*maturity(j+1); } //Max Age - 1 yr yield   += utmp * Btmp; Ntmp    /= (1-survtmp*(1.-utmp*seltmp(nages))); Btmp    += Ntmp*wt(nages)*seltmp(nages); Stmp    += 0.5 *wt(nages)* Ntmp *maturity(nages); yield   += utmp * Btmp; //cout<<yield<<" "<<Stmp<<" "<<Btmp<<" ";*/
  dvariable phi;
  dvariable Req;
  dvar_vector Ntmp(1,nages);
  dvar_vector Ctmp(1,nages);
  dvariable   yield;
  yield.initialize();

  dvar_matrix seltmp(1,nfsh,1,nages);
  for (k=1;k<=nfsh;k++)
   seltmp(k) = sel_fsh(k,endyr); // NOTE uses last-year of fishery selectivity for projections.

  dvar_matrix Fatmp(1,nfsh,1,nages);
  dvar_vector Ztmp(1,nages);

  Ztmp = natmort;
  for (k=1;k<=nfsh;k++)
  { 
    Fatmp(k) = Fratio(k) * Ftmp * seltmp(k);
    Ztmp    += Fatmp(k);
  } 
  dvar_vector survtmp = mfexp(-Ztmp);
  // cout<<Ftmp<<" ";

  Ntmp(1) = 1.;
  for ( j=1 ; j < nages; j++ )
    Ntmp(j+1)  =   Ntmp(j) * survtmp(j); // Begin numbers in the next year/age class
  Ntmp(nages)  /= (1.- survtmp(nages)); 

  for (k=1;k<=nfsh;k++)
  {
    Ctmp.initialize();
    for ( j=1 ; j <= nages; j++ )
      Ctmp(j)      = Ntmp(j) * Fatmp(k,j) * (1. - survtmp(j)) / Ztmp(j);

    yield  += wt_fsh(k,endyr) * Ctmp;
  }
  phi    = elem_prod( Ntmp , pow(survtmp,spmo_frac ) )* wt_mature;
  // Req    = Requil(phi) * mfexp(sigmarsq/2);
  Req    = Requil(phi) ;
  yield *= Req;

  RETURN_ARRAYS_DECREMENT();
  return yield;

FUNCTION dvariable yield(const dvar_vector& Fratio, dvariable& Ftmp, dvariable& Stmp,dvariable& Req)
  RETURN_ARRAYS_INCREMENT();
  dvariable phi;
  dvar_vector Ntmp(1,nages);
  dvar_vector Ctmp(1,nages);
  dvariable   yield   = 0.;

  dvar_matrix seltmp(1,nfsh,1,nages);
  for (k=1;k<=nfsh;k++)
   seltmp(k) = sel_fsh(k,endyr); // NOTE uses last-year of fishery selectivity for projections.

  dvar_matrix Fatmp(1,nfsh,1,nages);
  dvar_vector Ztmp(1,nages);

  Ztmp = natmort;
  for (k=1;k<=nfsh;k++)
  { 
    Fatmp(k) = Fratio(k) * Ftmp * seltmp(k);
    Ztmp    += Fatmp(k);
  } 
  dvar_vector survtmp = mfexp(-Ztmp);

  Ntmp(1) = 1.;
  for ( j=1 ; j < nages; j++ )
    Ntmp(j+1)  =   Ntmp(j) * survtmp(j); // Begin numbers in the next year/age class
  Ntmp(nages)  /= (1.- survtmp(nages)); 
  for (k=1;k<=nfsh;k++)
  {
    Ctmp.initialize();
    for ( j=1 ; j <= nages; j++ )
      Ctmp(j)      = Ntmp(j) * Fatmp(k,j) * (1. - survtmp(j)) / Ztmp(j);

    yield  += wt_fsh(k,endyr) * Ctmp;
  }
  phi    = elem_prod( Ntmp , pow(survtmp,spmo_frac ) )* wt_mature;
  // Req    = Requil(phi) * exp(sigmarsq/2);
  Req    = Requil(phi) ;
  yield *= Req;
  Stmp   = phi*Req;

  RETURN_ARRAYS_DECREMENT();
  return yield;

FUNCTION Profile_F
  cout << "Doing a profile over F...."<<endl;
  ofstream prof_F("Fprof.yld");
 /* NOTE THis will need to be conditional on SrType too Function calculates 
  used in calculating MSY and MSYL for a designated component of the
  population, given values for stock recruitment and selectivity...  
  Fmsy is the trial value of MSY example of the use of "funnel" to 
  reduce the amount of storage for derivative calculations */
 dvariable sumF=0.;
  for (k=1;k<=nfsh;k++)
    sumF += sum(F(k,endyr));
  for (k=1;k<=nfsh;k++)
    Fratio(k) = sum(F(k,endyr)) / sumF;
  dvariable Stmp;
  dvariable Rtmp;
  double df=1.e-7;
  dvariable F1=.05;
  dvariable F2;
  dvariable F3;
  dvariable yld1;
  dvariable yld2;
  dvariable yld3;
  dvariable dyld;
  dvariable dyldp;
  prof_F <<"Profile of stock, yield, and recruitment over F"<<endl;
  prof_F << model_name<<" "<<datafile_name<<endl;
  prof_F <<endl<<endl<<"F  Stock  Yld  Recruit SPR"<<endl;
  prof_F <<0.0<<" "<< Bzero <<" "<<0.0<<" "<<Rzero<< " 1.00"<<endl; 
  dvar_vector ttt(1,5);
  for (int ii=1;ii<=500;ii++)
  {
    F1    = double(ii)/500;
    yld1  = yield(Fratio,F1,Stmp,Rtmp);
    ttt   = yld(Fratio,F1);
    prof_F <<F1<<" "<< ttt << endl; 
  } 

FUNCTION dvar_vector SRecruit(const dvar_vector& Stmp)
  RETURN_ARRAYS_INCREMENT();
  dvar_vector RecTmp(Stmp.indexmin(),Stmp.indexmax());
  switch (SrType)
  {
    case 1:
      RecTmp = elem_prod((Stmp / phizero) , mfexp( alpha * ( 1. - Stmp / Bzero ))) ; //Ricker form from Dorn
      break;
    case 2:
      RecTmp = elem_prod(Stmp , 1. / ( alpha + beta * Stmp));        //Beverton-Holt form
      break;
    case 3:
      RecTmp = mfexp(mean_log_rec);                    //Avg recruitment
      break;
    case 4:
      RecTmp = elem_prod(Stmp , mfexp( alpha  - Stmp * beta)) ; //Old Ricker form
      break;
  }
  RETURN_ARRAYS_DECREMENT();
  return RecTmp;

FUNCTION dvariable SRecruit(const double& Stmp)
  RETURN_ARRAYS_INCREMENT();
  dvariable RecTmp;
  switch (SrType)
  {
    case 1:
      RecTmp = (Stmp / phizero) * mfexp( alpha * ( 1. - Stmp / Bzero )) ; //Ricker form from Dorn
      break;
    case 2:
      RecTmp = Stmp / ( alpha + beta * Stmp);        //Beverton-Holt form
      break;
    case 3:
      RecTmp = mfexp(mean_log_rec);                    //Avg recruitment
      break;
    case 4:
      RecTmp = Stmp * mfexp( alpha  - Stmp * beta) ; //old Ricker form
      break;
  }
  RETURN_ARRAYS_DECREMENT();
  return RecTmp;

FUNCTION dvariable SRecruit(_CONST dvariable& Stmp)
  RETURN_ARRAYS_INCREMENT();
  dvariable RecTmp;
  switch (SrType)
  {
    case 1:
      RecTmp = (Stmp / phizero) * mfexp( alpha * ( 1. - Stmp / Bzero )) ; //Ricker form from Dorn
      break;
    case 2:
      RecTmp = Stmp / ( alpha + beta * Stmp);        //Beverton-Holt form
      break;
    case 3:
      RecTmp = mfexp(mean_log_rec );                    //Avg recruitment
      break;
    case 4:
      RecTmp = Stmp * mfexp( alpha  - Stmp * beta) ; //old Ricker form
      break;
  }
  RETURN_ARRAYS_DECREMENT();
  return RecTmp;

 //=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
FUNCTION Get_Bzero
  Bzero.initialize();
  Rzero    =  mfexp(log_Rzero); 

  dvariable survtmp = mfexp(-natmort);
  surv = survtmp;

  dvar_matrix natagetmp(styr_rec,styr,1,nages);
  natagetmp.initialize();

  natagetmp(styr_rec,1) = Rzero;
  for (j=2; j<=nages; j++)
    natagetmp(styr_rec,j) = natagetmp(styr_rec,j-1) * survtmp;
  natagetmp(styr_rec,nages) /= (1.-survtmp); 

  Bzero = wt_mature * pow(survtmp,spmo_frac)*natagetmp(styr_rec) ;
  phizero = Bzero/Rzero;

  switch (SrType)
  {
    case 1:
      alpha = log(-4.*steepness/(steepness-1.));
      break;
    case 2:
    {
      alpha  =  Bzero * (1. - (steepness - 0.2) / (0.8*steepness) ) / Rzero;
      beta   = (5. * steepness - 1.) / (4. * steepness * Rzero);
    }
    break;
    case 4:
    {
      beta  = log(5.*steepness)/(0.8*Bzero) ;
      alpha = log(Rzero/Bzero)+beta*Bzero;
    }
      break;
  }
  Sp_Biom.initialize();
  Sp_Biom(styr_sp,styr_rec-1) = Bzero;
  for (i=styr_rec;i<styr;i++)
  {
    Sp_Biom(i) = natagetmp(i)*pow(surv,spmo_frac) * wt_mature; 
    // natagetmp(i,1)          = mfexp(rec_dev(i) + log_Rzero); // OjO numbers a function of mean not SR curve...
    natagetmp(i,1)          = mfexp(rec_dev(i) + mean_log_rec);
    natagetmp(i+1)(2,nages) = ++(natagetmp(i)(1,nages-1)*mfexp(-natmort ));
    natagetmp(i+1,nages)   += natagetmp(i,nages)*mfexp(-natmort);
  }
  natagetmp(styr,1)   = mfexp(rec_dev(styr) + mean_log_rec);
  mod_rec(styr_rec,styr) = column(natagetmp,1);
  natage(styr)  = natagetmp(styr); // OjO
  Sp_Biom(styr) = natagetmp(styr)*pow(surv,spmo_frac) * wt_mature; 
  // cout <<natagetmp<<endl;exit(1);

FUNCTION dvariable Requil(dvariable& phi)
  RETURN_ARRAYS_INCREMENT();
  dvariable RecTmp;
  switch (SrType)
  {
    case 1:
      RecTmp =  Bzero * (alpha + log(phi) - log(phizero) ) / (alpha*phi);
      break;
    case 2:
      RecTmp =  (phi-alpha)/(beta*phi);
      break;
    case 3:
      RecTmp =  mfexp(mean_log_rec);
      break;
    case 4:
      RecTmp =  (log(phi)+alpha) / (beta*phi); //RecTmp =  (log(phi)/alpha + 1.)*beta/phi;
      break;
  }
  // Req    = Requil(phi) * exp(sigmarsq/2);
  // return RecTmp* exp(sigmarsq/2);
  RETURN_ARRAYS_DECREMENT();
  return RecTmp;

FUNCTION write_mceval_hdr
    for (k=1;k<=nsrv;k++)
      mceval<< " q_srv_"<< k<< " ";
    mceval<<"M steepness depletion MSY MSYL Fmsy Fcur_Fmsy Bcur_Bmsy Bmsy totbiom_"<<endyr<<" "<< 
    " F35          "<< 
    " F40          "<< 
    " F50          "<< 
    " fut_SPB_Fmsy_"<< endyr_fut<<" "<< 
    " fut_SPB_F50%_"<< endyr_fut<<" "<< 
    " fut_SPB_F40%_"<< endyr_fut<<" "<< 
    " fut_SPB_F35%_"<< endyr_fut<<" "<< 
    " fut_SPB_F0_"  << endyr_fut<<" "<< 
    " fut_catch_Fmsy_"<<styr_fut<<" "<<  
    " fut_catch_F50%_"<<styr_fut<<" "<<  
    " fut_catch_F40%_"<<styr_fut<<" "<<  
    " fut_catch_F35%_"<<styr_fut<<" "<<  endl;

//+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+==+ 
REPORT_SECTION

  if (last_phase())
  {
    int nvar1=initial_params::nvarcalc(); // get the number of active parameters
    int ndvar=stddev_params::num_stddev_calc();
    int offset=1;
    dvector param_values(1,nvar1+ndvar);
    initial_params::copy_all_values(param_values,offset);
    for (int i=0;i<initial_params::num_initial_params;i++)
    {
      cout << "# " << initial_params::varsptr[i]->label() << "\n" << endl; 
      if (withinbound(0,(initial_params::varsptr[i])->phase_start, initial_params::current_phase))
      {
        int sc = (initial_params::varsptr[i])->size_count();
        if (sc>0)
        {
          write_input_log << "# " << initial_params::varsptr[i]->label() << endl<<param_values(i)<<"\n" << endl; 
        } 
      }
    }
        
    for (k=1;k<=nsrv;k++)
    {
      cout<<srvname(k)<<endl;
      cout<<get_AC(k)<<endl<<endl;
    }
    if (!Popes)
      for (k=1;k<=nfsh;k++)
        Ftot += F(k);
    log_param(M);
    log_param(mean_log_rec);
    log_param(steepness);
    log_param(log_Rzero);
    log_param(rec_dev);
    log_param(log_sigmar);
    log_param(fmort);
    // log_param(log_selcoffs_fsh);
    // log_param(log_sel_spl_fsh);
    // log_param(logsel_slope_fsh);
    // log_param(sel50_fsh);
    // log_param(logsel_dslope_fsh);
    // log_param(seld50_fsh);
    log_param(rec_dev_future);
    // log_param(log_q_srv);
    // log_param(log_q_power_srv);
    // log_param(log_selcoffs_srv);
    // log_param(logsel_slope_srv);
    // log_param(logsel_dslope_srv);
    // log_param(sel50_srv);
    // log_param(seld50_srv);
  }
    
  if (oper_mod)
    Oper_Model();

  cout <<"==============================================================="<<endl;
  if(last_phase())
    cout<<"||  ++++++ Completed phase "<<current_phase()<<" In last phase now +++++"<< endl<<"||"<<endl<<"||  "<<cntrlfile_name <<endl;
  else
    cout<<"||  ++++++ Completed phase "<<current_phase()<<" ++++++++++++++++"<< endl<<"||"<<endl<<"||  "<<cntrlfile_name <<endl;
  cout<<"||"<<endl<<"||"<<endl;
  cout <<"_______________________________________________________________"<<endl;
    adstring comma = adstring(","); 
    report << model_name<<" "<< endl<< endl;
    report << "Estimated annual F's " << endl;
    Fmort.initialize();
    for (k=1;k<=nfsh;k++)
      for (i=styr;i<=endyr;i++) 
        Fmort(i) += mean(F(k,i));
    report << Fmort<<endl;
    report << "Total mortality (Z)"<<endl;
    report << Z<<endl;
    report << "Estimated numbers of fish " << endl;
    for (i=styr;i<=endyr;i++) 
      report <<"       Year: "<< i << " "<< natage(i) << endl;
    report << endl<< "Estimated F mortality " << endl;
    for (k=1;k<=nfsh;k++)
    {
      report << "Fishery "<< k <<" : "<< endl ;
      for (i=styr;i<=endyr;i++) 
        report << "        Year: "<<i<<" "<<F(k,i)<<  " "<< endl;
    }

    report << endl<< "Observed survey values " << endl;
    for (k=1;k<=nsrv;k++)
    {
      int ii=1;
      report <<endl<< "Yr_Obs_Pred_Survey "<< k <<" : "<< endl ;
      for (i=styr;i<=endyr;i++)
      {
        if (ii<=yrs_srv(k).indexmax())
        {
          if (yrs_srv(k,ii)==i)
          {
            report << i<< " "<< obs_srv(k,ii) << " "<< pred_srv(k,i) <<endl;
            ii++;
          }
          else
            report << i<< " -1 "<< " "<< pred_srv(k,i)<<endl;
        }
        else
          report << i<< " -1 "<< " "<< pred_srv(k,i)<<endl;
      }
    }

    report << endl<< "Survey_Q:  "<<q_srv << endl;

    report << endl<< "Observed Prop " << endl;
    for (k=1;k<=nfsh;k++)
    {
      report << "ObsFishery "<< k <<" : "<< endl ;
      for (i=1;i<=nyrs_fsh_age(k);i++) 
        report << yrs_fsh_age(k,i)<< " "<< oac_fsh(k,i) << endl;
    }
    report << endl<< "Predicted prop  " << endl;
    for (k=1;k<=nfsh;k++)
    {
      report << "PredFishery "<< k <<" : "<< endl;
      for (i=1;i<=nyrs_fsh_age(k);i++) 
        report << yrs_fsh_age(k,i)<< " "<< eac_fsh(k,i) << endl;
    }
    report << endl<< "Observed prop Survey" << endl;
    for (k=1;k<=nsrv;k++)
    {
      report << "ObsSurvey "<<k<<" : "<<  endl;
      for (i=1;i<=nyrs_srv_age(k);i++) 
        report << yrs_srv_age(k,i)<< " "<< oac_srv(k,i) << endl;
    }
    report << endl<< "Predicted prop Survey" << endl;
    for (k=1;k<=nsrv;k++)
    {
      report << "PredSurvey "<<k<<" : "<<  endl;
      for (i=1;i<=nyrs_srv_age(k);i++) 
        report << yrs_srv_age(k,i)<< " "<< eac_srv(k,i) << endl;
    }
    report << endl<< "Observed catch biomass " << endl;
    report << catch_bio << endl;
    report << "predicted catch biomass " << endl;
    report << pred_catch << endl;

    report << endl<< "Estimated annual fishing mortality " << endl;
    for (k=1;k<=nfsh;k++)
      report << " Average_F_Fshry_"<<k<< " Full_selection_F_Fshry_"<<k;

    report << endl;
    for (i=styr;i<=endyr;i++)
    {
      report<< i<< " ";
      for (k=1;k<=nfsh;k++)
        report<< mean(F(k,i)) <<" "<< mean(F(k,i))*max(sel_fsh(k,i)) << " ";

      report<< endl;
    }
    report << endl<< "Selectivity" << endl;
    for (k=1;k<=nfsh;k++)
      for (i=styr;i<=endyr;i++)
        report << "Fishery "<< k <<"  "<< i<<" "<<sel_fsh(k,i) << endl;
    for (k=1;k<=nsrv;k++)
      for (i=styr;i<=endyr;i++)
        report << "Survey  "<< k <<"  "<< i<<" "<<sel_srv(k,i) << endl;

    report << endl<< "Stock Recruitment stuff "<< endl;
    for (i=styr_rec;i<=endyr;i++)
      if (active(log_Rzero))
        report << i<< " "<<Sp_Biom(i-rec_age)<< " "<< SRecruit(Sp_Biom(i-rec_age))<< " "<< mod_rec(i)<<endl;
      else 
        report << i<< " "<<Sp_Biom(i-rec_age)<< " "<< " 999" << " "<< mod_rec(i)<<endl;

    report << endl<< "Curve to plot "<< endl;
    report <<"stock Recruitment"<<endl;
    report <<"0 0 "<<endl;
    dvariable stock;
    for (i=1;i<=30;i++)
    {
      stock = double (i) * Bzero /25.;
      if (active(log_Rzero))
        report << stock <<" "<< SRecruit(stock)<<endl;
      else
        report << stock <<" 99 "<<endl;
    }

    report   << endl<<"Likelihood Components" <<endl;
    report   << "----------------------------------------- " <<endl;
    report   << "  catch_like  age_like_fsh sel_like_fsh surv_like age_like_srv sel_like_srv rec_like fpen post_priors_srvq post_priors residual total"<<endl;
    report   << " "<<obj_comps<<endl;

    obj_comps(11)= obj_fun - sum(obj_comps) ; // Residual 
    obj_comps(12)= obj_fun ;                  // Total
    report   <<"  catch_like       "<<setw(10)<<obj_comps(1) <<endl
             <<"  age_like_fsh     "<<setw(10)<<obj_comps(2) <<endl
             <<"  sel_like_fsh     "<<setw(10)<<obj_comps(3) <<endl
             <<"  surv_like        "<<setw(10)<<obj_comps(4) <<endl
             <<"  age_like_srv     "<<setw(10)<<obj_comps(5) <<endl
             <<"  sel_like_srv     "<<setw(10)<<obj_comps(6) <<endl
             <<"  rec_like         "<<setw(10)<<obj_comps(7) <<endl
             <<"  fpen             "<<setw(10)<<obj_comps(8) <<endl
             <<"  post_priors_srvq "<<setw(10)<<obj_comps(9) <<endl
             <<"  post_priors      "<<setw(10)<<obj_comps(10)<<endl
             <<"  residual         "<<setw(10)<<obj_comps(11)<<endl
             <<"  total            "<<setw(10)<<obj_comps(12)<<endl;

    report   << endl;
    report << "Fit to Catch Biomass "<<endl;
    report   << "-------------------------" <<endl;
    for (k=1;k<=nfsh;k++)
      report << "  Catch_like_Fshry_#"<< k <<"  "<< catch_like(k) <<endl;
    report   << endl;

    report << "Age likelihoods for fisheries :"<<endl;
    report   << "-------------------------" <<endl;
    for (k=1;k<=nfsh;k++)
      report << "  Age_like_Fshry_#"<< k <<"  "<< age_like_fsh(k) <<endl;
    report   << endl;

    report   << "Selectivity penalties for fisheries :"<<endl;
    report   << "-------------------------" <<endl;
    report   << "  Fishery Curvature_Age Change_Time Dome_Shaped"<<endl;
    for (k=1;k<=nfsh;k++)
      report << "  Sel_Fshry_#"<< k <<"  "<< sel_like_fsh(k,1) <<" "<<sel_like_fsh(k,2)<<" "<<sel_like_fsh(k,3)<< endl;
    report   << endl;
  
    report   << "survey Likelihood(s) " <<endl;
    report   << "-------------------------" <<endl;
    for (k=1;k<=nsrv;k++)
      report << "  Survey_Index_#"<< k <<"  " << surv_like(k)<<endl;
    report   << endl;

    report << setw(10)<< setfixed() << setprecision(5) <<endl;
    report   << "Age likelihoods for surveys :"<<endl;
    report   << "-------------------------" <<endl;
    for (k=1;k<=nsrv;k++)
      report << "  Age_Survey_#"<< k <<"  " << age_like_srv(k)<<endl;
    report   << endl;

    report   << "Selectivity penalties for surveys :"<<endl;
    report   << "-------------------------" <<endl;
    report   << "  Survey Curvature_Age Change_Time Dome_Shaped"<<endl;
    for (k=1;k<=nsrv;k++)
      report << "  Sel_Survey_#"<< k <<"  "<< sel_like_srv(k,1) <<" "<<sel_like_srv(k,2)<<" "<<sel_like_srv(k,3)<< endl;
    report   << endl;

    report << setw(10)<< setfixed() << setprecision(5) <<endl;
    report   << "Recruitment penalties: " <<rec_like<<endl;
    report   << "-------------------------" <<endl;
    report   << "  (sigmar)            " <<sigmar<<endl;
    report   << "  S-R_Curve           " <<rec_like(1)<< endl;
    report   << "  Regularity          " <<rec_like(2)<< endl;
    report   << "  Future_Recruits     " <<rec_like(3)<< endl;
    report   << endl;

    report   << "F penalties:          " <<endl;
    report   << "-------------------------" <<endl;
    report   << "  Avg_F               " <<fpen(1) <<endl;
    report   << "  Effort_Variability  " <<fpen(2) <<endl;
    report   << endl;

    report   << "Contribution of Priors:"<<endl;
    report   << "-------------------------" <<endl;
    report   << "Source                ";
    report   <<           " Posterior";
    report   <<           " Param_Val";
    report   <<           " Prior_Val";
    report   <<           "  CV_Prior"<<endl;
  // (*ad_printf)("f = %lf\n",value(f));
    for (k=1;k<=nsrv;k++)
    {
      report << "Q_Survey_#"<< k <<"           "
             << setw(10)<<post_priors_srvq(k) 
             << setw(10)<< q_srv(k)
             << setw(10)<< qprior(k)
             << setw(10)<< cvqprior(k)<<endl;

      report << "Q_power_Survey_#"<< k <<"           "
             << setw(10)<<post_priors_srvq(k) 
             << setw(10)<< q_power_srv(k)
             << setw(10)<< q_power_prior(k)
             << setw(10)<< cvq_power_prior(k)<<endl;
    }

    // writerep(post_priors(1),repstring);
    // cout <<repstring<<endl;
    report   << "Natural_Mortality     "
             << setw(10)<< post_priors(1)
             << setw(10)<< M
             << setw(10)<< natmortprior
             << setw(10)<< cvnatmortprior <<endl;
    report   << "Steepness             "
             << setw(10)<< post_priors(2)
             << setw(10)<< steepness
             << setw(10)<< steepnessprior
             << setw(10)<< cvsteepnessprior <<endl;
    report   << "SigmaR                "
             << setw(10)<< post_priors(3)
             << setw(10)<< sigmar
             << setw(10)<< sigmarprior
             << setw(10)<< cvsigmarprior <<endl;
    report   << endl;
    report<<"Num_parameters_Estimated "<<initial_params::nvarcalc()<<endl;
    
  report <<cntrlfile_name<<endl;
  report <<datafile_name<<endl;
  report <<model_name<<endl;
  if (SrType==2) 
    report<< "Beverton-Holt" <<endl;
  else
    report<< "Ricker" <<endl;
  report<<"Steepnessprior,_CV,_phase: " <<steepnessprior<<" "<<
    cvsteepnessprior<<" "<<
    phase_srec<<" "<< endl;

  report<<"sigmarprior,_CV,_phase: " <<sigmarprior<<" "<<  cvsigmarprior <<" "<<phase_sigmar<<endl;

  report<<"Rec_estimated_in_styr_endyr: " <<styr_rec    <<" "<<endyr        <<" "<<endl;
  report<<"SR_Curve_fit__in_styr_endyr: " <<styr_rec_est<<" "<<endyr_rec_est<<" "<<endl;
  report<<"Model_styr_endyr:            " <<styr        <<" "<<endyr        <<" "<<endl;

  report<<"M_prior,_CV,_phase "<< natmortprior<< " "<< cvnatmortprior<<" "<<phase_M<<endl;
  report<<"qprior,_CV,_phase " <<qprior<<" "<<cvqprior<<" "<< phase_q<<endl;
  report<<"q_power_prior,_CV,_phase " <<q_power_prior<<" "<<cvq_power_prior<<" "<< phase_q_power<<endl;

  report<<"cv_catchbiomass: " <<cv_catchbiomass<<" "<<endl;
  report<<"Projection_years "<< nproj_yrs<<endl;
  for (k=1;k<=nfsh;k++)
    report << "Fsh_sel_opt_fish: "<<k<<" "<<fsh_sel_opt(k)<<" "<<sel_change_in_fsh(k)<<endl;
  for (k=1;k<=nsrv;k++)
    report<<"Survey_Sel_Opt_Survey: " <<k<<" "<<(srv_sel_opt(k))<<endl;
    
  report <<"Phase_survey_Sel_Coffs: "<<phase_selcoff_srv<<endl; 
  report <<"Fshry_Selages: " << nselages_in_fsh  <<endl;
  report <<"Survy_Selages: " << nselages_in_srv <<endl;
  report << "Phase_for_age-spec_fishery "<<phase_selcoff_fsh<<endl;
  report << "Phase_for_logistic_fishery "<<phase_logist_fsh<<endl;
  report << "Phase_for_dble_logistic_fishery "<<phase_dlogist_fsh<<endl;
  report << "Phase_for_age-spec_survey  "<<phase_selcoff_srv<<endl;
  report << "Phase_for_logistic_survey  "<<phase_logist_srv<<endl;
  report << "Phase_for_dble_logistic_srvy "<<phase_dlogist_srv<<endl;

  for (k=1; k<=nfsh;k++)
  {
    report <<"Number_of_select_changes_fishery: "<<k<<" "<<n_sel_ch_fsh(k)<<endl;
    report<<"Yrs_fsh_sel_change: "<<yrs_sel_ch_fsh(k)<<endl;
    report << "sel_change_in: "<<sel_change_in_fsh(k) << endl;
  }
  for (k=1; k<=nsrv;k++)
  {
    report <<"Number_of_select_changes_survey: "<<k<<" "<<n_sel_ch_srv(k)<<endl;
    report<<"Yrs_srv_sel_change: "<<yrs_sel_ch_srv(k)<<endl;
    report << "sel_change_in: "<<sel_change_in_srv(k) << endl;
  }

FUNCTION write_msy_out
  ofstream msyout("msyout.dat");
  msyout << " # Natural Mortality       " <<endl;
  for (j=1;j<=nages;j++) 
    msyout <<M <<" ";
  msyout <<endl;
  msyout << spawnmo<< "  # Spawnmo                   " <<endl;
  msyout <<"# Wt spawn"<<endl<< wt_pop<< endl;
  msyout <<"# Wt fish"<<endl;
  for (k=1;k<=nfsh;k++) 
    msyout <<wt_fsh(k,endyr)<< " ";
  msyout <<endl;
  msyout <<"# Maturity"<<endl<< maturity<< endl;
  msyout <<"# selectivity"<<endl;
  for (k=1;k<=nfsh;k++) 
    msyout<< sel_fsh(k,endyr) <<" ";
  msyout<< endl;
  msyout<<"Srec_Option "<<SrType<< endl;
  msyout<<"Alpha "<<alpha<< endl;
  msyout<<"beta "<<beta<< endl;
  msyout<<"steepness "<<steepness<< endl;
  msyout<<"Bzero "<<Bzero<< endl;
  msyout<<"Rzero "<<Rzero<< endl;

FUNCTION write_projout
// Function to write out data file for projection model....
  ofstream projout( projfile_name );
  
  projout <<"# "<<model_name <<" "<< projfile_name<<endl;
  projout <<"123  # seed"<<endl;
  // Flag to tell if this is a SSL species                 
  projout << "1 # Flag to tell if this is a SSL forage species                 "<<endl;
  projout << "0 # Flag to Dorn's version of a constant buffer                  "<<endl;
  // Flag to solve for F in first year or not 0==don't solve
  projout<< " 1 # Flag to solve for F in first year or not 0==don't solve"<<endl;
  // Flag to use 2nd-year catch/TAC
  projout<< "0 # Flag to use 2nd-year catch/TAC"<<endl;
  projout << nfsh<<"   # Number of fisheries"<<endl;
  projout <<"14   # Number of projection years"<<endl;
  projout <<"1000 # Number of simulations"<<endl;
  projout <<endyr<< " # Begin year of projection" <<endl;
  projout <<nages<< " # Number of ages" <<endl;
  for (j=1;j<=nages;j++) 
    projout <<M <<" ";
  projout << " # Natural Mortality       " <<endl;
  double sumtmp;
  sumtmp = 0.;
  for (k=1;k<=nfsh;k++) 
    sumtmp += catch_bio(k,endyr);
  projout << sumtmp<< " # TAC in current year (assumed catch) " <<endl;
  projout << sumtmp<< " # TAC in current year+1 (assumed catch) " <<endl;

  for (k=1;k<=nfsh;k++) 
    projout <<  F(k,endyr)/mean((F(k,endyr)))<<" ";
   //  + fmort_dev(k,endyr)) /Fmort(endyr)<<" ";

  projout << "   # Fratio                  " <<endl;

  projout << mean(Fmort(endyr-4,endyr))<<"  # average f                  " <<endl;
  projout << " 1  # author f                  " <<endl;
  projout << spawnmo<< "  # Spawnmo                   " <<endl;
  projout <<"# Wt spawn"<<endl<< wt_pop<< endl;
  projout <<"# Wt fish"<<endl;
  for (k=1;k<=nfsh;k++) 
    projout <<wt_fsh(k,endyr)<< " ";
  projout <<endl;
  projout <<"# Maturity"<<endl<< maturity<< endl;
  projout <<"# selectivity"<<endl;
  for (k=1;k<=nfsh;k++) 
    projout<< sel_fsh(k,endyr) <<" ";
  projout<< endl;
  projout <<"# natage"<<endl<< natage(endyr) << endl;
  projout <<"#_N_recruitment_years (not including last 3 estimates)"<<endl<<endyr-(1977+rec_age+3) << endl;
  projout <<"#_Recruitment_start_at_1977_yearclass=1978_for_age_1_recruits"<<yy(1977+rec_age,endyr-3)<<endl<<mod_rec(1977+rec_age,endyr-3)<< endl;

FUNCTION write_proj
 ofstream newproj("proj.dat");
// Function to write out data file for new Ianelli 2005 projection model....
 newproj <<"#Species name here:"<<endl;
 newproj <<model_name+"_"+datafile_name<<endl;
 newproj <<"#SSL Species?"<<endl;
 newproj <<"1"<<endl;
 newproj <<"#Constant buffer of Dorn?"<<endl;
 newproj <<"0"<<endl;
 newproj <<"#Number of fisheries?"<<endl;
 newproj <<"1"<<endl;
 newproj <<"#Number of sexes?"<<endl;
 newproj <<"1"<<endl;
 newproj <<"#5year_Average_F(endyr-4,endyr_as_estimated_by_ADmodel)"<<endl;
 // Need to correct for maxf standardization 
 dvector seltmp = value(sel_fsh(1,endyr)) 
                 +value(sel_fsh(1,endyr-1))  
                 +value(sel_fsh(1,endyr-2))  
                 +value(sel_fsh(1,endyr-3))  
                 +value(sel_fsh(1,endyr-4));  
 newproj << mean(Fmort(endyr-4,endyr))<<endl;
 newproj <<"#_Author_F_as_fraction_F_40%"<<endl;
 newproj <<"1"<<endl;
 newproj <<"#ABC SPR" <<endl;
 newproj <<"0.4"<<endl;
 newproj <<"#MSY SPR" <<endl;
 newproj <<"0.35"<<endl;
 newproj <<"#_Spawn_month"<<endl;
 newproj << spmo_frac*12+1<<endl;
 newproj <<"#_Number_of_ages"<<endl;
 newproj <<nages<<endl;
 newproj <<"#_F_ratio(must_sum_to_one_only_one_fishery)"<<endl;
 newproj <<"1"<<endl;
 newproj <<"#_Natural_Mortality" << aa << endl;
   for (j=1;j<=nages;j++) newproj <<natmort<<" "; newproj<<endl;
 newproj <<"#_Maturity_divided_by_2(projection_program_uses_to_get_female_spawning_biomass_if_divide_by_2"<<aa<<endl<<2.*maturity<< endl;
 newproj <<"#_Wt_at_age_spawners"<<aa<<endl<<wt_pop<< endl;
 newproj <<"#_Wt_at_age_fishery" <<aa<<endl<<wt_fsh(1,endyr) << endl;
 newproj <<"#" <<endl;

 newproj <<"#_Selectivity_fishery_scaled_to_max_at_one"<<aa<<endl;
 seltmp = value(sel_fsh(1,endyr)) +value(sel_fsh(1,endyr-1))  +value(sel_fsh(1,endyr-2));  
 newproj << seltmp/max(seltmp)<<endl;
 newproj <<"#_Numbers_at_age_end_year"<<aa<<endl<<natage(endyr)<< endl;
 newproj <<"#_N_recruitment_years (not including last estimate)"<<endl<<endyr-(1977+rec_age) << endl;

 newproj <<"#_Recruitment_start_at_1977_yearclass=1978_for_age_1_recruits"<<yy(1977+rec_age,endyr-1)
         <<endl<<mod_rec(1977+rec_age,endyr-1)<< endl;

 newproj <<"#_Spawning biomass "<<endl<<Sp_Biom(styr-rec_age,endyr-rec_age)/1000<< endl;
 newproj.close();
 
RUNTIME_SECTION
  convergence_criteria 1.e-1,,1.e-01,1.e-03,1e-5,1e-5
  maximum_function_evaluations 1500

TOP_OF_MAIN_SECTION
  gradient_structure::set_MAX_NVAR_OFFSET(1500);
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(1500);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(1000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(10000000);
  arrmblsize=500000000;

FINAL_SECTION
  // Calc_Dependent_Vars();
  Write_R();
  write_proj();
  write_projout();
  write_msy_out();
  Profile_F();
  
GLOBALS_SECTION
  #include <admodel.h>
	/**
	\def log_input(object)
	Prints name and value of \a object on ADMB report %ofstream file.
	*/
	#undef log_input
	#define log_input(object) write_input_log << "# " #object "\n" << object << endl;
	#undef log_param
  // #define log_param(object) for(int i=0;i<initial_params::num_initial_params;i++) {if(withinbound(0,(initial_params::varsptr[i])->phase_start, initial_params::current_phase)) { int sc= (initial_params::varsptr[i])->size_count(); if (sc>0) { write_input_log << "# " << initial_params::varsptr[i] ->label() << "\n" << object<<endl; } }}
	#define log_param(object) if (active(object)) write_input_log << "# " #object "\n" << object << endl;
  ofstream write_input_log("write_input_log.log");
 // void get_sel_changes(int& k);
  adstring_array fshname;
  adstring_array srvname;
  adstring simname;
  adstring model_name;
  adstring projfile_name;
  adstring datafile_name;
  adstring cntrlfile_name;
  adstring tmpstring;
  adstring repstring;
  adstring version_info;

 typedef vcubic_spline_function * pvcubic_spline_function;
 class vcubic_spline_function_array
 {
 public:
   int indexmin(void) {return index_min;}
   int indexmax(void) {return index_max;}
   vcubic_spline_function_array(int,int,const dvector & x, 
           const dvar_matrix& _y);
   ~vcubic_spline_function_array();
   vcubic_spline_function & operator () (int i);
   vcubic_spline_function_array(const dvector & x);
   // so that this will fail define it if you need it.
   vcubic_spline_function_array(const vcubic_spline_function_array&);
   dvar_matrix operator( ) (const dvector & v);
   
 private:
   vcubic_spline_function ** ptr;
   int index_min;
   int index_max;
 };

  vcubic_spline_function & vcubic_spline_function_array::operator () (int i)
  {
    if (i<indexmin() || i> indexmax())
    {
      cerr << "index out of range in function"
           " vcubic_spline_function & operator () (int i)"
       << endl;
      ad_exit(1);
    }
    return *(ptr[i]);
  }
    
  dvar_matrix vcubic_spline_function_array::operator( ) (const dvector & v)
  {
    int mmin=indexmin();
    int mmax=indexmax();
    dvar_matrix tmp(mmin,mmax,v.indexmin(),v.indexmax());
    for (int i=mmin;i<=mmax;i++)
    {
      tmp(i)=(*this)(i)(v);   
    }
    return tmp;
  }

  vcubic_spline_function_array::vcubic_spline_function_array(int mmin,int mmax,
  const dvector & x, const dvar_matrix& y)
  {
    index_min=mmin;
    index_max=mmax;
    int n=mmax-mmin+1;
    ptr = new pvcubic_spline_function[n];
    ptr-=mmin;
    for (int i=mmin;i<=mmax;i++)
    {
      ptr[i]= new  vcubic_spline_function(x,y(i));
    }
  }
  vcubic_spline_function_array::~vcubic_spline_function_array()
  {
    int mmin=indexmin();
    int mmax=indexmax();
  
    for (int i=mmin;i<=mmax;i++)
    {
      delete ptr[i];
    }
    ptr+=mmin;
    delete ptr;
    ptr=0;
  }

FUNCTION dvariable get_spr_rates(double spr_percent)
  RETURN_ARRAYS_INCREMENT();
  dvar_matrix sel_tmp(1,nages,1,nfsh);
  sel_tmp.initialize();
  for (k=1;k<=nfsh;k++)
    for (j=1;j<=nages;j++)
      sel_tmp(j,k) = sel_fsh(k,endyr,j); // NOTE uses last-year of fishery selectivity for projections.
  dvariable sumF=0.;
  for (k=1;k<=nfsh;k++)
  {
    Fratio(k) = sum(F(k,endyr)) ;
    sumF += Fratio(k) ;
  }
  Fratio /= sumF;
  double df=1.e-3;
  dvariable F1 ;
  F1.initialize();
  F1 = .8*natmortprior;
  dvariable F2;
  dvariable F3;
  dvariable yld1;
  dvariable yld2;
  dvariable yld3;
  dvariable dyld;
  dvariable dyldp;
  // Newton Raphson stuff to go here
  for (int ii=1;ii<=6;ii++)
  {
    F2     = F1 + df;
    F3     = F1 - df;
    yld1   = -1000*square(log(spr_percent/spr_ratio(F1, sel_tmp)));
    yld2   = -1000*square(log(spr_percent/spr_ratio(F2, sel_tmp)));
    yld3   = -1000*square(log(spr_percent/spr_ratio(F3, sel_tmp)));
    dyld   = (yld2 - yld3)/(2*df);                          // First derivative (to find the root of this)
    dyldp  = (yld3-(2*yld1)+yld2)/(df*df);  // Newton-Raphson approximation for second derivitive
    F1    -= dyld/dyldp;
  }
  RETURN_ARRAYS_DECREMENT();
  return(F1);

FUNCTION dvariable spr_ratio(dvariable trial_F,dvar_matrix sel_tmp)
  RETURN_ARRAYS_INCREMENT();
  dvariable SBtmp;
  dvar_vector Ntmp(1,nages);
  dvar_vector srvtmp(1,nages);
  SBtmp.initialize();
  Ntmp.initialize();
  srvtmp.initialize();
  dvar_matrix Ftmp(1,nages,1,nfsh);
  Ftmp = sel_tmp;
  for (j=1;j<=nages;j++) 
  {
    Ftmp(j) = elem_prod(Ftmp(j), trial_F * Fratio);
    srvtmp(j)  = mfexp(-sum(Ftmp(j)) - natmort);
  }
  Ntmp(1)=1.;
  j=1;
  SBtmp  += Ntmp(j)*wt_mature(j)*pow(srvtmp(j),spmo_frac);
  for (j=2;j<nages;j++)
  {
    Ntmp(j) = Ntmp(j-1)*srvtmp(j-1);
    SBtmp  += Ntmp(j)*wt_mature(j)*pow(srvtmp(j),spmo_frac);
  }
  Ntmp(nages)=Ntmp(nages-1)*srvtmp(nages-1)/(1.-srvtmp(nages));

  SBtmp  += Ntmp(nages)*wt_mature(nages)*pow(srvtmp(nages),spmo_frac);
  RETURN_ARRAYS_DECREMENT();
  return(SBtmp/phizero);

FUNCTION dvariable spr_unfished()
  RETURN_ARRAYS_INCREMENT();
  dvariable Ntmp;
  dvariable SBtmp;
  SBtmp.initialize();
  Ntmp = 1.;
  for (j=1;j<nages;j++)
  {
    SBtmp += Ntmp*wt_mature(j)*exp(-spmo_frac * natmort);
    Ntmp  *= mfexp( -natmort);
  }
  Ntmp    /= (1.-exp(-natmort));
  SBtmp += Ntmp*wt_mature(j)*exp(-spmo_frac * natmort);
  RETURN_ARRAYS_DECREMENT();
  return(SBtmp);

FUNCTION compute_spr_rates
  //Compute SPR Rates and add them to the likelihood for Females 
  dvariable sumF=0.;
  for (k=1;k<=nfsh;k++)
  {
    Fratio(k) = sum(F(k,endyr)) ;
    sumF += Fratio(k) ;
  }
  Fratio /= sumF;

  F35_est = get_spr_rates(.35);
  F50_est = get_spr_rates(.50);
  F40_est = get_spr_rates(.40);

  for (k=1;k<=nfsh;k++)
  {
    F50(k) = F50_est * (Fratio(k));
    F40(k) = F40_est * (Fratio(k));
    F35(k) = F35_est * (Fratio(k));
  }

  /* FUNCTION get_agematrix
  // B3+($B$5-$B$4)/(COUNT($B$2:$L$2)-1)
  for (i=1;i<=nages;i++)
    {
     sd_age(i)=cv_age(i)*len_age(i);
     var_age(i)=sd_age(i)*sd_age(i);
     for (int j=1;j<=nlenbins;j++)
       {
        diff = len(j)-len_age(i);
        trans(i,j)=2/sd_age(i)*exp(-diff*diff/(2*var_age(i)));
       }
      trans(i)=trans(i)/sum(trans(i));
    }
    */ 


FUNCTION void writerep(dvariable& tmp,adstring& tmpstring)
  cout <<tmpstring<<endl<<endl;
  tmpstring = printf("3.5%f",value(tmp));

FUNCTION void do_Newton_Raphson_for_mortality(dvariable hrate)
  dvariable Fold ;
  Fold = hrate;
  for (int ii=1;ii<=4;ii++)
  {
      dvariable ZZ = Fold + natmort;
      dvariable SS = mfexp(-ZZ);
      dvariable AA = Fold * (1. - SS);
      dvariable BB = ZZ;
      dvariable CC = 1. + (Fold - 1) * SS;
      dvariable dd = 1.;
      dvariable FX = AA / BB - hrate;
      dvariable FPX = (BB * CC - AA * dd) / (BB * BB);
      Fnew = Fold - FX / FPX;
      Fold = Fnew;
  }

FUNCTION dvariable SolveF2(const int& iyr, const dvar_vector& N_tmp, const double&  TACin)
  RETURN_ARRAYS_INCREMENT();
  dvariable dd = 10.;
  dvariable cc; 
  dvar_matrix Fratsel(1,nfsh,1,nages);
  dvar_vector Z_tmp(1,nages) ;
  dvar_vector S_tmp(1,nages) ;
  dvar_vector Ftottmp(1,nages);
  dvariable btmp =  N_tmp * elem_prod(sel_fsh(1,iyr),wt_pop);
  dvariable ftmp;
  ftmp = TACin/btmp;
    for (k=1;k<=nfsh;k++)
      Fratsel(k) = Fratio(k)*sel_fsh(k,iyr);
    for (int ii=1;ii<=5;ii++)
    {
      Ftottmp.initialize();
      for (k=1;k<=nfsh;k++)
        Ftottmp += ftmp*Fratsel(k);
  
      Z_tmp = Ftottmp  + M; 
      S_tmp = mfexp( -Z_tmp );
      cc = 0.0;
      for (k=1;k<=nfsh;k++)
        cc += wt_fsh(k,endyr) * elem_prod(elem_div(ftmp*Fratsel(k),  Z_tmp),elem_prod(1.-S_tmp,N_tmp)); // Catch equation (vectors)
  
      dd = cc / TACin - 1.;
      if (dd<0.) dd *= -1.;
      ftmp += (TACin-cc) / btmp;
    }
  RETURN_ARRAYS_DECREMENT();
  return(ftmp);

FUNCTION dvar_vector SolveF2(const int& iyr, const dvector&  Catch)
  // Returns vector of F's (given year) by fleet
  // Requires: N and fleet specific wts & selectivities at age, catch 
  // iterate to get Z's right
  RETURN_ARRAYS_INCREMENT();
  dvariable dd = 10.;
  dvariable cc; 
  dvar_matrix  seltmp(1,nfsh,1,nages);
  dvar_matrix  wt_tmp(1,nfsh,1,nages);
  dvar_matrix Fratsel(1,nfsh,1,nages);
  dvar_vector N_tmp = natage(iyr);
  dvar_vector Z_tmp(1,nages) ;
  dvar_vector S_tmp(1,nages) ;
  dvar_vector Ftottmp(1,nages);
  dvar_vector Frat(1,nfsh);
  dvar_vector btmp(1,nfsh);
  dvar_vector ftmp(1,nfsh);
  dvar_vector hrate(1,nfsh);
  btmp.initialize(); 
  // Initial guess for Fratio
  for (k=1;k<=nfsh;k++)
  {
    seltmp(k)= sel_fsh(k,iyr); // Selectivity
    wt_tmp(k)= wt_fsh(k,iyr); // 
    btmp(k)  =  N_tmp * elem_prod(seltmp(k),wt_tmp(k));
    hrate(k) = Catch(k)/btmp(k);
    Frat(k)  = Catch(k)/sum(Catch);
    Fratsel(k) = Frat(k)*seltmp(k);
    ftmp(k) = 1.1*(1.- posfun(1.-hrate(k),.10,fpen(4)));
  }
  // Initial fleet-specific F
  // iterate to balance effect of multiple fisheries...........
  for (int kk=1;kk<=nfsh;kk++) 
  {
    for (k=1;k<=nfsh;k++)
    {
      if (hrate(k) <.9999) 
      {
        for (int ii=1;ii<=6;ii++)
        {
          Ftottmp.initialize();
          Ftottmp   = ftmp*Fratsel;
          Z_tmp = Ftottmp  + M; 
          S_tmp = mfexp( -Z_tmp );
          cc = wt_tmp(k) * elem_prod(elem_div(ftmp(k)*Fratsel(k),  Z_tmp),elem_prod(1.-S_tmp,N_tmp)); // Catch equation (vectors)
          ftmp(k) += ( Catch(k)-cc ) / btmp(k);
        }
        Frat(k)    = ftmp(k)/sum(ftmp);
        Fratsel(k) = Frat(k)*seltmp(k);
      }
    }
  }
  RETURN_ARRAYS_DECREMENT();
  return(ftmp);


FUNCTION Write_Datafile
  int nsims;
  // get the number of simulated datasets to create...
  ifstream sim_in("nsims.dat"); sim_in >> nsims; sim_in.close();
  char buffer [33];
  // compute the autocorrelation term for residuals of fit to indices...
  for (k=1;k<=nsrv;k++)
    ac(k) = get_AC(k);
  for (int isim=1;isim<=nsims;isim++)
  {
    // Create the name of the simulated dataset
    simname = "sim_"+ adstring(itoa(isim,buffer,10)) + ".dat";
    // Open the simulated dataset for writing
    ofstream simdat(simname);
    simdat << "# first year" <<endl;
    simdat << styr <<endl;
    simdat << "# Last  year" <<endl;
    simdat << endyr <<endl;
    simdat << "# age recruit" <<endl;
    simdat << rec_age <<endl;
    simdat << "# oldest age" <<endl;
    simdat << oldest_age <<endl;
    simdat << "# Number of fisheries " <<endl;
    simdat << nfsh <<endl;                                   
    simdat << fshnameread <<endl;                                   
    simdat << "# Catch biomass by fishery " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k) <<" " << k <<endl;
      simdat << catch_bio(k) <<endl;
    }
    simdat << "# Catch biomass uncertainty by fishery (std errors)" <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k) <<" " << k <<endl;
      simdat << catch_bio_sd(k) <<endl;   
    }
    simdat << "# number of years for fishery age data " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k)<< " " << k <<endl;
      simdat << nyrs_fsh_age(k) <<endl;
    }
    simdat << "# years for fishery age data " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k)<< " " << k <<endl;
      simdat << yrs_fsh_age(k)  <<endl;
    }
    simdat << "# sample sizes for fishery age data " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k)<< " " << k <<endl;
      simdat << n_sample_fsh_age(k)  <<endl;    
    }
    simdat << "# Observed age compositions for fishery" <<endl;
    for (k=1;k<=nfsh;k++)
    {
      dvector p(1,nages);
      double Ctmp; // total catch
      dvector freq(1,nages);
      simdat << "# " << fshname(k) <<endl;
      for (i=1;i<=nyrs_fsh_age(k);i++)
      {
        int iyr = yrs_fsh_age(k,i);
        // Add noise here
        freq.initialize();
        ivector bin(1,n_sample_fsh_age(k,i));
        p = value(catage(k,iyr));
        p /= sum(p);
        bin.fill_multinomial(rng,p); // fill a vector v
        for (int j=1;j<=n_sample_fsh_age(k,i);j++)
          freq(bin(j))++;
        // Apply ageing error to samples..............
        p = age_err *freq/sum(freq); 
        cout << p  <<endl;
        simdat << p  <<endl;
        // Compute total catch given this sample size
        Ctmp = catch_bio(k,iyr) / (p*wt_fsh(k,iyr)); 
        // Simulated catage = proportion sampled
        sim_catage(k,i) = p * Ctmp;
      }

    }
    simdat << "# Annual wt-at-age for fishery" <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k)<< " " << (k) <<endl;
      // Add noise here
      simdat << wt_fsh(k)  <<endl;  
    }
    simdat << "# number of indices" <<endl;
    simdat << nsrv <<endl;                                   
    simdat << srvnameread <<endl;                                   
    simdat << "# Number of years of index values (annual)" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " << srvname(k) <<endl;
      simdat << nyrs_srv(k)  <<endl;                   
    }
    simdat << "# Years of index values (annual)" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " << srvname(k) <<endl;
      simdat << yrs_srv(k)  <<endl;         
    }
    simdat << "# Month that index occurs "<<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " << srvname(k) <<endl;
      simdat << mo_srv(k) <<endl;
    }
    simdat << "# values for indices (annual)"<<endl;
    dmatrix new_srv(1,nsrv,1,nyrs_srv);
    new_srv.initialize();
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< " " << k <<endl;
      // Add noise here
      dvector ran_srv_vect(1,nyrs_srv(k));
      // fill vector with unit normal RVs
      ran_srv_vect.fill_randn(rng);
      // do first year uncorrelated
      i=1;
      int iyr=yrs_srv(k,i);
      corr_dev(k)  = ran_srv_vect;
      new_srv(k,i) = mfexp(corr_dev(k,i) * obs_lse_srv(k,i) ) * 
                     value(elem_prod(wt_srv(k,iyr),elem_prod(pow(S(iyr),srv_month_frac(k)), natage(iyr)))*
                     q_srv(k)*sel_srv(k,iyr)); 
      // do next years correlated with previous
      for (i=2;i<=nyrs_srv(k);i++)
      {
        iyr=yrs_srv(k,i);
        corr_dev(k,i) = ac(k) * corr_dev(k,i-1) + sqrt(1.-square(ac(k))) * corr_dev(k,i);
        new_srv(k,i) = mfexp(corr_dev(k,i) * obs_lse_srv(k,i) ) * 
                        value(elem_prod(wt_srv(k,iyr),elem_prod(pow(S(iyr),srv_month_frac(k)), 
                        natage(iyr))) * q_srv(k)*sel_srv(k,iyr)); 
      }
      simdat << new_srv(k)      <<endl;
    }
    simdat << "# standard errors for indices (by year) " <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< " " << k <<endl;
      simdat << obs_se_srv(k)  <<endl;
    }
    simdat << "# Number of years of age data available for index" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< " " << k <<endl;
      simdat << nyrs_srv_age(k)  <<endl;
    }
    simdat << "# Years of index values (annual)" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< endl;
      simdat << yrs_srv_age(k)  <<endl;
    }
    simdat << "# Sample sizes for age data from indices" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< endl;
      simdat << n_sample_srv_age(k)  <<endl;
    }
    simdat << "# values of proportions at age in index" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< endl;
      dvector p(1,nages);
      dvector freq(1,nages);
      for (i=1;i<=nyrs_srv_age(k);i++)
      {
        int iyr = yrs_srv_age(k,i);
        // Add noise here
        freq.initialize();
        ivector bin(1,n_sample_srv_age(k,i));
        p = age_err * value(elem_prod( elem_prod(pow(S(iyr),srv_month_frac(k)), natage(iyr))*q_srv(k) , sel_srv(k,iyr))); 
        p /= sum(p);
        // fill vector with multinomial samples
        bin.fill_multinomial(rng,p); // fill a vector v
        for (int j=1;j<=n_sample_srv_age(k,i);j++)
          freq(bin(j))++;
        simdat << "# " <<srvname(k)<< " year: "<< iyr<< endl;
        simdat << freq/sum(freq) <<endl;
      }
    }
    simdat << "# Mean wts at age for indices" <<endl;
    for (k=1;k<=nsrv;k++)
    {
      simdat << "# " <<srvname(k)<< endl;
      // Could add noise here
      simdat <<  wt_srv(k)  <<endl;
    }
  
    simdat << "# Population mean wt at age" <<endl;
    simdat << wt_pop <<endl;
  
    simdat << "# Population maturity at age" <<endl;
    simdat << maturity  <<endl;
  
    simdat << "# Peak spawning month" <<endl;
    simdat << spawnmo <<endl;
  
    simdat << "# ageing error " <<endl;
    simdat << age_err <<endl;

    simdat <<endl<<endl<<"Additional output"<<endl;
    simdat << "# Fishery_Effort " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      dvector ran_fsh_vect(styr,endyr);
      // fill vector with unit normal RVs
      ran_fsh_vect.fill_randn(rng);
      // Sigma on effort is ~15% white noise (add red noise later)
      ran_fsh_vect *= 0.15; 
      dvector avail_biom(styr,endyr);
      for (i=styr;i<=endyr;i++)
      {
        avail_biom(i) = wt_fsh(k,i)*value(elem_prod(natage(i),sel_fsh(k,i))); 
      }
      act_eff(k) = elem_prod(exp(ran_fsh_vect), (elem_div(catch_bio(k), avail_biom)) );
      // Normalize effort
      act_eff(k) /= mean(act_eff(k));
      for (i=styr;i<=endyr;i++)
        simdat<<fshname(k)<<" "<<i<<" "<<act_eff(k,i) <<endl;
    }
    simdat << "# Fishery catch-at-age " <<endl;
    for (k=1;k<=nfsh;k++)
    {
      simdat << "# " <<fshname(k)<< " " << k <<endl;
      simdat << "Fishery Year "<<age_vector << endl;
      for (i=1;i<=nyrs_fsh_age(k);i++)
        simdat<<fshname(k)<<" "<<yrs_fsh_age(k,i)<<" "<<sim_catage(k,i) <<endl;
    }
  }
  exit(1);
  // End of simulating datasets...................
  

FUNCTION Write_R
  ofstream R_report("For_R.dat");
  R_report<<"$Yr"<<endl; for (i=styr;i<=endyr;i++) R_report<<i<<" "; R_report<<endl;

  R_report<<"$TotF"<<endl << Ftot<<endl;

  R_report<<"$TotBiom_NoFish"<<endl; for (i=styr;i<=endyr;i++) 
  {
    double lb=value(totbiom_NoFish(i)/exp(2.*sqrt(log(1+square(totbiom_NoFish.sd(i))/square(totbiom_NoFish(i))))));
    double ub=value(totbiom_NoFish(i)*exp(2.*sqrt(log(1+square(totbiom_NoFish.sd(i))/square(totbiom_NoFish(i))))));
    R_report<<i<<" "<<totbiom_NoFish(i)<<" "<<totbiom_NoFish.sd(i)<<" "<<lb<<" "<<ub<<endl;
  }
  R_report<<"$SSB_NoFishR"<<endl; for (i=styr;i<=endyr;i++) 
  {
    double lb=value(Sp_Biom_NoFishR(i)/exp(2.*sqrt(log(1+square(Sp_Biom_NoFishR.sd(i))/square(Sp_Biom_NoFishR(i))))));
    double ub=value(Sp_Biom_NoFishR(i)*exp(2.*sqrt(log(1+square(Sp_Biom_NoFishR.sd(i))/square(Sp_Biom_NoFishR(i))))));
    R_report<<i<<" "<<Sp_Biom_NoFishR(i)<<" "<<Sp_Biom_NoFishR.sd(i)<<" "<<lb<<" "<<ub<<endl;
  }

  R_report<<"$TotBiom"<<endl; 
  for (i=styr;i<=endyr;i++) 
  {
    double lb=value(totbiom(i)/exp(2.*sqrt(log(1+square(totbiom.sd(i))/square(totbiom(i))))));
    double ub=value(totbiom(i)*exp(2.*sqrt(log(1+square(totbiom.sd(i))/square(totbiom(i))))));
    R_report<<i<<" "<<totbiom(i)<<" "<<totbiom.sd(i)<<" "<<lb<<" "<<ub<<endl;
  }

  for (int k=1;k<=5;k++){
    R_report<<"$SSB_fut_"<<k<<endl; 
    for (i=styr_fut;i<=endyr_fut;i++) 
    {
      double lb=value(future_biomass(k,i)/exp(2.*sqrt(log(1+square(future_biomass.sd(k,i))/square(future_biomass(k,i))))));
      double ub=value(future_biomass(k,i)*exp(2.*sqrt(log(1+square(future_biomass.sd(k,i))/square(future_biomass(k,i))))));
      R_report<<i<<" "<<future_biomass(k,i)<<" "<<future_biomass.sd(k,i)<<" "<<lb<<" "<<ub<<endl;
    }
  }
  double ctmp;
  for (int k=1;k<=5;k++){
    R_report<<"$Catch_fut_"<<k<<endl; 
    for (i=styr_fut;i<=endyr_fut;i++) 
    {
      if (k==5) ctmp=0.;else ctmp=value(catch_future(k,i));
      R_report<<i<<" "<<ctmp<<endl;
    }
  }

  R_report<<"$SSB"<<endl; for (i=styr_sp;i<=endyr+1;i++) 
  {
    double lb=value(Sp_Biom(i)/exp(2.*sqrt(log(1+square(Sp_Biom.sd(i))/square(Sp_Biom(i))))));
    double ub=value(Sp_Biom(i)*exp(2.*sqrt(log(1+square(Sp_Biom.sd(i))/square(Sp_Biom(i))))));
    R_report<<i<<" "<<Sp_Biom(i)<<" "<<Sp_Biom.sd(i)<<" "<<lb<<" "<<ub<<endl;
  }

  R_report<<"$R"<<endl; for (i=styr;i<=endyr;i++) 
  {
    double lb=value(recruits(i)/exp(2.*sqrt(log(1+square(recruits.sd(i))/square(recruits(i))))));
    double ub=value(recruits(i)*exp(2.*sqrt(log(1+square(recruits.sd(i))/square(recruits(i))))));
    R_report<<i<<" "<<recruits(i)<<" "<<recruits.sd(i)<<" "<<lb<<" "<<ub<<endl;
  }
    R_report << "$N"<<endl;
    for (i=styr;i<=endyr;i++) 
      R_report <<   i << " "<< natage(i) << endl;
      R_report   << endl;

    for (int k=1;k<=nfsh;k++)
    {
      R_report << "$F_age_"<< (k) <<""<< endl ;
      for (i=styr;i<=endyr;i++) 
        R_report <<i<<" "<<F(k,i)<<" "<< endl;
        R_report   << endl;
    }

    R_report <<endl<< "$Fshry_names"<< endl;
    for (int k=1;k<=nfsh;k++)
      R_report << fshname(k) << endl ;

    R_report <<endl<< "$Index_names"<< endl;
    for (int k=1;k<=nsrv;k++)
      R_report << srvname(k) << endl ;

    for (int k=1;k<=nsrv;k++)
    {
      int ii=1;
      R_report <<endl<< "$Obs_Survey_"<< k <<""<< endl ;
      for (i=styr;i<=endyr;i++)
      {
        if (ii<=yrs_srv(k).indexmax())
        {
          if (yrs_srv(k,ii)==i)
          {
            R_report << i<< " "<< obs_srv(k,ii) << " "<< pred_srv(k,i) <<" "<< obs_se_srv(k,ii) <<endl; //values of survey index value (annual)
            ii++;
          }
          else
            R_report << i<< " -1 "<< " "<< pred_srv(k,i)<<" -1 "<<endl;
        }
        else
          R_report << i<< " -1 "<< " "<< pred_srv(k,i)<<" -1 "<<endl;
      }
      R_report   << endl;
    }

    R_report << endl<< "$Survey_Q"<<endl;
    R_report<< q_srv << endl;
    R_report   << endl;
    for (int k=1;k<=nfsh;k++)
    {
      if (nyrs_fsh_age(k)>0) 
      { 
        R_report << "$pobs_fsh_"<< (k) <<""<< endl;
        for (i=1;i<=nyrs_fsh_age(k);i++) 
          R_report << yrs_fsh_age(k,i)<< " "<< oac_fsh(k,i) << endl;
        R_report   << endl;
      }
    }
    for (int k=1;k<=nfsh;k++)
    {
      if (nyrs_fsh_age(k)>0) 
      { 
        R_report << "$phat_fsh_"<< (k) <<""<< endl;
        for (i=1;i<=nyrs_fsh_age(k);i++) 
          R_report << yrs_fsh_age(k,i)<< " "<< eac_fsh(k,i) << endl;
          R_report   << endl;
      }
    }
    for (int k=1;k<=nsrv;k++)
    {
      if (nyrs_srv_age(k)>0) 
      { 
        R_report << "$pobs_srv_"<<(k)<<""<<  endl;
        for (i=1;i<=nyrs_srv_age(k);i++) 
          R_report << yrs_srv_age(k,i)<< " "<< oac_srv(k,i) << endl;
          R_report   << endl;
      }
    }
    for (int k=1;k<=nsrv;k++)
    {
      if (nyrs_srv_age(k)>0) 
      { 
        R_report << "$phat_srv_"<<(k)<<""<<  endl;
        for (i=1;i<=nyrs_srv_age(k);i++) 
          R_report << yrs_srv_age(k,i)<< " "<< eac_srv(k,i) << endl;
          R_report   << endl;
      }
    }
    for (int k=1;k<=nfsh;k++)
    {
      R_report << endl<< "$Obs_catch_"<<(k) << endl;
      R_report << catch_bio(k) << endl;
      R_report   << endl;
      R_report << "$Pred_catch_" <<(k) << endl;
      R_report << pred_catch(k) << endl;
      R_report   << endl;
    }

    for (int k=1;k<=nfsh;k++)
    {
      R_report << "$F_fsh_"<<(k)<<" "<<endl;
      for (i=styr;i<=endyr;i++)
      {
        R_report<< i<< " ";
        R_report<< mean(F(k,i)) <<" "<< mean(F(k,i))*max(sel_fsh(k,i)) << " ";
        R_report<< endl;
      }
    }

    for (int k=1;k<=nfsh;k++)
    {
      R_report << endl<< "$sel_fsh_"<<(k)<<"" << endl;
      for (i=styr;i<=endyr;i++)
        R_report << k <<"  "<< i<<" "<<sel_fsh(k,i) << endl; 
      R_report   << endl;
    }

    for (int k=1;k<=nsrv;k++)
    {
      R_report << endl<< "$sel_srv_"<<(k)<<"" << endl;
      for (i=styr;i<=endyr;i++)
        R_report << k <<"  "<< i<<" "<<sel_srv(k,i) << endl;
        R_report << endl;

    }
    R_report << endl<< "$Stock_Rec"<< endl;
    for (i=styr_rec;i<=endyr;i++)
      if (active(log_Rzero))
        R_report << i<< " "<<Sp_Biom(i-rec_age)<< " "<< SRecruit(Sp_Biom(i-rec_age))<< " "<< mod_rec(i)<<endl;
      else 
        R_report << i<< " "<<Sp_Biom(i-rec_age)<< " "<< " 999" << " "<< mod_rec(i)<<endl;
        
        R_report   << endl;

    R_report <<"$stock_Rec_Curve"<<endl;
    R_report <<"0 0"<<endl;
    dvariable stock;
    for (i=1;i<=30;i++)
    {
      stock = double (i) * Bzero /25.;
      if (active(log_Rzero))
        R_report << stock <<" "<< SRecruit(stock)<<endl;
      else
        R_report << stock <<" 99 "<<endl;
    }
    R_report   << endl;

    R_report   << endl<<"$Like_Comp" <<endl;
    obj_comps(11)= obj_fun - sum(obj_comps) ; // Residual
    obj_comps(12)= obj_fun ;

    R_report   <<obj_comps<<endl;
    R_report   << endl;
    R_report   << endl<<"$Like_Comp_names" <<endl;
    R_report   <<"catch_like     "<<endl
             <<"age_like_fsh     "<<endl
             <<"sel_like_fsh     "<<endl
             <<"surv_like        "<<endl
             <<"age_like_srv     "<<endl
             <<"sel_like_srv     "<<endl
             <<"rec_like         "<<endl
             <<"fpen             "<<endl
             <<"post_priors_srvq "<<endl
             <<"post_priors      "<<endl
             <<"residual         "<<endl
             <<"total            "<<endl;
    for (int k=1;k<=nfsh;k++)
    {
      R_report << "$Sel_Fshry_"<< (k) <<""<<endl;
      R_report << sel_like_fsh(k) << endl;
    }
    R_report   << endl;
  
    for (int k=1;k<=nsrv;k++)
    {
      R_report << "$Survey_Index_"<< (k) <<"" <<endl;
      R_report<< surv_like(k)<<endl;
    }
    R_report   << endl;

    R_report << setw(10)<< setfixed() << setprecision(5) <<endl;
    for (int k=1;k<=nsrv;k++)
    {
      R_report << "$Age_Survey_"<< (k) <<"" <<endl;
      R_report << age_like_srv(k)<<endl;
    }
    R_report   << endl;

    for (int k=1;k<=nsrv;k++)
    {
      R_report << "$Sel_Survey_"<< (k) <<""<<endl;
      R_report<< sel_like_srv(k,1) <<" "<<sel_like_srv(k,2)<<" "<<sel_like_srv(k,3)<< endl;
    }
    R_report   << endl;

    R_report << setw(10)<< setfixed() << setprecision(5) <<endl;
    R_report   << "$Rec_Pen" <<endl<<sigmar<<"  "<<rec_like<<endl;
    R_report   << endl;

    R_report   << "$F_Pen" <<endl;
    R_report<<fpen(1)<<"  "<<fpen(2)<<endl;
    R_report   << endl;
    for (int k=1;k<=nsrv;k++)
    {
      R_report << "$Q_Survey_"<< (k) <<""<<endl
             << " "<<post_priors_srvq(k)
             << " "<< q_srv(k)
             << " "<< qprior(k)
             << " "<< cvqprior(k)<<endl;
      R_report << "$Q_power_Survey_"<< (k) <<""<<endl
             << " "<<post_priors_srvq(k)
             << " "<< q_power_srv(k)
             << " "<< q_power_prior(k)
             << " "<< cvq_power_prior(k)<<endl;
    }
             R_report   << endl;
    R_report << "$M"<<endl;
    R_report << " "<< post_priors(1)
             << " "<< M
             << " "<< natmortprior
             << " "<< cvnatmortprior <<endl;
    R_report   << endl;
    R_report << "$Steep"<<endl;
    R_report << " "<< post_priors(2)
             << " "<< steepness
             << " "<< steepnessprior
             << " "<< cvsteepnessprior <<endl;
    R_report   << endl;
    R_report << "$Sigmar"<<endl;
    R_report << " "<< post_priors(3)
             << " "<< sigmar
             << " "<< sigmarprior
             << " "<< cvsigmarprior <<endl;
    R_report   << endl;
    R_report<<"$Num_parameters_Est"<<endl;
    R_report<<initial_params::nvarcalc()<<endl;
    R_report   << endl;
    
  R_report<<"$Steep_Prior" <<endl;
  R_report<<steepnessprior<<" "<<
    cvsteepnessprior<<" "<<
    phase_srec<<" "<< endl;
    R_report   << endl;

  R_report<<"$sigmarPrior " <<endl;
  R_report<<sigmarprior<<" "<<  cvsigmarprior <<" "<<phase_sigmar<<endl;
  R_report   << endl;

  R_report<<"$Rec_estimated_in_styr_endyr " <<endl;
  R_report<<styr_rec    <<" "<<endyr        <<" "<<endl;
  R_report   << endl;
  R_report<<"$SR_Curve_fit__in_styr_endyr " <<endl;
  R_report<<styr_rec_est<<" "<<endyr_rec_est<<" "<<endl;
  R_report   << endl;
  R_report<<"$Model_styr_endyr" <<endl;
  R_report<<styr        <<" "<<endyr        <<" "<<endl;
  R_report   << endl;

  R_report<<"$M_prior "<<endl;
  R_report<< natmortprior<< " "<< cvnatmortprior<<" "<<phase_M<<endl;
  R_report   << endl;
  R_report<<"$qprior " <<endl;
  R_report<< qprior<<" "<<cvqprior<<" "<< phase_q<<endl;
  R_report<<"$q_power_prior " <<endl;
  R_report<< q_power_prior<<" "<<cvq_power_prior<<" "<< phase_q_power<<endl;
  R_report   << endl;

  R_report<<"$cv_catchbiomass " <<endl;
  R_report<<cv_catchbiomass<<" "<<endl;
  R_report   << endl;
  R_report<<"$Projection_years"<<endl;
  R_report<< nproj_yrs<<endl;
  R_report   << endl;
  
  R_report << "$Fsh_sel_opt_fish "<<endl;
  for (int k=1;k<=nfsh;k++)
    R_report<<k<<" "<<fsh_sel_opt(k)<<" "<<sel_change_in_fsh(k)<<endl;
    R_report   << endl;
   R_report<<"$Survey_Sel_Opt_Survey " <<endl;
  for (int k=1;k<=nsrv;k++)
  R_report<<k<<" "<<(srv_sel_opt(k))<<endl;
  R_report   << endl;
    
  R_report <<"$Phase_survey_Sel_Coffs "<<endl;
  R_report <<phase_selcoff_srv<<endl;
  R_report   << endl;
  R_report <<"$Fshry_Selages " << endl;
  R_report << nselages_in_fsh  <<endl;
  R_report   << endl;
  R_report <<"$Survy_Selages " <<endl;
  R_report <<nselages_in_srv <<endl;
  R_report   << endl;

  R_report << "$Phase_for_age_spec_fishery"<<endl;
  R_report <<phase_selcoff_fsh<<endl;
  R_report   << endl;
  R_report << "$Phase_for_logistic_fishery"<<endl;
  R_report <<phase_logist_fsh<<endl;
  R_report   << endl;
  R_report << "$Phase_for_dble_logistic_fishery "<<endl;
  R_report <<phase_dlogist_fsh<<endl;
  R_report   << endl;

  R_report << "$Phase_for_age_spec_survey  "<<endl;
  R_report <<phase_selcoff_srv<<endl;
  R_report   << endl;
  R_report << "$Phase_for_logistic_survey  "<<endl;
  R_report <<phase_logist_srv<<endl;
  R_report   << endl;
  R_report << "$Phase_for_dble_logistic_srvy "<<endl;
  R_report <<phase_dlogist_srv<<endl;
  R_report   << endl;
  
  for (int k=1;k<=nfsh;k++)
  {
    if (nyrs_fsh_age(k)>0)
    {
      R_report <<"$EffN_Fsh_"<<(k)<<""<<endl;
      for (i=1;i<=nyrs_fsh_age(k);i++)
      {
        double sda_tmp = Sd_age(oac_fsh(k,i));
        R_report << yrs_fsh_age(k,i);
        R_report<< " "<<Eff_N(oac_fsh(k,i),eac_fsh(k,i)) ;
        R_report << " "<<Eff_N2(oac_fsh(k,i),eac_fsh(k,i));
        R_report << " "<<mn_age(oac_fsh(k,i));
        R_report << " "<<mn_age(eac_fsh(k,i));
        R_report << " "<<sda_tmp;
        R_report << " "<<mn_age(oac_fsh(k,i)) - sda_tmp *2. / sqrt(n_sample_fsh_age(k,i));
        R_report << " "<<mn_age(oac_fsh(k,i)) + sda_tmp *2. / sqrt(n_sample_fsh_age(k,i));
        R_report <<endl;
      }
    }
  }

  for (int k=1;k<=nfsh;k++)
  {
    R_report <<"$C_fsh_" <<(k)<<"" << endl; 
    for (i=styr;i<=endyr;i++)
      R_report <<i<<" "<<catage(k,i)<< endl;
  }

  R_report <<"$wt_a_pop" << endl<< wt_pop  <<endl;
  R_report <<"$mature_a" << endl<< maturity<<endl;
  for (int k=1;k<=nfsh;k++)
  {
    R_report <<"$wt_fsh_"<<(k)<<""<<endl;
    for (i=styr;i<=endyr;i++)
      R_report <<i<<" "<<wt_fsh(k,i)<< endl;
  }
  
  for (int k=1;k<=nsrv;k++)
  {
    R_report <<"$wt_srv_"<<(k)<<""<<endl;
    for (i=styr;i<=endyr;i++)
      R_report <<i<<" "<<wt_srv(k,i)<< endl;
  }
  for (int k=1;k<=nsrv;k++)
  {
    if (nyrs_srv_age(k)>0)
    {
      R_report <<"$EffN_Survey_"<<(k)<<""<<endl;
      for (i=1;i<=nyrs_srv_age(k);i++)
      {
        double sda_tmp = Sd_age(oac_srv(k,i));
        R_report << yrs_srv_age(k,i)
                << " "<<Eff_N(oac_srv(k,i),eac_srv(k,i)) 
                 << " "<<Eff_N2(oac_srv(k,i),eac_srv(k,i))
                 << " "<<mn_age(oac_srv(k,i))
                 << " "<<mn_age(eac_srv(k,i))
                 << " "<<sda_tmp
                 << " "<<mn_age(oac_srv(k,i)) - sda_tmp *2. / sqrt(n_sample_srv_age(k,i))
                 << " "<<mn_age(oac_srv(k,i)) + sda_tmp *2. / sqrt(n_sample_srv_age(k,i))
                 <<endl;
      }
    }
  }

  R_report.close();


FUNCTION double mn_age(_CONST dvector& pobs)
  // int lb1 = pobs.indexmin();
  // int ub1 = pobs.indexmax();
  // dvector av = age_vector(lb1,ub1)  ;
  // double mobs = value(pobs.shift(rec_age)*age_vector);
  double mobs = (pobs*age_vector);
  return mobs;

FUNCTION double mn_age(_CONST dvar_vector& pobs)
  // int lb1 = pobs.indexmin();
  // int ub1 = pobs.indexmax();
  // dvector av = age_vector(lb1,ub1)  ;
  // double mobs = value(pobs.shift(rec_age)*age_vector);
  double mobs = value(pobs*age_vector);
  return mobs;

FUNCTION double Sd_age(_CONST dvector& pobs)
  // double mobs = (pobs.shift(rec_age)*age_vector);
  // double stmp = (sqrt(elem_prod(age_vector,age_vector)*pobs.shift(rec_age) - mobs*mobs));
  double mobs = (pobs*age_vector);
  double stmp = sqrt((elem_prod(age_vector,age_vector)*pobs) - mobs*mobs);
  return stmp;

FUNCTION double Eff_N_adj(_CONST double, _CONST dvar_vector& pobs, _CONST dvar_vector& phat)
  int lb1 = pobs.indexmin();
  int ub1 = pobs.indexmax();
  dvector av = age_vector(lb1,ub1)  ;
  double mobs = value(pobs*av);
  double mhat = value(phat*av );
  double rtmp = mobs-mhat;
  double stmp = value(sqrt(elem_prod(av,av)*pobs - mobs*mobs));
  return square(stmp)/square(rtmp);

FUNCTION double Eff_N2(_CONST dvector& pobs, _CONST dvar_vector& phat)
  int lb1 = pobs.indexmin();
  int ub1 = pobs.indexmax();
  dvector av = age_vector(lb1,ub1)  ;
  double mobs =      (pobs*av);
  double mhat = value(phat*av );
  double rtmp = mobs-mhat;
  double stmp = (sqrt(elem_prod(av,av)*pobs - mobs*mobs));
  return square(stmp)/square(rtmp);

FUNCTION double Eff_N(_CONST dvector& pobs, _CONST dvar_vector& phat)
  dvar_vector rtmp = elem_div((pobs-phat),sqrt(elem_prod(phat,(1-phat))));
  double vtmp;
  vtmp = value(norm2(rtmp)/size_count(rtmp));
  return 1./vtmp;

FUNCTION double get_AC(_CONST int& srvind)
  // Functions to compute autocorrelation in residuals 
  int i1,i2,yy;
  i1 = 1;
  i2 = nyrs_srv(srvind);
  double actmp;
  dvector res(1,i2);
  for (i=1;i<=i2;i++)
  {
    yy = int(yrs_srv(srvind,i));
    cout<<yy<<" "<<obs_srv(srvind,i)<<" " <<pred_srv(srvind,yy)<<endl;
    res(i) = log(obs_srv(srvind,i)) - value(log(pred_srv(srvind,yy)));
  }
  double m1 = (mean(res(i1,i2-1)));
  double m2 = (mean(res(i1+1,i2))); 
  actmp = mean( elem_prod( ++res(i1,i2-1) - m1, res(i1+1,i2) - m2)) /
          (sqrt(mean( square(res(i1,i2-1) - m1 )))  * sqrt(mean(square(res(i1+1,i2) - m2 ))) );
  return(actmp);
