StatSeasAtlHurr <- function(atlano=NULL, tropano=NULL, hrvar='HR') {
  #
  # Verify that variables are either TC, HR or PDI.
  # -----------------------------------------------
  #
  if (hrvar != "HR" && hrvar != "TC" && hrvar != "PDI") {
    stop("Hurricane variable not recognized.")
  }
  #
  # Verify that both Atl and Trop SSTA are present.
  # -----------------------------------------------
  #
  if (is.null(atlano)) {
    stop("Atlantic SST missing.")
  }
  if (is.null(tropano)) {
    stop("Tropical SST missing.")
  }
  #
  # Verify that Atl and Trop SSTA are of the same dimensions.
  # ---------------------------------------------------------
  #
  if (length(dim(atlano)) != length(dim(tropano))) {
    stop("Input arrays are of different dimensions.")
  } else {
    for (i in 1:length(dim(atlano))) {
      if (dim(atlano)[i] != dim(tropano)[i]) {
        stop("Input arrays are of different sizes.")
      }
    }
  }
  #
  # Get the values of the betas according to the hurricane activity measure we specified.
  # -------------------------------------------------------------------------------------
  #
  if (hrvar == 'HR') {
    # beta's are derived from Villarini et al. (2012), Mon Wea Rev, 140, 44-65.
    # beta's are for corrected hurricane data + ERSST with SBC criteria (table 2)
    beta0    <-  1.85    
    betaAtl  <-  1.05 
    betaTrop <- -1.17 
  } else if (hrvar == 'TC') {
    # beta's are from Villarini et al. (2010), Mon Wea Rev, 138, 2681-2705. 
    # beta's are for corrected TC data (lifetime >= 48h) + ERSST (table 5) 
    beta0    <-  2.10
    betaAtl  <-  1.02
    betaTrop <- -1.05
  } else if (hrvar == 'PDI') {
    #  beta's are from Villarini et al. (2012), J Clim, 25, 625-637.
    #  beta's are from ERSST, with SBC penalty criterion (table 1)
    beta0    <-  0.76  
    betaAtl  <-  1.94  
    betaTrop <- -1.78
  } 
  #
  # Create matrix of similar dimension as atlano for beta0.
  # -------------------------------------------------------
  #
  intercept <- array(beta0, dim(atlano))
  #
  # Compute statistical relationship b/w SSTAs and mean hurricane activity. 
  # -----------------------------------------------------------------------
  #
  atl  <- betaAtl  * atlano
  trop <- betaTrop * tropano
  #
  temp <- intercept + atl + trop 
  #
  statval <- list(mean=array(NaN,dim(atl)),var=array(NaN,dim(atl)))
  statval$mean[] <- vapply(X=temp, FUN=exp, numeric(1))
  #
  # Compute the variance of the distribution.
  # TC and HR follow a Poisson distribution, so the variance is equal to the mean.
  # PDI follows a gamma distribution, with sigma = -0.57.  (variance = sigma^2 * mean^2).
  # ----------------------------------------------------------------------------------
  #
  if (hrvar == "HR" && hrvar == "TC") {
   statval$var <- statval$mean
  } else {
   sigma <- -0.57
   statval$var[] <- sigma^2 * vapply(X=statval$mean, FUN=function(x) x^2, numeric(1))
  }
  #  
  #  Output
  # ~~~~~~~~
  #
  statval
}