#' @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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.