From fcbc604068005cd7cce6c3f329376c2b4ef7b540 Mon Sep 17 00:00:00 2001
From: timjmiller <timothy.j.miller@noaa.gov>
Date: Wed, 29 Sep 2021 11:56:48 -0400
Subject: [PATCH] q_model branch prior to testing. Need to merge in changes in
 devel

---
 R/asap3_to_wham.R         |   2 +-
 R/prepare_projection.R    |  23 +++-
 R/prepare_wham_input.R    |  40 ++++---
 R/set_ecov.R              | 206 ++++++++++++++++++++++-----------
 R/set_q.R                 |  49 ++++----
 R/set_random.R            |   2 +
 R/wham_plots_tables.R     |  58 +++++++++-
 man/prepare_wham_input.Rd |  51 ++++++--
 src/helper_functions.hpp  |  30 ++---
 src/wham_v0.cpp           | 236 ++++++++++++++++++++------------------
 10 files changed, 443 insertions(+), 254 deletions(-)

diff --git a/R/asap3_to_wham.R b/R/asap3_to_wham.R
index 778a50f1..cf4f2212 100644
--- a/R/asap3_to_wham.R
+++ b/R/asap3_to_wham.R
@@ -151,7 +151,7 @@ asap3_to_wham <- function(asap3, recruit_model = 2, model_name)
   data$bias_correct_pe = 0 #bias correct log-normal process errors?
   data$bias_correct_oe = 0 #bias correct log-normal observation errors?
   data$Fbar_ages = 1:data$n_ages
-  data$simulate_state = rep(1,4) #simulate state variables (NAA, M, sel, Ecov)
+  data$simulate_state = rep(1,5) #simulate state variables (NAA, M, sel, Ecov, q)
   data$simulate_data = rep(1,3) #simulate data types (catch, indices, Ecov)
   data$simulate_period = rep(1,2) #simulate above items for (model years, projection years)
   data$percentSPR = 40 #percentage of unfished SSB/R to use for SPR-based reference points
