R/metabolism.R

In gholtgrieve/gassyPants: Functions Useful for Working With Gasses Dissolved in Water

#' @title Calculate ecosytem metabolism and its effect on dissolved O2 and CO2 pools
#' @description Calculates ecosystem metabolic rates for a given time step and applies them to dissolved O2 and
#' CO2 pools in a water body.  Rates are given in umol/L/time step.
#' @return  A numeric list of all carbonate system components, O2 concentration, & instantaneous P, R, G (O2 and CO2)
#' @note Function is incomplete. Needs runge-kutta 4th order approximation applied to dif eqs.
#' @author Gordon W. Holtgrieve (gholt@uw.edu)
#' @export

metabolism <- function(TCO2_start, TAlk, pCO2water, pCO2air, O2conc,
                       alphaPI, Rstd, kstd, referenceT, PAR, O2_start,
                       temperature, salinity, pressure, depth){

  # Error handling


 # Model
  P <- get.P(PAR = PAR, alphaPI = alphaPI)
  R <- get.R(temperature = temperature, Rstd = Rstd, referenceT = referenceT)
  G_CO2 <- get.G(temperature = temperature, salinity = salinity,
             pressure = pressure, pCO2water = pCO2water, gas = "CO2",
             pCO2air = pCO2air, kstd = kstd, referenceT = referenceT,
             depth = depth)
  G_O2 <- get.G(temperature = temperature, salinity = salinity, O2conc = O2_start,
                pressure = pressure, gas = "O2", kstd = kstd,
                referenceT = referenceT, depth = depth)

  # Apply CO2 mass balance
  # Units of umol/L/timestep
  # Needs R-K in the future
  dCO2dt <- R - P + G_CO2

  # Add/substract for TCO2
  TCO2_final <- TCO2_start + dCO2dt

  # Re-equilibrate the carbonate system. Returns a list of all carbonate parameters.
  output <- get.CO2Sys.TAlk_TCO2(TAlk = TAlk, TCO2 = TCO2_final,
                                 temperature = temperature,
                                 salinity = salinity,
                                 pressure = pressure)

  # Now do O2 mass balance
  # Units of umol/L/timestep
  # Needs R-K in the future
  dO2dt <- P - R + G_O2
  # print(O2_start); print(dO2dt)
  output$O2 <- O2_start + dO2dt

  output$P <- P
  output$R <- R
  output$G_CO2 <- G_CO2
  output$G_O2 <- G_O2

  return(output)
}