mainF <- function(v) {#This gets run once for every river cell
  #Set discharges lower than 0 to NA
  v[v<=minQ] <- NA

  df <- na.omit(data.table(discharge=v, year=tyears))
  rm(v)

  alloc.col(df, length(durations)+2) #allocate sufficient columns
  #assign rolling means in a loop
  for (i in durations) {
    set(df,
        j = paste0("D", i),
        value = roll_mean(df[["discharge"]], i, na.rm=TRUE, fill=NA))
  }
  df <- na.omit(df)

  dfsumm <- df[, lapply(.SD, max, na.rm = TRUE), by = year]
  #Now clean away infinites from dfsumm
  invisible(lapply(names(dfsumm),function(.name) set(dfsumm, which(is.infinite(dfsumm[[.name]])), j = .name,value =NA)))
  d <- dfsumm[["D1"]]
  gscale <- sd(d, na.rm=TRUE) * sqrt(6) / pi          #Alpha
  glocation <- mean(d, na.rm=TRUE) - emgamma * gscale #Mu
  fit <- fitdist(d, #default options: MLE estimation
               'gumbel',
               start = list(location = glocation, scale = gscale)
               )
  fitmu <- fit$estimate[1]
  fitalpha <- fit$estimate[2]
  eps<-NULL
  for (i in seq_along(durations)) {
    eps[i] <- mean(dfsumm[[paste0("D", durations[i])]], na.rm=TRUE) / mean(d, na.rm=TRUE)
  }

  #Fit for epsilon values (least squares method)
  fit <- nlsLM(epsilon ~ I((1 + fitbeta * D)^(-fitgamma)), #Function from NERC(1975)
            data = data.frame(epsilon=eps, D=durations),
            start = list(fitbeta = startbeta, fitgamma=startgamma),
            lower=c(lowerbeta, lowergamma),
            upper=c(upperbeta, uppergamma))
  #Get the coefficients of the fit
  fitbeta  <- coef(fit)[1]
  fitgamma <- coef(fit)[2]

  return(c(eps, fitmu, fitalpha, fitbeta, fitgamma))
}