diff --git a/R/prepare_projection.R b/R/prepare_projection.R
index 555c94dd..7958036c 100644
--- a/R/prepare_projection.R
+++ b/R/prepare_projection.R
@@ -87,16 +87,17 @@ prepare_projection = function(model, proj.opts)
 
   # check ecov options are valid
   if(any(input1$data$Ecov_model > 0)){
+    n_effects = dim(input1$par$Ecov_beta)[1]
     # if Ecovs already extend beyond population model, Ecov proj options are unnecessary
-    n.beyond <- rep(NA, data$n_Ecov)
-    end.beyond <- rep(NA, data$n_Ecov)
-    for(i in 1:data$n_Ecov){
-      n.beyond[i] <- data$n_years_Ecov-1-data$ind_Ecov_out_end[i]
+    n.beyond <- end.beyond <- integer()
+    for(i in 1:data$n_Ecov) {
+      n.beyond[i] = data$n_years_Ecov-1-max(data$ind_Ecov_out_end[i,])
       end.beyond[i] <- min(n.beyond[i], proj.opts$n.yrs)     
       if(end.beyond[i] == proj.opts$n.yrs) print(paste0("ecov ",i," already fit through projection years. Using fit ecov ",i," for projections..."))
     }
     Ecov.proj <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
-    Ecov.map <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
+    #this isn't used anywhere
+    #Ecov.map <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
     if(all(end.beyond < proj.opts$n.yrs)){ # if Ecov proj options ARE necessary, check they are valid
       Ecov.opt.ct <- sum(proj.opts$cont.ecov, proj.opts$use.last.ecov, !is.null(proj.opts$avg.ecov.yrs), !is.null(proj.opts$proj.ecov))
       if(Ecov.opt.ct == 0) proj.opts$cont.ecov = TRUE; Ecov.opt.ct = 1;
@@ -210,6 +211,16 @@ prepare_projection = function(model, proj.opts)
     map$log_NAA = factor(tmp)
   }
 
+    # pad q_re
+    par$q_re <- rbind(par$q_re[1:data$n_years_model,,drop=F], matrix(0, data$n_years_proj, data$n_indices))
+    map$q_re <- rbind(matrix(as.integer(map$q_re), data$n_years_model, data$n_indices), matrix(NA, data$n_years_proj, data$n_indices))
+    if(any(data$use_q_re == 1)){
+      ind = which(data$use_q_re)
+      map$q_re[,ind] <- 1:(length(ind) * (data$n_years_model + data$n_years_proj)) #turn on appropriate columns of q_re
+    }
+    map$q_re <- factor(map$q_re)
+
+
   # pad Ecov_re for projections
   if(any(data$Ecov_model > 0)){
     if(any(end.beyond < proj.opts$n.yrs)){ # need to pad Ecov_re
@@ -217,7 +228,7 @@ prepare_projection = function(model, proj.opts)
         # for(i in 1:(proj.opts$n.yrs-end.beyond[j])){
         for(i in 1:max(proj.opts$n.yrs-end.beyond)){
           if(!is.null(proj.opts$use.last.ecov)) if(proj.opts$use.last.ecov){ # use last Ecov (pad Ecov_re but map to NA)
-            Ecov.proj[i,j] <- model$rep$Ecov_re[data$ind_Ecov_out_end[j]+1+end.beyond[j],j]
+            Ecov.proj[i,j] <- model$rep$Ecov_re[max(data$ind_Ecov_out_end[j,])+1+end.beyond[j],j]
           }
           if(!is.null(proj.opts$avg.ecov.yrs)){ # use average Ecov (pad Ecov_re but map to NA)
             ecov.yrs <- data$year1_Ecov+0:(data$n_years_Ecov-1+data$n_years_proj_Ecov)
diff --git a/R/prepare_wham_input.R b/R/prepare_wham_input.R
index cd735b71..2c54e54b 100644
--- a/R/prepare_wham_input.R
+++ b/R/prepare_wham_input.R
@@ -38,12 +38,18 @@
 #'     \item{$year}{Years corresponding to observations (vector of same length as \code{$mean} and \code{$logsigma})}
 #'     \item{$use_obs}{T/F (or 0/1) vector/matrix of the same dimension as \code{$mean} and \code{$logsigma}.
 #'     Use the observation? Can be used to ignore subsets of the ecov without changing data files.}
-#'     \item{$lag}{Offset between the ecov observations and their affect on the stock.
-#'     I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}.}
+#'     \item{$lag}{integer vector of offsets between the ecov observation/process and their affect on the stock.
+#'     I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}. May also be a list (length=n_Ecov) of vectors (lenght = 2+n_indices) if multiple effects of one or more Ecovs are modeled.}
 #'     \item{$process_model}{Process model for the ecov time-series. \code{"rw"} = random walk, \code{"ar1"} = 1st order autoregressive, \code{NA} = do not fit}
-#'     \item{$where}{Where does the ecov affect the population? \code{"recuit"} = recruitment,
-#'     \code{"M"} = natural mortality, \code{"none"} = fit ecov but without an effect on the population.}
-#'     \item{$how}{How does the ecov affect the \code{$where} process? These options are
+#'     \item{$where}{Character vector for where each ecov affects the population. \code{"recuit"} = recruitment,
+#'     \code{"M"} = natural mortality, \code{"q"} = catchability for indices, \code{"none"} = fit ecov but without an effect on the 
+#'     population. May also be a list (element for each ecov) of multiple character vectors so each ecov can multiple effects.}
+#'     \item{$indices}{indices that each ecov affects. Must be a list of length n.ecov, where each element is a vector of indices (1:n_indices). Must be provided when any of \code{where} = "q"}
+#'     \item{$link_model}{vector or (orthogonal polynomial order) models for effects or each ecov on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
+#'     where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b_0 + b_1*ecov + b_2*ecov^2 + ...}). A list (length = n_Ecov) of character vectors for modeling
+#'      effects on recruitment, M, and catchabilities for index 1,..., index n_indices (length of the vector is 2 + n_indices).}
+#'     \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
+#'     \item{$how}{How does the ecov affect the \code{$where} process? Currently this only isThese options are
 #'     specific to the \code{$where} process.
 #' Options for recruitment are described in \href{https://www.sciencedirect.com/science/article/pii/S1385110197000221}{Iles & Beverton (1998)}:
 #'   \describe{
@@ -58,10 +64,12 @@
 #'   \describe{
 #'     \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
 #'     \item{= 1}{effect on mean M (shared across ages)}
+#'   }
+#' Options for catchability:
+#'   \describe{
+#'     \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
+#'     \item{= 1}{effect on one or more catchabilities (see \code{$indices)})}
 #'   }}
-#'     \item{$link_model}{Model describing ecov effect on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
-#'     where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b0 + b1*ecov + b2*ecov^2 + ...}).}
-#'     \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
 #'   }
 #'
 #' \code{selectivity} specifies options for selectivity, to overwrite existing options specified in the ASAP data file.
@@ -166,11 +174,10 @@
 #'		 }
 #'   }
 #'
-#' \code{catchability} specifies options for catchability. If \code{NULL} and \code{asap3} is not NULL, a single catchability parameter used with initial values specified in ASAP file. If both are NULL, initial catchabilities for all indices = 0.3.
+#' \code{catchability} specifies options for catchability. If \code{NULL} and \code{asap3} is not NULL, a single catchability parameter for each index is used with initial values specified in ASAP file. If both are NULL, initial catchabilities for all indices = 0.3.
 #' Otherwise, it is a list with the following optional entries:
 #'   \describe{
 #'     \item{$re}{Time-varying (random effects) for each index. Vector with length = number of indices.
-#'                  If \code{NULL}, catchability parameters are constant over time.
 #'                  Each entry of \code{catchability$re} must be one of the following options:
 #'                  \describe{
 #'                    \item{"none"}{(default) are constant}
@@ -178,11 +185,14 @@
 #'                    \item{"ar1"}{correlated by year (AR1)}
 #'                  }
 #'                 }
-#'     \item{$initial_q}{Initial catchabilities for each index. vector length = number of indices. Will override values provided in \code{basic_info$q}. If \code{basic_info$q} and \code{asap3} are not provided default q values are 0.3.}
+#'     \item{$initial_q}{Initial catchabilities for each index. vector length = number of indices. Will override values provided in \code{basic_info$q}. 
+#'        If \code{basic_info$q} and \code{asap3} are not provided, default q values are 0.3.}
 #'     \item{$q_lower}{Lower bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 0.}
 #'     \item{$q_upper}{Upper bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 1000.}
 #'     \item{$fix_q}{Which catcabilities not estimated. vector of integers indicating which indices have catchability fixed. E.g. c(2,3) when the second and third catchabilties are fixed.}
-#'     \item{$prior}{vector of NULL and sigmas to use for gaussian prior on logit transform of catchability parameter. Length = number of indices. Indices with non-NULL values will have mean logit q as a random effect with mean determined by logit transform of \code{catchability$initial_q} and sigma as standard error.}
+#'     \item{$prior_sd}{vector of NA and standard deviations to use for gaussian prior on logit transform of catchability parameter. Length = number of indices. 
+#'       Indices with non-NA values will have mean logit q as a random effect with mean determined by logit transform of \code{catchability$initial_q} and  
+#'       sigma as standard error.}
 #'   }
 #'
 #' \code{age_comp} specifies the age composition models for fleet(s) and indices.
@@ -240,7 +250,7 @@
 #'     \item{$percentFMSY}{(0-100) percentage of Fmsy to use in projections.}
 #'     \item{$XSPR_R_avg_yrs}{which years to average recruitments for calculating SPR-based SSB reference points.}
 #'     \item{$XSPR_R_opt}{1(3): use annual R estimates(predictions) for annual SSB_XSPR, 2(4): use average R estimates(predictions).}
-#'     \item{$simulate_process_error}{T/F vector (length = 4). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, and Ecov. Only used if applicable.}
+#'     \item{$simulate_process_error}{T/F vector (length = 5). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, Ecov, and q. Only used if applicable.}
 #'     \item{$simulate_observation_error}{T/F vector (length = 3). When simulating from the model, whether to simulate  catch, index, and ecov observations.}
 #'     \item{$simulate_period}{T/F vector (length = 2). When simulating from the model, whether to simulate base period (model years) and projection period.}
 #'     \item{$bias_correct_process}{T/F. Perform bias correction of log-normal random effects for NAA.}
@@ -330,7 +340,7 @@ prepare_wham_input <- function(asap3 = NULL, model_name="WHAM for unnamed stock"
 	if(is.null(catchability)){
 		q_opts = NULL
 		q_names = c("q")
-		if(any(names(basic_info) %in% q_names)) q_opts$initial_q = basic_info[q_names]
+		if(any(names(basic_info) %in% q_names)) q_opts$initial_q = basic_info[[q_names]]
 	} else {
 		q_opts = catchability
 	}
@@ -476,7 +486,7 @@ initial_input_fn = function(input, basic_info){
 
   input$data$bias_correct_pe = 1 #bias correct log-normal process errors?
   input$data$bias_correct_oe = 1 #bias correct log-normal observation errors?
-  input$data$simulate_state = rep(1,4) #simulate state variables (NAA, M, sel, Ecov)
+  input$data$simulate_state = rep(1,5) #simulate state variables (NAA, M, sel, Ecov)
   input$data$simulate_data = rep(1,3) #simulate data types (catch, indices, Ecov)
   input$data$simulate_period = rep(1,2) #simulate above items for (model years, projection years)
   input$data$percentSPR = 40 #percentage of unfished SSB/R to use for SPR-based reference points
diff --git a/R/set_ecov.R b/R/set_ecov.R
index 0e83b56c..331b5eac 100644
--- a/R/set_ecov.R
+++ b/R/set_ecov.R
@@ -5,6 +5,11 @@ set_ecov = function(input, ecov)
   map = input$map
   #random = input$random
 
+  #define new dimensions for all effects a given ecov can have on assessment model
+  #currently the order is recruitment, M, index 1, ..., n_indices
+  n_effects = 2 + data$n_indices
+  index_effects = 2+1:data$n_indices
+
   # --------------------------------------------------------------------------------
   # Environmental covariate data
   if(is.null(ecov)){
@@ -19,15 +24,17 @@ set_ecov = function(input, ecov)
     data$Ecov_year <- matrix(0, nrow=1, ncol=1)
     data$year1_Ecov <- 0
     data$year1_model <- input$years[1]
-    data$Ecov_lag <- 0
+    #data$Ecov_lag <- 0 #This is not used anywhere
     data$Ecov_model <- 0
-    data$Ecov_where <- 0
+    data$Ecov_where = matrix(0, data$n_Ecov, n_effects)
     data$Ecov_how <- 0
-    data$Ecov_poly <- 1
+    #not used
+    #data$Ecov_poly <- 1
     data$n_years_Ecov <- 1
-    data$ind_Ecov_out_start <- data$ind_Ecov_out_end <- 0
+    data$ind_Ecov_out_start <- data$ind_Ecov_out_end <- matrix(0, data$n_Ecov, n_effects)
     data$Ecov_label <- "none"
     data$Ecov_use_re <- matrix(0, nrow=1, ncol=1)
+    Ecov_poly <- matrix(1,data$n_Ecov, n_effects)
   } else {
     if(class(ecov$mean)[1] == "matrix") {data$Ecov_obs <- ecov$mean} else{
       warning("ecov$mean is not a matrix. Coercing to a matrix...")
@@ -101,14 +108,18 @@ set_ecov = function(input, ecov)
 
     # # check that Ecov year vector doesn't have missing gaps
     # pad Ecov if it starts after model year1 - max(lag)
-    if(data$year1_Ecov > data$year1_model - max(ecov$lag)){
-      print("ecov does not start by model year 1 - max(lag). Padding ecov...")
-      data$Ecov_obs <- rbind(matrix(0, nrow = data$year1_Ecov-(data$year1_model-max(ecov$lag)), ncol = data$n_Ecov), data$Ecov_obs)
-      par.Ecov.obs.logsigma <- rbind(matrix(-1.3, nrow = data$year1_Ecov-(data$year1_model-max(ecov$lag)), ncol = data$n_Ecov), par.Ecov.obs.logsigma)
-      map.Ecov.obs.logsigma <- rbind(matrix(NA, nrow = data$year1_Ecov-(data$year1_model-max(ecov$lag)), ncol = data$n_Ecov), map.Ecov.obs.logsigma)
-      data$Ecov_use_obs <- rbind(matrix(0, nrow = data$year1_Ecov-(data$year1_model-max(ecov$lag)), ncol = data$n_Ecov), data$Ecov_use_obs)
-      data$Ecov_year <- c(seq(data$year1_model - max(ecov$lag), data$year1_Ecov-1), data$Ecov_year)
-      data$year1_Ecov <- data$year1_model - max(ecov$lag)
+    #make ecov$lag compatible with multiple effects
+    if(is.null(ecov$lag)) stop("ecov$lag needs to be provided for each ecov")
+    if(!is.list(ecov$lag)) ecov$lag = lapply(ecov$lag, function(x) rep(x,n_effects))
+    max.lag = max(unlist(ecov$lag))
+    if(data$year1_Ecov > data$year1_model - max.lag){
+      print("one or more ecov does not start by model year 1 - max(lag). Padding ecov...")
+      data$Ecov_obs <- rbind(matrix(0, nrow = data$year1_Ecov-(data$year1_model-max.lag), ncol = data$n_Ecov), data$Ecov_obs)
+      par.Ecov.obs.logsigma <- rbind(matrix(-1.3, nrow = data$year1_Ecov-(data$year1_model-max.lag), ncol = data$n_Ecov), par.Ecov.obs.logsigma)
+      map.Ecov.obs.logsigma <- rbind(matrix(NA, nrow = data$year1_Ecov-(data$year1_model-max.lag), ncol = data$n_Ecov), map.Ecov.obs.logsigma)
+      data$Ecov_use_obs <- rbind(matrix(0, nrow = data$year1_Ecov-(data$year1_model-max.lag), ncol = data$n_Ecov), data$Ecov_use_obs)
+      data$Ecov_year <- c(seq(data$year1_model - max.lag, data$year1_Ecov-1), data$Ecov_year)
+      data$year1_Ecov <- data$year1_model - max.lag
     }
 
     # pad Ecov if it ends before last model year
@@ -125,73 +136,100 @@ set_ecov = function(input, ecov)
     data$Ecov_use_re <- matrix(1, nrow=data$n_years_Ecov, ncol=data$n_Ecov)
 
     # get index of Ecov_x to use for Ecov_out (Ecovs can have diff lag)
-    data$ind_Ecov_out_start <- data$ind_Ecov_out_end <- rep(NA, data$n_Ecov)
-    for(i in 1:data$n_Ecov){
-      data$ind_Ecov_out_start[i] <- which(data$Ecov_year==data$year1_model)-ecov$lag[i]-1 # -1 is for cpp indexing
-      data$ind_Ecov_out_end[i] <- which(data$Ecov_year==end_model)-ecov$lag[i]-1 # -1 is for cpp indexing
+    data$ind_Ecov_out_start <- data$ind_Ecov_out_end <- matrix(NA, data$n_Ecov, n_effects)
+    for(i in 1:data$n_Ecov) for(j in 1:n_effects) {
+      data$ind_Ecov_out_start[i,j] <- which(data$Ecov_year==data$year1_model)-ecov$lag[[i]][j]-1 # -1 is for cpp indexing
+      data$ind_Ecov_out_end[i,j] <- which(data$Ecov_year==end_model)-ecov$lag[[i]][j]-1 # -1 is for cpp indexing
     }
 
-    if(!identical(length(ecov$lag), length(ecov$label), data$n_Ecov)) stop("Length of Ecov_lag and Ecov_label vectors not equal to # Ecov")
-    data$Ecov_lag <- ecov$lag
-    if(!all(ecov$process_model %in% c(NA,"rw", "ar1"))){
-      stop("ecov$process_model must be 'rw' (random walk), 'ar1', or NA (do not fit)")
-    }
-    if(is.na(ecov$process_model) && ecov$how !=0){
-      stop("ecov$process_model not chosen (NA) but ecov$how specified.
-       Either 1) choose an ecov process model ('rw' or 'ar1'),
-              2) turn off ecov (set ecov$how = 0 and ecov$process_model = NA),
-           or 3) fit ecov but with no effect on population (ecov$how = 0, ecov$process_model 'rw' or 'ar1').")
-    }
-    data$Ecov_model <- sapply(ecov$process_model, match, c("rw", "ar1"))
+    if(!identical(length(ecov$lag), length(ecov$label), data$n_Ecov)) stop("Length of ecov$lag and ecov$label not equal to # Ecov")
+    ecov$lag <- t(matrix(unlist(ecov$lag), n_effects, data$n_Ecov)) #just used on R side (e.g., plotting)
+    
+    data$Ecov_where = matrix(0, data$n_Ecov, n_effects)
+    if(is.character(ecov$where)) ecov$where = as.list(ecov$where) #put in new allowed format so each ecov can have mulitple effects.
 
-    if(any(ecov$where == 'recruit') & data$n_NAA_sigma == 0){
+    if(any(sapply(ecov$where, function(x) any(x == 'recruit'))) & data$n_NAA_sigma == 0){
       stop("Cannot estimate ecov effect on recruitment when
       recruitment in each year is estimated freely as fixed effect parameters.
       Either remove ecov-recruit effect or estimate recruitment
       (or all numbers-at-age) as random effects.")
     }
-    if(is.null(ecov$where)) stop("ecov$where must be specified, 'recruit' or 'M'")
-    if(!any(ecov$where %in% c('recruit','M','none'))){
-      stop("Only ecov effects on recruitment and M currently implemented.
-      Set ecov$where = 'recruit', 'M', or 'none'.")
+    if(is.null(ecov$where)) stop("ecov$where must be specified, 'recruit', 'M', or 'q'")
+    if(any(sapply(ecov$where, function(x) any(!(x %in% c('recruit','M','q','none')))))){
+      stop("Only ecov effects on recruitment, M, and catchability (q) currently implemented.
+      Set ecov$where = 'recruit', 'M', 'q', or 'none'.")
     }
-    data$Ecov_where <- sapply(ecov$where, match, c('none','recruit','M')) - 1
 
-    if(is.null(ecov$how)) stop("ecov$how must be specified")
-    if(length(ecov$how) != data$n_Ecov) stop("ecov$how must be a vector of length(n.ecov)")
+    if(!all(ecov$process_model %in% c(NA,"rw", "ar1"))){
+      stop("ecov$process_model must be 'rw' (random walk), 'ar1', or NA (do not fit)")
+    }
+    for(i in 1:data$n_Ecov) {
+      if(is.na(ecov$process_model[i]) & any(ecov$where[[i]] %in% c("recruit", "M", "q"))){ #ecov$how !=0){
+      stop(paste0("ecov$process_model ", i, " is turned off (NA) but ecov$where includes 'M', 'recruit', or 'q'.
+       Either 1) choose an ecov process model ('rw' or 'ar1'),
+              2) turn off ecov (set ecov$where[i] = 'none' and ecov$process_model = NA),
+           or 3) fit ecov but with no effect on population (ecov$where[i] = 'none', ecov$process_model[i] = 'rw' or 'ar1')."))
+      }
+    }
+    data$Ecov_model <- sapply(ecov$process_model, match, c("rw", "ar1"))
     for(i in 1:data$n_Ecov){
-      if(data$Ecov_where[i] == 2) if(!ecov$how[i] %in% c(0,1)){
-        stop("Sorry, only the following ecov effects on M are currently implemented.
-        Set ecov$how = 0 (no effect) or 1 (effect on mean M, shared across ages).")
+      if(any(ecov$where[[i]] == "recruit")) data$Ecov_where[i,1] = 1 
+      if(any(ecov$where[[i]] == "M")) data$Ecov_where[i,2] = 1 
+      if(any(ecov$where[[i]] == "q")) {
+        if(is.null(ecov$indices)) stop("ecov$indices must be specified if any ecov$where == 'q'") 
+        if(!is.list(ecov$indices)) stop("ecov$indices must be a specified as a list (length = n_Ecov) of vectors of any indices each Ecov affects") 
+        data$Ecov_where[i,2 + ecov$indices[[i]]] = 1 
       }
-      if(data$Ecov_where[i] == 1) if(!ecov$how[i] %in% c(0,1,2,4)){
+    }
+    #data$Ecov_where <- sapply(ecov$where, match, c('none','recruit','M','q')) - 1
+
+    if(is.null(ecov$how) & ('recruit' %in% unlist(ecov$where))) stop("ecov$how must be specified when any ecov is affecting recruitment")
+    if(length(ecov$how) != data$n_Ecov) stop("ecov$how must be a vector of length(n.ecov)")
+    data$Ecov_how = rep(0, data$n_Ecov) 
+    for(i in 1:data$n_Ecov){
+      #if(data$Ecov_where[i,2] == 1) if(!ecov$how[i] %in% c(0,1)){
+      #  stop("Sorry, only the following ecov effects on M are currently implemented.
+      #  Set ecov$how = 0 (no effect) or 1 (effect on mean M, shared across ages).")
+      #}
+      #if(any(data$Ecov_where[i,index_effects] == 1)) if(!ecov$how[i] %in% c(0,1)){
+      #  stop("Sorry, only the following ecov effects on q are currently implemented.
+      #  Set ecov$how = 0 (no effect) or 1 (effect on q).")
+      #}
+      if(data$Ecov_where[i,1] == 1) if(!ecov$how[i] %in% c(0,1,2,4)){
         stop("Sorry, only the following ecov effects on recruitment are currently implemented.
         Set ecov$how = 0 (no effect), 1 (controlling), 2 (limiting, Bev-Holt only), or 4 (masking).")
       }
-      if(data$Ecov_where[i] == 1 & data$recruit_model == 1){
+      if(data$Ecov_where[i,1] == 1 & data$recruit_model == 1){
         stop("Random walk recruitment cannot have an ecov effect on recruitment.
         Either choose a different recruit_model (2, 3, or 4), or remove the Ecov effect.")
       }
-      if(data$Ecov_where[i] == 1) if(data$recruit_model == 4 & ecov$how[i] == 2){
+      if(data$Ecov_where[i,1] == 1) if(data$recruit_model == 4 & ecov$how[i] == 2){
         stop("'Limiting' ecov effect on Ricker recruitment not implemented.
         Either set ecov$how = 0 (no effect), 1 (controlling), or 4 (masking)...
         Or set recruit_model = 3 (Bev-Holt).")
       }      
+      #currently only need this if recruitment effects modeled.
+      if(data$Ecov_where[i,1] == 1) data$Ecov_how[i] = ecov$how[i]
     }
-    data$Ecov_how <- ecov$how
-    data$Ecov_poly <- rep(1,data$n_Ecov)
-    ecov_str <- as.list(rep('linear',data$n_Ecov))
+
+    #data$Ecov_how <- ecov$how
+    #data$Ecov_poly is not used in wham_v0.cpp
+    Ecov_poly <- matrix(1,data$n_Ecov, n_effects)
+    #ecov_str <- list()
+    #ecov_str <- as.list(rep('linear',data$n_Ecov))
     if(!is.null(ecov$link_model)){
       if(!is.na(ecov$link_model)){
-        for(i in 1:data$n_Ecov){
-          ecov_str[[i]] <- strsplit(ecov$link_model[i],"-")[[1]]
-          if(!ecov_str[[i]][1] %in% c('linear','poly')) stop("Only 'linear' or 'poly-x' (x = 1, 2, ...) ecov link models implemented.")
-          if(ecov_str[[i]][1]=='linear') data$Ecov_poly[i] <- 1
-          if(ecov_str[[i]][1]=='poly') data$Ecov_poly[i] <- as.numeric(ecov_str[[i]][2])
-          ecov_str[[i]] = ecov$link_model[i]
+        if(!is.list(ecov$link_model)) ecov$link_model = lapply(ecov$link_model, rep, n_effects)
+        for(i in 1:data$n_Ecov) for(j in 1:n_effects){
+          ecov_str <- strsplit(ecov$link_model[[i]][j],"-")[[1]]
+          if(!ecov_str[1] %in% c('linear','poly')) stop("Only 'linear' or 'poly-x' (x = 1, 2, ...) ecov link models implemented.")
+          #if(ecov_str[1]=='linear') Ecov_poly[i,j] <- 1 #already supllied above
+          if(ecov_str[1]=='poly') Ecov_poly[i,j] <- as.numeric(ecov_str[2])
+          #ecov_str[[i]] = ecov$link_model[i]
         }
       }
     }
+
     if(!is.null(ecov$ages)){
       if(length(ecov$ages) != data$n_Ecov) stop("ecov$ages must be a list of length n.ecov")
       for(i in 1:data$n_Ecov){
@@ -202,46 +240,76 @@ set_ecov = function(input, ecov)
       for(i in 1:data$n_Ecov) ecov$ages[[i]] <- 1:data$n_ages # default: ecov affects all ages
     }
 
-    cat(paste0("Please check that the environmental covariates have been loaded
-and interpreted correctly.
+    cat(paste0("Please check that the environmental covariates have been loaded and interpreted correctly.
 
 Model years: ", data$year1_model, " to ", end_model,"
 Ecov years: ", data$year1_Ecov, " to ", end_Ecov,"
 
 "))
+
     for(i in 1:data$n_Ecov){
       years <- data$Ecov_year[as.logical(data$Ecov_use_obs[,i])]
+      lastyr <- tail(years,1)
 
-      if(data$Ecov_where[i] == 1){ # recruitment
+      if(data$Ecov_where[i,1] == 1){ # recruitment
         cat(paste0("Ecov ",i,": ",ecov$label[i],"
-",c('*NO*','Controlling','Limiting','Lethal','Masking','Directive')[data$Ecov_how[i]+1]," (",ecov_str[[i]],") effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
+",c('*NO*','Controlling','Limiting','Lethal','Masking','Directive')[data$Ecov_how[i]+1]," (",ecov$link[[i]][1],") effect on: recruitment
 
-In model years:
+Model years:
 "))
+
+        cat(years, fill=TRUE)
+
+        cat(paste0("Lag: ",ecov$lag[i,1],"
+        Ex: ",ecov$label[i]," in ",years[1]," affects recruitment in ",years[1+ecov$lag[i,1]],"
+            ",ecov$label[i]," in ",lastyr," affects recruitment in ",lastyr+ecov$lag[i,1],"
+
+        "))
       }
-      if(data$Ecov_where[i] == 2){ # M
-        cat(paste0("Ecov ",i,": ",ecov$label[i],"
-",c('*NO*',ecov_str[[i]])[data$Ecov_how[i]+1]," effect on: ", c('recruitment','M')[data$Ecov_where[i]],"
 
-In model years:
+      if(data$Ecov_where[i,2] == 1){ # M
+        cat(paste0("Ecov ",i,": ",ecov$label[i]," effect (", ecov$link[[i]][2], ") on: M 
+
+Model years:
 "))
+        cat(years, fill=TRUE)
+
+        cat(paste0("Lag: ",ecov$lag[i,2],"
+        Ex: ",ecov$label[i]," in ",years[1]," affects M in ",years[1+ecov$lag[i,2]],"
+            ",ecov$label[i]," in ",lastyr," affects M in ",lastyr+ecov$lag[i,2],"
+
+        "))
       }
 
-cat(years, fill=TRUE)
-lastyr <- tail(years,1)
-cat(paste0("Lag: ",data$Ecov_lag[i],"
-Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",years[1+data$Ecov_lag[i]],"
-    ",ecov$label[i]," in ",lastyr," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",lastyr+data$Ecov_lag[i],"
+      for(j in index_effects) if(data$Ecov_where[i,j] == 1){ # q
+        cat(paste0("Ecov ",i,": ",ecov$label[i]," effect (", ecov$link[[i]][j], ") on: q for index ", j + 1 - min(index_effects), " 
 
+Model years:
 "))
+        cat(years, fill=TRUE)
+
+        cat(paste0("Lag: ",ecov$lag[i,j],"
+        Ex: ",ecov$label[i]," in ",years[1]," affects index ", j + 1 - min(index_effects), " in ",years[1+ecov$lag[i,j]],"
+            ",ecov$label[i]," in ",lastyr," affects M index ", j + 1 - min(index_effects), " in ",lastyr+ecov$lag[i,j],"
+
+        "))
+      }
+
+#cat(years, fill=TRUE)
+
+#cat(paste0("Lag: ",ecov$lag[i],"
+#Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",years[1+ecov$lag[i]],"
+#    ",ecov$label[i]," in ",lastyr," affects ", c('recruitment','M')[data$Ecov_where[i]]," in ",lastyr+ecov$lag[i],"
+#
+#"))
     }
     data$Ecov_label <- list(data$Ecov_label)
   } # end load Ecov
 
   # Ecov pars
   par$Ecov_re = matrix(rnorm(data$n_years_Ecov*data$n_Ecov), data$n_years_Ecov, data$n_Ecov)
-  max.poly <- max(data$Ecov_poly)
-  par$Ecov_beta = array(0, dim=c(max.poly, data$n_Ecov, data$n_ages)) # beta_R in eqns 4-5, Miller et al. (2016)
+  max.poly <- max(Ecov_poly)
+  par$Ecov_beta = array(0, dim=c(n_effects, max.poly, data$n_Ecov, data$n_ages)) # beta_R in eqns 4-5, Miller et al. (2016)
   par$Ecov_process_pars = matrix(0, 3, data$n_Ecov) # nrows = RW: 2 par (Ecov1, log_sig), AR1: 3 par (mu, log_sig, phi); ncol = N_ecov
   par$Ecov_process_pars[1,] = -1.3 # start sig_ecov at 0.27
   par$Ecov_obs_logsigma <- par.Ecov.obs.logsigma
@@ -256,10 +324,10 @@ Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_
   # turn off Ecov_beta to fit Ecov process model without effect on population
   tmp <- array(NA, dim=dim(par$Ecov_beta))
   ct = 1
-  for(j in 1:data$n_Ecov){
+  for(n in 1:n_effects) for(j in 1:data$n_Ecov){
     for(i in 1:max.poly){
       for(a in 1:data$n_ages){
-        if(data$Ecov_how[j] > 0 & i <= data$Ecov_poly[j] & a %in% ecov$ages[[j]]) tmp[i,j,a] = ct # default share ecov_beta across ages
+        if(data$Ecov_where[j,n] == 1 & i <= Ecov_poly[j,n] & a %in% ecov$ages[[j]]) tmp[n,i,j,a] = ct # default share ecov_beta across ages
       }
       ct = ct+1
     }
diff --git a/R/set_q.R b/R/set_q.R
index dbd96f7d..7a9a4b3e 100644
--- a/R/set_q.R
+++ b/R/set_q.R
@@ -10,45 +10,50 @@ set_q = function(input, q_opts = NULL){
   data$q_lower <- rep(0,data$n_indices)
   data$q_upper <- rep(1000,data$n_indices)
   if(!is.null(q_opts$q_lower)) {
-  	ind = which(!is.null(q_opts$q_lower))
-  	data$q_lower[ind] = q_opts$q_lower[ind]
+		if(length(q_opts$q_lower) != data$n_indices) stop("the length of q_opts$q_lower is not equal to the number of indices")
+  	data$q_lower = q_opts$q_lower
   }
   if(!is.null(q_opts$q_upper)) {
-  	ind = which(!is.null(q_opts$q_upper))
-  	data$q_upper[ind] = q_opts$q_upper[ind]
+		if(length(q_opts$q_upper) != data$n_indices) stop("the length of q_opts$q_upper is not equal to the number of indices")
+  	data$q_upper = q_opts$q_upper
   }
-	par$logit_q = rep(gen.logit(0.3, data$q_lower, data$q_upper), data$n_indices)
+	par$logit_q = gen.logit(0.3, data$q_lower, data$q_upper)
   if(!is.null(asap3)) par$logit_q = gen.logit(asap3$q_ini[which(asap3$use_index ==1)], data$q_lower, data$q_upper) # use q_ini values from asap3 file
 	if(!is.null(q_opts$initial_q)) {
-  	ind = which(!is.null(q_opts$initial_q))
-		par$logit_q[ind] = gen.logit(q_opts$initial_q[ind], data$q_lower[ind], data$q_upper[ind])
+		if(length(q_opts$initial_q) != data$n_indices) stop("the length of q_opts$initial_q is not equal to the number of indices")
+		par$logit_q = gen.logit(q_opts$initial_q, data$q_lower, data$q_upper)
 	}
-	data$q_re_model = rep(0, data$n_indices)
-	par$q_re = matrix(0, data$n_years_model + data$n_years_proj, data$n_indices)
-	map$q_re = matrix(NA, data$n_years_model + data$n_years_proj, data$n_indices)
-	par$q_re_pars = matrix(0, data$n_indices, 2)
-	map$q_re_pars = matrix(NA, data$n_indices, 2)
-	if(!is.null(q_opts$model)){
-  	ind = which(q_opts$model %in% c("iid", "ar1")
-		map$q_re[,ind] = 1:(length(ind) * (data$n_years_model + data$n_years_proj)) #turn on appropriate columns of q_re
-		data$q_re_model[ind] = 1
-		map$q_re_pars[ind,1] = 1:length(ind)
-  	ind = which(q_opts$model == "ar1")
-		map$q_re_pars[ind,2] = sum(!is.na(map$q_re_pars[,1])) + 1:length(ind)
+	data$use_q_re = rep(0, data$n_indices)
+	par$q_re = matrix(0, data$n_years_model, data$n_indices)
+	map$q_re = matrix(NA, data$n_years_model, data$n_indices)
+	par$q_repars = matrix(0, data$n_indices, 2)
+	map$q_repars = matrix(NA, data$n_indices, 2)
+	if(!is.null(q_opts$re)){
+  	ind = which(q_opts$re %in% c("iid", "ar1"))
+		map$q_re[,ind] = 1:(length(ind) * (data$n_years_model)) #turn on appropriate columns of q_re
+		data$use_q_re[ind] = 1
+		map$q_repars[ind,1] = 1:length(ind)
+  	ind = which(q_opts$re == "ar1")
+		map$q_repars[ind,2] = sum(!is.na(map$q_repars[,1])) + 1:length(ind)
 	}
-	map$q_re_pars = factor(map$q_re_pars)
+	map$q_repars = factor(map$q_repars)
+  map$q_re = factor(map$q_re)
 
 	data$use_q_prior = rep(0, data$n_indices)
 	par$q_prior_re = rep(0, data$n_indices)
 	map$q_prior_re = rep(NA, data$n_indices)
 	data$logit_q_prior_sigma = rep(1, data$n_indices)
+	map$logit_q = 1:data$n_indices
 	if(!is.null(q_opts$prior)){
-  	ind = which(!is.null(q_opts$prior))
+  	ind = which(!is.na(q_opts$prior_sd))
 		data$use_q_prior[ind] = 1
-		data$logit_q_prior_sigma[ind] = q_opts$prior[ind]
+		data$logit_q_prior_sigma[ind] = q_opts$prior_sd[ind]
+		par$q_prior_re[ind] = par$logit_q[ind]
 		map$q_prior_re[ind] = 1:length(ind)
+		map$logit_q[ind] = NA #turn off logit q (mean of the prior) when random effects are used
 	}
 	map$q_prior_re = factor(map$q_prior_re)
+	map$logit_q = factor(map$logit_q)
 
   input$data = data
   input$par = par
diff --git a/R/set_random.R b/R/set_random.R
index 092b0ae5..e105ff74 100644
--- a/R/set_random.R
+++ b/R/set_random.R
@@ -5,6 +5,8 @@ set_random <- function(input){
   if(input$data$M_re_model > 1) random = c(random, "M_re")
   if(sum(input$data$Ecov_model) > 0) random = c(random, "Ecov_re")
   if(input$data$n_NAA_sigma > 0) random = c(random, "log_NAA")
+  if(sum(input$data$use_q_prior)) random = c(random, "q_prior_re")
+  if(sum(input$data$use_q_re)) random = c(random, "q_re")
   input$random = random	
   return(input)
 }
\ No newline at end of file
diff --git a/R/wham_plots_tables.R b/R/wham_plots_tables.R
index 55a9a3f0..2079f698 100644
--- a/R/wham_plots_tables.R
+++ b/R/wham_plots_tables.R
@@ -258,7 +258,7 @@ fit.summary.text.plot.fn <- function(mod){
   recs <- c("Random walk","Random about mean","Bev-Holt","Ricker")
   env.mod <- c("RW", "AR1")
   # env.where <- c('Recruitment','Growth','Mortality')
-  env.where <- c('Recruitment','Mortality')  
+  env.where <- c('Recruitment','Mortality',paste0("q for index ", 1:mod$env$data$n_indices))  
   env.how <- c("controlling", "limiting", "lethal", "masking", "directive")
   fleet_selblocks = lapply(1:mod$env$data$n_fleets, function(x) unique(mod$env$data$selblock_pointer_fleets[,x]))
   index_selblocks = lapply(1:mod$env$data$n_indices, function(x) unique(mod$env$data$selblock_pointer_indices[,x]))
@@ -275,9 +275,19 @@ fit.summary.text.plot.fn <- function(mod){
   text(5,nl <- nl-0.5, paste0("Index Age Comp Models: ", paste(acm[mod$env$data$age_comp_model_indices], collapse = ", ")))
   text(5,nl <- nl-0.5,paste0("Recruitment model: ", recs[mod$env$data$recruit_model]))
   if(!all(mod$env$data$Ecov_model == 0)){
-    for(ec in 1:mod$env$data$n_Ecov) text(5,nl <- nl-0.5, paste0("Environmental effect ", ec,": ", mod$env$data$Ecov_label[[1]][ec]," (",env.mod[mod$env$data$Ecov_model[ec]],") on ",env.where[mod$env$data$Ecov_where[ec]], " (", env.how[mod$env$data$Ecov_how[ec]],")"))
+    for(ec in 1:mod$env$data$n_Ecov) {
+      ec.where = env.where[which(mod$env$data$Ecov_where[ec,]==1)]
+      ec.mod = env.mod[mod$env$data$Ecov_model[ec]]
+      ec.label = mod$env$data$Ecov_label[[1]][ec]
+      if(length(ec.mod)) {
+        out = paste0("Environmental covariate ", ec,": ", ec.label," modeled as ",ec.mod,".")
+        if(length(ec.where)) out = paste0(out, paste0(" Effects on ", paste(ec.where,collapse = ','), " estimated."))
+        else out = paste(out, " No effects estimated.")
+        text(5,nl <- nl-0.5, out)
+      }
+    }
   } else {
-    text(5,nl <- nl-0.5, "Environmental effects: none")
+    text(5,nl <- nl-0.5, "No Environmental covariates.")
   }
   text(5,nl <- nl-0.5,paste0("Number of Selectivity blocks: ", mod$env$data$n_selblocks))
   text(5,nl <- nl-0.5, paste0("Selectivity Block Types: ", paste(selmods[mod$env$data$selblock_models], collapse = ", ")))
@@ -1356,6 +1366,48 @@ plot.index.sel.blocks <- function(mod, ages, ages.lab, plot.colors, do.tex = FAL
   # par(origpar)
 }
 
+plot.q.trend<-function(mod, alpha = 0.05)
+{
+  origpar <- par(no.readonly = TRUE)
+  years_full <- mod$years_full
+  years <- mod$years
+  tcol <- col2rgb('black')
+  tcol <- paste(rgb(tcol[1,],tcol[2,], tcol[3,], maxColorValue = 255), "55", sep = '')
+  if(class(mod$sdrep)[1] == "sdreport"){
+    std = summary(mod$sdrep)
+  } else {
+    std = mod$sdrep
+  }
+  par(mfrow=c(2,1), mar=c(1,1,1,1), oma = c(4,4,0,0))
+
+  ssb.ind <- which(rownames(std) == "log_SSB")
+  log.ssb <- std[ssb.ind,1]
+  ssb = exp(log.ssb)/1000
+  ssb.cv <- std[ssb.ind,2]
+  log.ssb.ci <- log.ssb + cbind(qnorm(1-alpha/2)*ssb.cv, -qnorm(1-alpha/2)*ssb.cv)
+  ssb.ci = exp(log.ssb.ci)/1000
+  no.ssb.ci <- all(is.na(ssb.ci))
+  if(!no.ssb.ci){ # have CI
+    plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,max(ssb.ci)), axes = FALSE)
+    axis(1, labels = FALSE)
+    axis(2)
+    box()
+    mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
+    grid(col = gray(0.7))
+    polygon(c(years_full,rev(years_full)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
+  } else { # no CI but plot SSB trend
+    plot(years_full, ssb, type='l', lwd=2, xlab="", ylab="", ylim=c(0,max(ssb)), axes = FALSE)
+    axis(1, labels = FALSE)
+    axis(2)
+    box()
+    mtext(side = 2, "SSB (kmt)", outer = FALSE, line = 3)
+    grid(col = gray(0.7))
+    # polygon(c(years,rev(years)), c(ssb.ci[,1],rev(ssb.ci[,2])), col = tcol, border = tcol, lwd = 1)
+  }
+  if(length(years_full) > length(years)) abline(v=tail(years,1), lty=2, lwd=1)
+  par(origpar)
+}
+
 plot.SSB.F.trend<-function(mod, alpha = 0.05)
 {
   origpar <- par(no.readonly = TRUE)
diff --git a/man/prepare_wham_input.Rd b/man/prepare_wham_input.Rd
index 4cd4d51b..192a94ae 100644
--- a/man/prepare_wham_input.Rd
+++ b/man/prepare_wham_input.Rd
@@ -12,6 +12,7 @@ prepare_wham_input(
   selectivity = NULL,
   M = NULL,
   NAA_re = NULL,
+  catchability = NULL,
   age_comp = NULL,
   basic_info = NULL
 )
@@ -31,6 +32,8 @@ prepare_wham_input(
 
 \item{NAA_re}{(optional) list specifying options for random effect on numbers-at-age, initial numbers at age, and recruitment model (see details)}
 
+\item{catchability}{(optional) list specifying options for priors and random effects on catchability (see details)}
+
 \item{age_comp}{(optional) character or named list, specifies age composition model for fleet(s) and indices (see details)}
 
 \item{basic_info}{(optional) list of basic population information for use when asap3 is not provided (see details)}
@@ -87,12 +90,18 @@ Otherwise, it must be a named list with the following components:
     \item{$year}{Years corresponding to observations (vector of same length as \code{$mean} and \code{$logsigma})}
     \item{$use_obs}{T/F (or 0/1) vector/matrix of the same dimension as \code{$mean} and \code{$logsigma}.
     Use the observation? Can be used to ignore subsets of the ecov without changing data files.}
-    \item{$lag}{Offset between the ecov observations and their affect on the stock.
-    I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}.}
+    \item{$lag}{integer vector of offsets between the ecov observation/process and their affect on the stock.
+    I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}. May also be a list (length=n_Ecov) of vectors (lenght = 2+n_indices) if multiple effects of one or more Ecovs are modeled.}
     \item{$process_model}{Process model for the ecov time-series. \code{"rw"} = random walk, \code{"ar1"} = 1st order autoregressive, \code{NA} = do not fit}
-    \item{$where}{Where does the ecov affect the population? \code{"recuit"} = recruitment,
-    \code{"M"} = natural mortality, \code{"none"} = fit ecov but without an effect on the population.}
-    \item{$how}{How does the ecov affect the \code{$where} process? These options are
+    \item{$where}{Character vector for where each ecov affects the population. \code{"recuit"} = recruitment,
+    \code{"M"} = natural mortality, \code{"q"} = catchability for indices, \code{"none"} = fit ecov but without an effect on the 
+    population. May also be a list (element for each ecov) of multiple character vectors so each ecov can multiple effects.}
+    \item{$indices}{indices that each ecov affects. Must be a list of length n.ecov, where each element is a vector of indices (1:n_indices). Must be provided when any of \code{where} = "q"}
+    \item{$link_model}{vector or (orthogonal polynomial order) models for effects or each ecov on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
+    where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b_0 + b_1*ecov + b_2*ecov^2 + ...}). A list (length = n_Ecov) of character vectors for modeling
+     effects on recruitment, M, and catchabilities for index 1,..., index n_indices (length of the vector is 2 + n_indices).}
+    \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
+    \item{$how}{How does the ecov affect the \code{$where} process? Currently this only isThese options are
     specific to the \code{$where} process.
 Options for recruitment are described in \href{https://www.sciencedirect.com/science/article/pii/S1385110197000221}{Iles & Beverton (1998)}:
   \describe{
@@ -107,10 +116,12 @@ Options for M:
   \describe{
     \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
     \item{= 1}{effect on mean M (shared across ages)}
+  }
+Options for catchability:
+  \describe{
+    \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
+    \item{= 1}{effect on one or more catchabilities (see \code{$indices)})}
   }}
-    \item{$link_model}{Model describing ecov effect on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
-    where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b0 + b1*ecov + b2*ecov^2 + ...}).}
-    \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
   }
 
 \code{selectivity} specifies options for selectivity, to overwrite existing options specified in the ASAP data file.
@@ -215,6 +226,28 @@ To fit a state-space model, specify \code{NAA_re}, a list with the following ent
 	 }
   }
 
+\code{catchability} specifies options for catchability. If \code{NULL} and \code{asap3} is not NULL, a single catchability parameter for each index is used with initial values specified in ASAP file. If both are NULL, initial catchabilities for all indices = 0.3.
+Otherwise, it is a list with the following optional entries:
+  \describe{
+    \item{$re}{Time-varying (random effects) for each index. Vector with length = number of indices.
+                 If \code{NULL}, catchability parameters are constant over time.
+                 Each entry of \code{catchability$re} must be one of the following options:
+                 \describe{
+                   \item{"none"}{(default) are constant}
+                   \item{"iid"}{vary by year and age/par, but uncorrelated}
+                   \item{"ar1"}{correlated by year (AR1)}
+                 }
+                }
+    \item{$initial_q}{Initial catchabilities for each index. vector length = number of indices. Will override values provided in \code{basic_info$q}. 
+       If \code{basic_info$q} and \code{asap3} are not provided default q values are 0.3.}
+    \item{$q_lower}{Lower bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 0.}
+    \item{$q_upper}{Upper bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 1000.}
+    \item{$fix_q}{Which catcabilities not estimated. vector of integers indicating which indices have catchability fixed. E.g. c(2,3) when the second and third catchabilties are fixed.}
+    \item{$prior_sd}{vector of NULL and standard deviations to use for gaussian prior on logit transform of catchability parameter. Length = number of indices. 
+      Indices with non-NULL values will have mean logit q as a random effect with mean determined by logit transform of \code{catchability$initial_q} and  
+      sigma as standard error.}
+  }
+
 \code{age_comp} specifies the age composition models for fleet(s) and indices.
 If \code{NULL}, the multinomial is used because this was the only option in ASAP.
 The age composition models available are:
@@ -270,7 +303,7 @@ The current options are:
     \item{$percentFMSY}{(0-100) percentage of Fmsy to use in projections.}
     \item{$XSPR_R_avg_yrs}{which years to average recruitments for calculating SPR-based SSB reference points.}
     \item{$XSPR_R_opt}{1(3): use annual R estimates(predictions) for annual SSB_XSPR, 2(4): use average R estimates(predictions).}
-    \item{$simulate_process_error}{T/F vector (length = 4). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, and Ecov. Only used if applicable.}
+    \item{$simulate_process_error}{T/F vector (length = 5). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, Ecov, and q. Only used if applicable.}
     \item{$simulate_observation_error}{T/F vector (length = 3). When simulating from the model, whether to simulate  catch, index, and ecov observations.}
     \item{$simulate_period}{T/F vector (length = 2). When simulating from the model, whether to simulate base period (model years) and projection period.}
     \item{$bias_correct_process}{T/F. Perform bias correction of log-normal random effects for NAA.}
diff --git a/src/helper_functions.hpp b/src/helper_functions.hpp
index cc2423b1..8420de9a 100644
--- a/src/helper_functions.hpp
+++ b/src/helper_functions.hpp
@@ -135,15 +135,9 @@ Type ddirmultinom_osa(vector<Type> obs, vector<Type> p,  Type phi, int do_log, v
 template<class Type>
 vector<Type> rdirmultinom(Type N, vector<Type> p, Type phi) //dirichlet generated from iid gammas
 {
-  int Nint = CppAD::Integer(N);
-  int dim = p.size();
-  vector<Type> obs(dim);
-  obs.setZero();
-  for(int i = 0; i < Nint; i++)
-  {
-    vector<Type> dp = rdirichlet(p, phi);
-    obs = obs + rmultinom(Type(1),dp);
-  }
+  //int Nint = CppAD::Integer(N);
+  vector<Type> dp = rdirichlet(p, phi);
+  vector<Type> obs = rmultinom(N,dp);
   return(obs);
 }
 
@@ -1178,7 +1172,7 @@ matrix<Type> poly_trans(vector<Type> x, int degree, int n_years_model, int n_yea
 
 template <class Type>
 Type get_pred_recruit_y(int y, int recruit_model, vector<Type> mean_rec_pars, vector<Type> SSB, matrix<Type> NAA, vector<Type> log_SR_a, 
-  vector<Type> log_SR_b, vector<int> Ecov_where, vector<int> Ecov_how, vector<matrix<Type> > Ecov_lm){
+  vector<Type> log_SR_b, matrix<int> Ecov_where, vector<int> Ecov_how, array<Type> Ecov_lm){
 
   /*
    * y: year (between 1 and n_years_model+n_years_proj)
@@ -1188,9 +1182,9 @@ Type get_pred_recruit_y(int y, int recruit_model, vector<Type> mean_rec_pars, ve
    * NAA: matrix of numbers at age
    * log_SR_a: yearly "a" parameters for SR function
    * log_SR_b: yearly "b" parameters for SR function
-   * Ecov_where: integer determining if Ecov is affecting recruitment
+   * Ecov_where: matrix of 0/1 with first column determining if Ecov is affecting recruitment
    * Ecov_how: integer vector with an element that tells how the Ecov is affecting recruitment
-   * Ecov_lm: matrix that holds linear predictor for Ecov
+   * Ecov_lm: array that holds linear predictor for Ecov
    */
   //recruit_model == 1, random walk
   Type pred_recruit = NAA(y-1,0);
@@ -1203,9 +1197,9 @@ Type get_pred_recruit_y(int y, int recruit_model, vector<Type> mean_rec_pars, ve
     if(recruit_model == 2) // random about mean
     {
       pred_recruit = exp(mean_rec_pars(0));
-      int nE = Ecov_where.size();
+      int nE = Ecov_where.rows();
       for(int i=0; i < nE; i++){
-        if(Ecov_where(i) == 1) if(Ecov_how(i) == 1) pred_recruit *= exp(Ecov_lm(i)(y,0));
+        if(Ecov_where(i,0) == 1) if(Ecov_how(i) == 1) pred_recruit *= exp(Ecov_lm(i,0,y,0));
       }
       //pred_NAA(y,0) = exp(mean_rec_pars(0));
       //if(Ecov_recruit > 0) if(Ecov_how(Ecov_recruit-1) == 1) pred_NAA(y,0) *= exp(Ecov_lm(y,Ecov_recruit-1));
@@ -1229,7 +1223,7 @@ Type get_pred_recruit_y(int y, int recruit_model, vector<Type> mean_rec_pars, ve
 
 template <class Type>
 vector<Type> get_pred_NAA_y(int y, int recruit_model, vector<Type> mean_rec_pars, vector<Type> SSB, matrix<Type> NAA, vector<Type> log_SR_a, 
-  vector<Type> log_SR_b, vector<int> Ecov_where, vector<int> Ecov_how, vector<matrix<Type> > Ecov_lm, matrix<Type> ZAA){
+  vector<Type> log_SR_b, matrix<int> Ecov_where, vector<int> Ecov_how, array<Type> Ecov_lm, matrix<Type> ZAA){
 
   /*
    * y: year (between 1 and n_years_model+n_years_proj)
@@ -1239,9 +1233,9 @@ vector<Type> get_pred_NAA_y(int y, int recruit_model, vector<Type> mean_rec_pars
    * NAA: matrix of numbers at age
    * log_SR_a: yearly "a" parameters for SR function
    * log_SR_b: yearly "b" parameters for SR function
-   * Ecov_where: integer vector determining if Ecov i affects recruitment (= 1)
+   * Ecov_where: matrix of 0/1 integer with first column determining if Ecov i affects recruitment (= 1)
    * Ecov_how: integer vector with an element that tells how the Ecov is affecting recruitment
-   * Ecov_lm: matrix that holds linear predictor for Ecov
+   * Ecov_lm: array that holds linear predictor for Ecov
    * ZAA: matrix of total mortality rate by year and age
    */
   int n_ages = NAA.cols();
@@ -1383,7 +1377,7 @@ matrix<Type> get_F_proj(int y, int n_fleets, vector<int> proj_F_opt, array<Type>
 
 template <class Type>
 matrix<Type> sim_pop(array<Type> NAA_devs, int recruit_model, vector<Type> mean_rec_pars, vector<Type> SSBin, matrix<Type> NAAin, vector<Type> log_SR_a, 
-  vector<Type> log_SR_b, vector<int> Ecov_where, vector<int> Ecov_how, vector<matrix<Type> > Ecov_lm, int n_NAA_sigma, 
+  vector<Type> log_SR_b, matrix<int> Ecov_where, vector<int> Ecov_how, array<Type> Ecov_lm, int n_NAA_sigma, 
   int do_proj, vector<int> proj_F_opt, array<Type> FAA, matrix<Type> FAA_tot, matrix<Type> MAA, matrix<Type> mature, array<Type> waa, 
   int waa_pointer_totcatch, int waa_pointer_ssb, vector<Type> fracyr_SSB, vector<Type> log_SPR0, vector<int> avg_years_ind, 
   int n_years_model, int n_fleets, vector<int> which_F_age, Type percentSPR, vector<Type> proj_Fcatch, Type percentFXSPR, vector<Type> F_proj_init, Type percentFMSY){
diff --git a/src/wham_v0.cpp b/src/wham_v0.cpp
index dc33bee8..bb9a29e1 100644
--- a/src/wham_v0.cpp
+++ b/src/wham_v0.cpp
@@ -57,7 +57,7 @@ Type objective_function<Type>::operator() ()
   DATA_VECTOR(q_upper);
   DATA_IVECTOR(use_q_prior);
   DATA_VECTOR(logit_q_prior_sigma);
-  DATA_IVECTOR(q_re_model); //n_indices, 0= no re, >0 = use re 
+  DATA_IVECTOR(use_q_re); //n_indices, 0= no re, >0 = use re 
   DATA_MATRIX(selpars_lower);
   DATA_MATRIX(selpars_upper);
   DATA_INTEGER(n_NAA_sigma); // 0 = SCAA, 1 = logR only, 2 = full state-space with shared sig_a for a > 1
@@ -79,7 +79,6 @@ Type objective_function<Type>::operator() ()
   DATA_INTEGER(bias_correct_pe); //bias correct lognormal process error?
   DATA_INTEGER(bias_correct_oe); //bias correct lognormal observation error?
   DATA_IVECTOR(Fbar_ages);
-  //DATA_INTEGER(simulate_state); //if 1 then state parameters will be simulated
   DATA_IVECTOR(simulate_state); //vector (0/1) if 1 then state parameters (NAA, MAA, sel, Ecov, q) in that order) will be simulated.
   DATA_IVECTOR(simulate_data); //vector (0/1) if 1 then data type (catch, indices, Ecov obs) will be simulated.
   DATA_IVECTOR(simulate_period); //vector (0/1) if 1 then period (model years, projection years) will be simulated.
@@ -105,16 +104,17 @@ Type objective_function<Type>::operator() ()
   DATA_INTEGER(n_years_Ecov); // num years in Ecov  process model
   DATA_IMATRIX(Ecov_use_obs); // all 0 if no Ecov
   DATA_MATRIX(Ecov_obs);
-  DATA_IVECTOR(Ecov_lag);
-  DATA_IVECTOR(Ecov_how); // 0 = no effect, 1 = controlling, 2 = limiting, 3 = lethal, 4 = masking, 5 = directive
-  DATA_IVECTOR(Ecov_poly); // polynomial order for ecov effects (1 = linear, 2 = quadratic, 3 = cubic, ...)
-  DATA_IVECTOR(Ecov_where); // 0 = no Ecov, 1 = recruit, 2 = mortality, 3 = q
+  //Below is not used anymore
+  //DATA_IVECTOR(Ecov_lag);
+  DATA_IVECTOR(Ecov_how); // specific to recruitment effects. 0 = no effect, 1 = controlling, 2 = limiting, 3 = lethal, 4 = masking, 5 = directive
+  //Below is not used anymore
+  //DATA_IMATRIX(Ecov_poly); // n_Ecov x 2+n_indices. polynomial order for ecov effects (1 = linear, 2 = quadratic, 3 = cubic, ...)
+  DATA_IMATRIX(Ecov_where); // n_Ecov x 2+n_indices. 0/1 values with columns corresponding to recruit, mortality, indices in that order
   DATA_IVECTOR(Ecov_model); // 0 = no Ecov, 1 = RW, 2 = AR1
-  DATA_IVECTOR(Ecov_for_q); // (n_indices) 0 = no Ecov, else which Ecov to use
   DATA_INTEGER(year1_Ecov); // first year Ecov
   DATA_INTEGER(year1_model); // first year model
-  DATA_IVECTOR(ind_Ecov_out_start); // index of Ecov_x to use for Ecov_out (operates on pop model, lagged)
-  DATA_IVECTOR(ind_Ecov_out_end); // index of Ecov_x to use for Ecov_out (operates on pop model, lagged)
+  DATA_IMATRIX(ind_Ecov_out_start); // n_Ecov x (2 + n_indices) index of Ecov_x to use for Ecov_out (operates on pop model, lagged effects specific the multiple types of effects each Ecov can have)
+  DATA_IMATRIX(ind_Ecov_out_end); // n_Ecov x (2 + n_indices) index of Ecov_x to use for Ecov_out (operates on pop model, lagged effects specific the multiple types of effects each Ecov can have)
   DATA_INTEGER(Ecov_obs_sigma_opt); // 1 = given, 2 = estimate 1 value, shared among obs, 3 = estimate for each obs, 4 = estimate for each obs as random effects
   DATA_IMATRIX(Ecov_use_re); // 0/1: use Ecov_re? If yes, add to nll. (n_years_Ecov + n_years_proj_Ecov) x n_Ecov
 
@@ -158,7 +158,7 @@ Type objective_function<Type>::operator() ()
 
   // parameters - environmental covariate ("Ecov")
   PARAMETER_MATRIX(Ecov_re); // nrows = n_years_Ecov, ncol = N_Ecov
-  PARAMETER_ARRAY(Ecov_beta); // dim = n.poly x n.ecov x n.ages, beta_R in eqns 4-5, Miller et al. (2016)
+  PARAMETER_ARRAY(Ecov_beta); // dim = (2 + n_indices) x n_poly x n_ecov x n_ages, beta_R in eqns 4-5, Miller et al. (2016)
   PARAMETER_MATRIX(Ecov_process_pars); // nrows = RW: 2 par (Ecov1, sig), AR1: 3 par (mu, sig, phi); ncol = N_ecov
   PARAMETER_MATRIX(Ecov_obs_logsigma); // options: given (data), fixed effect(s), or random effects
   PARAMETER_MATRIX(Ecov_obs_sigma_par); // ncol = N_Ecov, nrows = 2 (mean, sigma of random effects)
@@ -186,7 +186,7 @@ Type objective_function<Type>::operator() ()
   matrix<Type> ZAA(n_years_model + n_years_proj,n_ages);
   array<Type> QAA(n_years_model+n_years_proj,n_indices,n_ages);
   vector<matrix<Type> > selAA(n_selblocks); // selAA(b)(y,a) gives selectivity by block, year, age; selAA(b) is matrix with dim = n_years x n_ages;
-  vector<Type> q(n_indices);
+  matrix<Type> q(n_years_model+n_years_proj,n_indices);
   vector<Type> t_paa(n_ages);
   vector<Type> t_pred_paa(n_ages);
   int n_toavg = avg_years_ind.size();
@@ -312,77 +312,76 @@ Type objective_function<Type>::operator() ()
 
   // Environmental covariate process model --------------------------------------
   matrix<Type> Ecov_x(n_years_Ecov + n_years_proj_Ecov, n_Ecov); // 'true' estimated Ecov (x_t in Miller et al. 2016 CJFAS)
-  matrix<Type> Ecov_out(n_years_model + n_years_proj, n_Ecov); // Pop model uses Ecov_out(t) for processes in year t (Ecov_x shifted by lag and padded)
   matrix<Type> nll_Ecov(n_years_Ecov + n_years_proj_Ecov, n_Ecov); // nll contribution each Ecov_re
   nll_Ecov.setZero();
 
-  if(Ecov_model(0) == 0){ // no Ecov
-    Ecov_out.setZero(); // set Ecov_out = 0
-  } else { // yes Ecov
-    for(int i = 0; i < n_Ecov; i++){ // loop over Ecovs
-      // Ecov model option 1: RW
-      if(Ecov_model(i) == 1){
-        Type Ecov_sig; // sd (sig_x in Eq1, pg 1262, Miller et al. 2016)
-        Ecov_sig = exp(Ecov_process_pars(1,i));
-        Type Ecov1; // Ecov_x in year 1 (fixed effect)
-        Ecov1 = Ecov_process_pars(0,i);
+// These next two lines don't allow for Ecov_model >0 after the first column of Ecov.
+//  if(Ecov_model(0) == 0){ // no Ecov
+//  } else { // yes Ecov
+  for(int i = 0; i < n_Ecov; i++){ // loop over Ecovs
+    // Ecov model option 1: RW
+    if(Ecov_model(i) == 1){
+      Type Ecov_sig; // sd (sig_x in Eq1, pg 1262, Miller et al. 2016)
+      Ecov_sig = exp(Ecov_process_pars(1,i));
+      Type Ecov1; // Ecov_x in year 1 (fixed effect)
+      Ecov1 = Ecov_process_pars(0,i);
 
-        Ecov_x(0,i) = Ecov1;
-        nll_Ecov(1,i) -= dnorm(Ecov_re(1,i), Ecov1, Ecov_sig, 1); // Ecov_re(0,i) set to NA
-        SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(1,i) == 1) {
-          if(simulate_period(0) == 1) {
-            Ecov_re(1,i) = rnorm(Ecov1, Ecov_sig);
-          }
+      Ecov_x(0,i) = Ecov1;
+      nll_Ecov(1,i) -= dnorm(Ecov_re(1,i), Ecov1, Ecov_sig, 1); // Ecov_re(0,i) set to NA
+      SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(1,i) == 1) {
+        if(simulate_period(0) == 1) {
+          Ecov_re(1,i) = rnorm(Ecov1, Ecov_sig);
         }
-        Ecov_x(1,i) = Ecov_re(1,i);
-        for(int y = 2; y < n_years_Ecov + n_years_proj_Ecov; y++){
-          nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_re(y-1,i), Ecov_sig, 1);
-          SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(y,i) == 1) {
-            if((simulate_period(0) == 1 & y < n_years_Ecov) | (simulate_period(1) == 1 & y > n_years_Ecov-1)) {
-              Ecov_re(y,i) = rnorm(Ecov_re(y-1,i), Ecov_sig);
-            }
+      }
+      Ecov_x(1,i) = Ecov_re(1,i);
+      for(int y = 2; y < n_years_Ecov + n_years_proj_Ecov; y++){
+        nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_re(y-1,i), Ecov_sig, 1);
+        SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(y,i) == 1) {
+          if((simulate_period(0) == 1 & y < n_years_Ecov) | (simulate_period(1) == 1 & y > n_years_Ecov-1)) {
+            Ecov_re(y,i) = rnorm(Ecov_re(y-1,i), Ecov_sig);
           }
-          Ecov_x(y,i) = Ecov_re(y,i);
         }
+        Ecov_x(y,i) = Ecov_re(y,i);
       }
+    }
 
-      // Ecov model option 2: AR1
-      if(Ecov_model(i) == 2){
-        Type Ecov_mu; // mean
-        Type Ecov_phi; // autocorrelation
-        Type Ecov_sig; // conditional sd
-        Ecov_mu = Ecov_process_pars(0,i);
-        Ecov_phi = -Type(1) + Type(2)/(Type(1) + exp(-Ecov_process_pars(2,i)));
-        Ecov_sig = exp(Ecov_process_pars(1,i));
+    // Ecov model option 2: AR1
+    if(Ecov_model(i) == 2){
+      Type Ecov_mu; // mean
+      Type Ecov_phi; // autocorrelation
+      Type Ecov_sig; // conditional sd
+      Ecov_mu = Ecov_process_pars(0,i);
+      Ecov_phi = -Type(1) + Type(2)/(Type(1) + exp(-Ecov_process_pars(2,i)));
+      Ecov_sig = exp(Ecov_process_pars(1,i));
 
-        nll_Ecov(0,i) -= dnorm(Ecov_re(0,i), Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))), 1);
-        SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(0,i) == 1) {
-          if(simulate_period(0) == 1) {
-            Ecov_re(0,i) = rnorm(Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))));
-          }
+      nll_Ecov(0,i) -= dnorm(Ecov_re(0,i), Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))), 1);
+      SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(0,i) == 1) {
+        if(simulate_period(0) == 1) {
+          Ecov_re(0,i) = rnorm(Type(0), Ecov_sig*exp(-Type(0.5) * log(Type(1) - pow(Ecov_phi,Type(2)))));
         }
-        for(int y = 1; y < n_years_Ecov + n_years_proj_Ecov; y++)
-        {
-          //FIXED on q_model branch 
-          nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_phi * Ecov_re(y-1,i), Ecov_sig, 1);
-          SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(y,i) == 1) {
-            if((simulate_period(0) == 1 & y < n_years_Ecov) | (simulate_period(1) == 1 & y > n_years_Ecov-1)) {
-              Ecov_re(y,i) = rnorm(Ecov_phi * Ecov_re(y-1,i), Ecov_sig);
-            }
+      }
+      for(int y = 1; y < n_years_Ecov + n_years_proj_Ecov; y++)
+      {
+        //FIXED on q_model branch 
+        nll_Ecov(y,i) -= dnorm(Ecov_re(y,i), Ecov_phi * Ecov_re(y-1,i), Ecov_sig, 1);
+        SIMULATE if(simulate_state(3) == 1 & Ecov_use_re(y,i) == 1) {
+          if((simulate_period(0) == 1 & y < n_years_Ecov) | (simulate_period(1) == 1 & y > n_years_Ecov-1)) {
+            Ecov_re(y,i) = rnorm(Ecov_phi * Ecov_re(y-1,i), Ecov_sig);
           }
         }
-        for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++) Ecov_x(y,i) = Ecov_mu + Ecov_re(y,i);
       }
+      for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++) Ecov_x(y,i) = Ecov_mu + Ecov_re(y,i);
+    }
 
-      // add to nll if estimated (option in projection years to fix Ecov at last or average value)
-      for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++){
-        if(Ecov_use_re(y,i) == 1){
-          nll += nll_Ecov(y,i);
-        }
+    // add to nll if estimated (option in projection years to fix Ecov at last or average value)
+    for(int y = 0; y < n_years_Ecov + n_years_proj_Ecov; y++){
+      if(Ecov_use_re(y,i) == 1){
+        nll += nll_Ecov(y,i);
       }
-    } // end loop over Ecovs
-    SIMULATE if(simulate_state(3) == 1) if(sum(simulate_period) > 0) REPORT(Ecov_re);
-  }
+    }
+  } // end loop over Ecovs
+  SIMULATE if(simulate_state(3) == 1) if(sum(simulate_period) > 0) REPORT(Ecov_re);
+//  }
 
   // Environmental covariate observation model -------------------------------------
   //TODO: Ecov obs are not yet simulated in projection years!!!!!!!!
@@ -417,37 +416,50 @@ Type objective_function<Type>::operator() ()
     REPORT(Ecov_obs_logsigma);
   }
 
+ 
   // Lag environmental covariates -------------------------------------
   // Then use Ecov_out(t) for processes in year t, instead of Ecov_x
+  int n_effects = Ecov_beta.dim(0); // 2 + n_indices (recruitment, mortality and any catchabilities)
+  array<Type> Ecov_out(n_years_model + n_years_proj, n_effects, n_Ecov); // Pop model uses Ecov_out(t) for processes in year t (Ecov_x shifted by lag and padded)
+  Ecov_out.setZero(); // set Ecov_out = 0
   for(int i = 0; i < n_Ecov; i++){
-    int ct = 0;
-    for(int y = ind_Ecov_out_start(i); y < ind_Ecov_out_end(i) + 1 + n_years_proj; y++){
-      Ecov_out(ct,i) = Ecov_x(y,i);
-      ct++;
+    for(int t = 0; t < n_effects; t++){
+      int ct = 0;
+      for(int y = ind_Ecov_out_start(i,t); y < ind_Ecov_out_end(i,t) + 1 + n_years_proj; y++){
+        Ecov_out(ct,t,i) = Ecov_x(y,i);
+        ct++;
+      }
     }
   }
 
   // Calculate ecov link model (b1*ecov + b2*ecov^2 + ...) --------------------
-  vector<matrix<Type> > Ecov_lm(n_Ecov); // ecov linear model for each Ecov, dim = n_years_model + n_years_proj, n_ages
+  // ecov_beta is now 4D array, dim = (2 + n_indices) x n_poly x n_ecov x n_ages
+  int n_poly = Ecov_beta.dim(1); // now a 4D array dim: (n_effects,n_poly,n_Ecov,n_ages) is second dimension
+  //vector<matrix<Type>> Ecov_lm(n_Ecov)(n_effects); // ecov linear model for each Ecov, dim = n_years_model + n_years_proj, n_ages
   // Ecov_lm.setZero();
-  // ecov_beta is now 3D array, dim = n.poly x n.ecov x n.ages
-  int n_poly = Ecov_beta.rows(); // still works on 3D array bc npoly is first dimension
+  // Ecov_lm stores the linear models for each Ecov and where it is used. dim = n_Ecov, n_effects, n_years_model + n_years_proj, n_ages
+  // n_effects dimension is: 0: recruitment, 1: M, 2-1+n_indices: which catchability it affects
+  array<Type> Ecov_lm(n_Ecov, n_effects,n_years_model + n_years_proj, n_ages); 
+  //vector<matrix<Type> > Ecov_lm(n_Ecov); // ecov linear model for each Ecov, dim = n_years_model + n_years_proj, n_ages
   for(int i = 0; i < n_Ecov; i++){
-    vector<Type> thecol = Ecov_out.col(i);
-    matrix<Type> X_poly(n_years_model + n_years_proj, n_poly);
-    X_poly.setZero();
-    if(n_poly == 1){ // n_poly = 1 if ecov effect is none or linear
-      X_poly = thecol.matrix();
-    } else { // n_poly > 1, get poly transformation for ith ecov
-      X_poly = poly_trans(thecol, n_poly, n_years_model, n_years_proj);
-    }
-    matrix<Type> tmp(n_years_model + n_years_proj, n_ages);
-    tmp.setZero();
-    Ecov_lm(i) = tmp;
-    for(int y = 0; y < n_years_model + n_years_proj; y++){
-      for(int a = 0; a < n_ages; a++){
-        for(int j = 0; j < n_poly; j++){
-          Ecov_lm(i)(y,a) += Ecov_beta(j,i,a) * X_poly(y,j); // poly transformation returns design matrix, don't need to take powers
+    for(int t = 0; t < n_effects; t++){
+      vector<Type> thecol(n_years_model + n_years_proj);
+      for(int y = 0; y < n_years_model + n_years_proj; y++) thecol(y) = Ecov_out(y,t,i);
+      matrix<Type> X_poly(n_years_model + n_years_proj, n_poly);
+      X_poly.setZero();
+      if(n_poly == 1){ // n_poly = 1 if ecov effect is none or linear
+        X_poly = thecol.matrix();
+      } else { // n_poly > 1, get poly transformation for ith ecov
+        X_poly = poly_trans(thecol, n_poly, n_years_model, n_years_proj);
+      }
+      //matrix<Type> tmp(n_years_model + n_years_proj, n_ages);
+      //tmp.setZero();
+      //Ecov_lm(i)(t) = tmp;
+      for(int y = 0; y < n_years_model + n_years_proj; y++){
+        for(int a = 0; a < n_ages; a++){
+          for(int j = 0; j < n_poly; j++){
+            Ecov_lm(i,t,y,a) += Ecov_beta(t,j,i,a) * X_poly(y,j); // poly transformation returns design matrix, don't need to take powers
+          }
         }
       }
     }
@@ -560,9 +572,9 @@ Type objective_function<Type>::operator() ()
   }
   // add ecov effect on M (by year, shared across ages)
   for(int i=0; i < n_Ecov; i++){
-    if(Ecov_where(i) == 2) if(Ecov_how(i) == 1){ // if ecov i affects M
+    if(Ecov_where(i,1) == 1) { //not sure why Ecov_how is needed for M. if(Ecov_how(i) == 1){ // if ecov i affects M
       for(int a = 0; a < n_ages; a++){
-        for(int y = 0; y < n_years_model + n_years_proj; y++) MAA(y,a) *= exp(Ecov_lm(i)(y,a));
+        for(int y = 0; y < n_years_model + n_years_proj; y++) MAA(y,a) *= exp(Ecov_lm(i,1,y,a));
       }
     }
   }
@@ -597,17 +609,17 @@ Type objective_function<Type>::operator() ()
   logit_q_mat.setZero();
   for(int i = 0; i < n_indices; i++) {
     
-    //use prior for q? q_prior_re are random effects with mean logit_q and sd = logit_q_prior_sigma.
+    //use prior for q? q_prior_re are random effects with mean logit_q (fixed) and sd = logit_q_prior_sigma.
     if(use_q_prior(i) == 1){ 
       nll_q_prior(i) -= dnorm(q_prior_re(i), logit_q(i), logit_q_prior_sigma(i), 1);
       SIMULATE if(simulate_state(4) == 1) if(sum(simulate_period) > 0){
         q_prior_re(i) = rnorm(logit_q(i), logit_q_prior_sigma(i));
       } 
-      for(int y = 1; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += q_prior_re(i);
+      for(int y = 0; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += q_prior_re(i);
     }
-    else for(int y = 1; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += logit_q(i);
+    else for(int y = 0; y < n_years_model + n_years_proj; y++) logit_q_mat(y,i) += logit_q(i);
     
-    if(q_re_model(i) > 0) // random effects on q, q_re = AR1 deviations on (year,age), dim = n_years x n_M_a
+    if(use_q_re(i) > 0) // random effects on q, q_re = AR1 deviations on (year,age), dim = n_years x n_M_a
     {
       sigma_q(i) = exp(q_repars(i,0)); // conditional sd
       rho_q(i) = rho_trans(q_repars(i,1)); // autocorrelation
@@ -616,6 +628,7 @@ Type objective_function<Type>::operator() ()
       SIMULATE if(simulate_state(4) == 1 & simulate_period(0) == 1) {
         q_re(0,i) = rnorm(Type(0), sigma_q(i)*exp(-0.5 * log(1 - pow(rho_q(i),Type(2)))));
       }
+      logit_q_mat(0,i) += q_re(0,i); //add in q random effects.
       for(int y = 1; y < n_years_model + n_years_proj; y++)
       {
         nll_q(y,i) -= dnorm(q_re(y,i), rho_q(i) * q_re(y-1,i), sigma_q(i), 1);
@@ -629,20 +642,13 @@ Type objective_function<Type>::operator() ()
     }
   }
 
-  REPORT(logit_q_mat);
-  REPORT(sigma_q);
-  REPORT(rho_q);
-  REPORT(nll_q);
-  REPORT(nll_q_prior);
-  REPORT(q_prior_re); //even if q_prior_re not simulated
-  REPORT(q_re); //even if q_re not simulated.
   nll += nll_q.sum() + nll_q_prior.sum();
   
   for(int y = 0; y < n_years_model + n_years_proj; y++) {
     for(int ind = 0; ind < n_indices; ind++) {
       for(int i=0; i < n_Ecov; i++){
-        if(Ecov_where(i) == 3) if(Ecov_for_q(ind) == i + 1){ // if ecov i affects q and which index
-          logit_q_mat(y,ind) += Ecov_lm(i)(y,0);
+        if(Ecov_where(i,1+ind) == 1){ // if(Ecov_for_q(ind) == i + 1){ // if ecov i affects q and which index
+          logit_q_mat(y,ind) += Ecov_lm(i,1+ind,y,0);
         }
       }
       q(y,ind) = q_lower(ind) + (q_upper(ind) - q_lower(ind))/(1 + exp(-logit_q_mat(y,ind)));
@@ -664,6 +670,16 @@ Type objective_function<Type>::operator() ()
       QAA(y,i,a) = q(y,i) * selAA(selblock_pointer_indices(n_years_model-1,i)-1)(n_years_model-1,a);
     }
   }
+  REPORT(logit_q_mat);
+  ADREPORT(logit_q_mat);
+  REPORT(sigma_q);
+  REPORT(rho_q);
+  REPORT(nll_q);
+  REPORT(nll_q_prior);
+  REPORT(q_prior_re); //even if q_prior_re not simulated
+  REPORT(q_re); //even if q_re not simulated.
+  REPORT(q);
+  REPORT(QAA);
 
   // Construct fishing mortality-at-age (FAA)
   FAA_tot.setZero();
@@ -748,18 +764,18 @@ Type objective_function<Type>::operator() ()
         log_SR_b.fill(mean_rec_pars(1));
       }
       for(int i=0; i < n_Ecov; i++){
-        if(Ecov_where(i) == 1){ // if ecov i affects recruitment
+        if(Ecov_where(i,0) == 1){ // if ecov i affects recruitment
           for(int y = 0; y < n_years_model + n_years_proj; y++)
           {
             // (1) "controlling" = dens-indep mortality or (4) "masking" = metabolic/growth (decreases dR/dS)
             if(Ecov_how(i) == 1 | Ecov_how(i) == 4)
             {
-              log_SR_a(y) += Ecov_lm(i)(y,0);
+              log_SR_a(y) += Ecov_lm(i,0,y,0);
             }
             // (2) "limiting" = carrying capacity or (4) "masking" = metabolic/growth (decreases dR/dS)
             if(Ecov_how(i) == 2 | Ecov_how(i) == 4)
             {
-              log_SR_b(y) += Ecov_lm(i)(y,0);
+              log_SR_b(y) += Ecov_lm(i,0,y,0);
             }
           }
         }
@@ -786,16 +802,16 @@ Type objective_function<Type>::operator() ()
         log_SR_b.fill(mean_rec_pars(1));
       }
       for(int i=0; i < n_Ecov; i++){
-        if(Ecov_where(i) == 1){ // if ecov i affects recruitment
+        if(Ecov_where(i,0) == 1){ // if ecov i affects recruitment
           for(int y = 0; y < n_years_model + n_years_proj; y++)
           {
             if(Ecov_how(i) == 1) // "controlling" = dens-indep mortality
             {
-              log_SR_a(y) += Ecov_lm(i)(y,0);
+              log_SR_a(y) += Ecov_lm(i,0,y,0);
             }
             if(Ecov_how(i) == 4) // "masking" = metabolic/growth (decreases dR/dS)
             { //NB: this is not identical to Iles and Beverton (1998), but their definition can give negative values of "b"
-              log_SR_b(y) += 1.0 + Ecov_lm(i)(y,0);
+              log_SR_b(y) += 1.0 + Ecov_lm(i,0,y,0);
             }
           }
         }
@@ -1279,8 +1295,6 @@ Type objective_function<Type>::operator() ()
   REPORT(SSB);
   REPORT(selAA);
   REPORT(MAA);
-  REPORT(q);
-  REPORT(QAA);
   REPORT(F);
   REPORT(FAA);
   REPORT(FAA_tot);