GLM examples

Warm Lake

No FABM

sim_folder <- "../release/warmlake/aed2"
# sim_folder <- "release/warmlake/aed2"
glm_nml_path <- file.path(sim_folder, "glm2.nml")
glm_nml <- read_nml(glm_nml_path)

names(glm_nml)

##  [1] "glm_setup"     "wq_setup"      "morphometry"   "time"         
##  [5] "output"        "init_profiles" "meteorology"   "bird_model"   
##  [9] "inflow"        "outflow"

(input_files <- c(glm_nml$meteorology$meteo_fl, 
                  glm_nml$inflow$inflow_fl, 
                  glm_nml$outflow$outflow_fl))

## [1] "met_hourly.csv"            "inflow_1.csv,inflow_2.csv"
## [3] "outflow.csv"

(start_end <- glm_nml$time[c("start", "stop")])

## $start
## [1] "1997-01-02 00:00:00"
## 
## $stop
## [1] "2001-09-30 00:00:00"

run_glm(sim_folder)

lake_csv <- read.csv(file.path(sim_folder, "lake.csv"), 
                     stringsAsFactors = FALSE)
lake_csv$time <- as.POSIXct(lake_csv$time)
lake_csv <- gather(lake_csv, variable, value, -time)
lake_csv <- filter(lake_csv, 
               variable %in% c("Tot.Inflow.Vol", "Tot.Outflow.Vol", 
                               "Rain", "Evaporation"))

ggplot(data = lake_csv) + 
  geom_line(aes(x = time, y = value, color = variable))

out_file <- file.path(sim_folder, "output.nc")
(s_vars <- sim_vars(out_file))

sim_metrics(with_nml = TRUE)

## Warning in layerD - Zcv: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.

##  [1] "buoyancy.freq"          "center.buoyancy"       
##  [3] "internal.energy"        "schmidt.stability"     
##  [5] "thermo.depth"           "water.density"         
##  [7] "epi.temperature"        "hypo.temperature"      
##  [9] "water.temperature"      "whole.lake.temperature"

heat_fluxes <- get_var(out_file, var_name = c("I_0"))

ggplot(data = heat_fluxes) + 
  geom_point(aes(x = DateTime, y = I_0))

epi_temp <- get_var(out_file, var_name = "temp")
epi_temp$DateTime <- as.POSIXct(substring(epi_temp$DateTime, 0, 20))

epi_temp_measured <- readxl::read_excel(
  file.path(
    paste(pathological::split_path(sim_folder)[[1]][c(-10)], 
        collapse = "/"), 
    "data", "ObservedTempData.xlsx"), 
  col_names = FALSE)[,c(2:4)]
names(epi_temp_measured) <- c("hours", "date", "temp_measured")
epi_temp_measured$DateTime <- as.POSIXct(paste0(epi_temp_measured$date, 
                                                " 00:00:00"))

epi_temp <- dplyr::left_join(epi_temp, epi_temp_measured)

## Joining, by = "DateTime"

lake_csv <- read.csv(file.path(sim_folder, "lake.csv"), 
                     stringsAsFactors = FALSE)
lake_csv$time <- as.POSIXct(lake_csv$time)
epi_temp <- dplyr::left_join(epi_temp, lake_csv, 
                             by = c("DateTime" = "time"))
lake_csv <- tidyr::gather(lake_csv, variable, value, -time)
lake_csv <- dplyr::filter(lake_csv, 
               variable %in% c("Tot.Inflow.Vol", "Tot.Outflow.Vol", 
                               "Rain", "Evaporation"))

# plot(epi_temp$DateTime, epi_temp$temp.elv_0)
# plot(epi_temp$DateTime, epi_temp$temp.elv_13.2350650316971)
# plot(epi_temp$DateTime, epi_temp$temp.elv_37.4993509231419)
# plot(epi_temp$DateTime, epi_temp$temp_measured)

epi_temp <- select(epi_temp, -date, -hours)
names(epi_temp) <- substring(names(epi_temp), 0, 12)
epi_temp <- remove_empty_cols(epi_temp)

cor_mat <- cor(
  select(epi_temp, -DateTime)[complete.cases(epi_temp),]
  )

## Warning in cor(select(epi_temp, -DateTime)[complete.cases(epi_temp), ]):
## the standard deviation is zero

cor_mat[lower.tri(cor_mat)] <- NA

epi_temp <- dplyr::select(epi_temp, 
                          c("DateTime", "temp.elv_37.", 
                            "temp_measure", "Surface.Temp"))

res <- tidyr::gather(epi_temp, variable, value, -DateTime)

# epi_temp$measured <- "modeled"
# epi_temp[epi_temp$variable == "temp_measured", "measured"] <- "measured"

ggplot(data = res) + 
  geom_point(aes(x = DateTime, y = value, color = variable))

## Warning: Removed 3583 rows containing missing values (geom_point).

Yes FABM

sim_folder <- "../release/warmlake/fabm"
# sim_folder <- "release/warmlake/fabm"
glm_nml_path <- file.path(sim_folder, "glm2.nml")
glm_nml <- read_nml(glm_nml_path)
fabm_nml_path <- file.path(sim_folder, "fabm.nml")
fabm_nml <- read_nml(fabm_nml_path)

names(fabm_nml)

##  [1] "fabm_nml"           "gotm_npzd"          "aed_sedflux"       
##  [4] "aed_sed_constant"   "aed_oxygen"         "aed_carbon"        
##  [7] "aed_silica"         "aed_nitrogen"       "aed_phosphorus"    
## [10] "aed_organic_matter" "aed_phytoplankton"  "aed_zooplankton"   
## [13] "aed_tracer"

get_nml_value(glm_nml, "csv_point_vars")

## [1] "temp,salt,aed_oxygen_oxy,gotm_npzd_nut"

glm_nml <- set_nml(glm_nml, "csv_point_vars", 
            arg_val = "temp,salt,aed_oxygen_oxy,gotm_npzd_nut")
write_nml(glm_nml, file.path(sim_folder, "glm2.nml"))

run_glm(sim_folder)

lake_csv_fabm <- read.csv(file.path(sim_folder, "lake.csv"), 
                     stringsAsFactors = FALSE)
lake_csv_fabm$time <- as.POSIXct(lake_csv_fabm$time)
lake_csv_fabm <- tidyr::gather(lake_csv_fabm, variable, value, -time)
lake_csv_fabm <- dplyr::filter(lake_csv_fabm, 
               variable %in% c("Tot.Inflow.Vol", "Tot.Outflow.Vol", 
                               "Rain", "Evaporation"))

ggplot(data = lake_csv_fabm) + 
  geom_line(aes(x = time, y = value, color = variable))

out_file_fabm <- file.path(sim_folder, "output.nc")
(s_vars_fabm <- sim_vars(out_file))

var_names <- paste0("SDF_Fsed_", 
                    c("amm", "ch4", "dic", 
                      "doc", "don", "dop",                                                            "feii", "frp"))

for(i in seq_len(length(var_names))){
  p_sed_rate <- get_var(out_file, var_name = var_names[i])
  plot(p_sed_rate)  
}

plot(get_var(out_file, "PHS_frp")[,c(1, 20)])

get_var(out_file, "OGM_doc")