R/main.R

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

#' @title Main control function for runnning the model
#' @description Control function to run the ecosystem metabolism model using a defined data source.
#' Performs necessary unit conversions, calls the metabolism() function, and stroes results.
#' @importFrom zoo rollsum
#' @return Returns a data frame of model results over the full diel period.
#' @note Currently works with the installed Fishtrap Creek data set only.
#' @author Gordon W. Holtgrieve (gholt@uw.edu)
#' @export

main <- function(temperature, localHour, DOY, PAR = 'n',
                 TAlk, pCO2_init, pCO2_air = 410, O2_init,
                 alphaPI_hr, Rstd_day, kstd_hr, referenceT = 20,
                 lat, long, depth, salinity = 0, masl=0){


    # tm <<- temperature; l <<- localHour; d <<- DOY;
    # tt <<- TAlk; pp <<- pCO2_init; o <<- O2_init;  d <<- depth;
    # aa <<- alphaPI_hr; r <<- Rstd_day; k <<- kstd_hr; lt <<- lat; ln <<- long;
    # PAR <<- 'n'; pCO2_air <<- 410; referenceT <<- 20; salinity <<- 0; masl <<- 0
#
#     temperature = tm; localHour = l; DOY = d;
#     TAlk = tt; pCO2_init = pp; O2_init = o;
#     alphaPI_hr = aa; Rstd_day = r; kstd_hr = k
#     lat = lt; long = ln; depth = d
#
#     x.lt = as.POSIXlt(FishtrapCr$dateTime)
#     temperature = FishtrapCr$Temp_C
#     PAR = FishtrapCr$PAR_uE
#     localHour = x.lt$hour + x.lt$min/60 + x.lt$sec/3600
#     DOY = x.lt$yday
#     lat = 48.93861111
#     long = -122.47861111
#
#     tz = -8
#     masl = 1
#     salinity = 0


  # Temperature, localHour, and DOY must be of the same length.
  # TAlk is a constant (bogus assumption??)

  ####  Need to add O2 side ####

  #Determine time step
  timeStep <- as.numeric(diff.POSIXt(localHour)[1])  #Should be units of hour as a number
  nObs <- length(localHour)
  x <- rep(NA, nObs)

  pressure <- 1  # atm.  Need to add atmoshperic pressure ajustment (based on altitude) function later.

  if(any(PAR == "n")){
    PAR <- get.PAR(lat = lat, long = long, localHour = localHour, DOY = DOY, tz = tz, masl = masl)
  }

  output <- data.frame(localHour = localHour, temperature = temperature, PAR = PAR,
                 pCO2 = NA, CO2_ = NA, HCO3 = NA, CO3 = NA, TCO2 = NA, TAlk = NA,
                 pH = NA, O2 = NA, P = NA, R = NA, G_CO2 = NA, G_O2 = NA)

  System_init <- get.CO2Sys.TAlk_pCO2(TAlk = TAlk, pCO2 = pCO2_init,
                                      temperature = temperature[1],
                                      salinity = salinity,
                                      pressure = pressure)

  # Store in output data frame
  output[1, 4:10] = as.data.frame(System_init)
  output$O2[1] <- O2_init

  # Convert R, P, and k to units of per timestep
  alphaPI_m2 <- alphaPI_hr * timeStep
  Rstd_m2 <- Rstd_day * timeStep / 24
  kstd <- kstd_hr * timeStep

  # Convert alphaPI and Rstd from m-2 to m-3
  alphaPI_m3 <- alphaPI_m2 / depth
  Rstd_m3 <- Rstd_m2 / depth

  # Convert alphaPI and Rstd from m-3 to L
  alphaPI <- alphaPI_m3 / 1000
  Rstd_mol <- Rstd_m3 / 1000

  # Convert Rstd from mol to umol
  Rstd <- Rstd_mol * 10^6

  #
  for (i in 2:(nObs)) {
      oneTime <- metabolism(TCO2_start = output$TCO2[i-1], TAlk = TAlk,
            pCO2water = output$pCO2[i-1], pCO2air = pCO2_air,
            O2_start = output$O2[i-1],
            alphaPI = alphaPI, Rstd = Rstd, kstd = kstd, referenceT = referenceT,
            PAR = PAR[i-1], temperature = temperature[i-1],
            salinity = salinity, pressure = pressure, depth = depth)

      output[i, 4:ncol(output)] = as.data.frame(oneTime)
  }

  #Calculate 24hr GPP, ER, and G (units of umol/L/day)
  GPP_umolL <- mean(zoo::rollsum(output$P, k = 24/timeStep, align = "right", na.rm = T))
  ER_umolL <- mean(zoo::rollsum(output$R, k = 24/timeStep, align = "right", na.rm = T))
  G_CO2_24hr_umolL <- mean(zoo::rollsum(output$G_CO2, k = 24/timeStep, align = "right", na.rm = T))
  G_O2_24hr_umolL <- mean(zoo::rollsum(output$G_O2, k = 24/timeStep, align = "right", na.rm = T))

  # Convert from umol/L/d to umol/m3/d
  GPP_m3 <- GPP_umolL * 1000
  ER_m3 <- ER_umolL * 1000
  G_CO2_24hr_m3 <- G_CO2_24hr_umolL * 1000
  G_O2_24hr_m3 <- G_O2_24hr_umolL * 1000

  # Convert from umol/m3/d to umol/m2/d
  GPP_m2 <- GPP_m3 * depth
  ER_m2 <- ER_m3 * depth
  G_CO2_24hr_m2 <- G_CO2_24hr_m3 * depth
  G_O2_24hr_m2 <- G_O2_24hr_m3 * depth

  # Convert from umol/m2/d to mol/m2/d
  GPP <- GPP_m2 / 10^6
  ER <- ER_m2 / 10^6
  G_CO2_24hr <- G_CO2_24hr_m2 / 10^6
  G_O2_24hr <- G_O2_24hr_m2 / 10^6

  out <- list(timeSeries = output, GPP = GPP, ER = ER, G_CO2_24hr = G_CO2_24hr, G_O2_24hr = G_O2_24hr)

  return(out)
}