Update to v2

Author: jaseeverett

README.md

@@ -1,14 +1,55 @@
-# Model code for Zooplankton Model of Size Spectrum (ZooMSS v1)
+# Model code for Zooplankton Model of Size Spectrum (ZooMSS v2)
 
-This repo holds the code for Zooplankton Model of Size Spectrum (ZooMSS v1) as published in Heneghan et al. 2016 (pdf attached).   
 
-The size-spectrum model consists of a size-spectrum comprised of three communities: phytoplankton, zooplankton, and fish. The phytoplankton component covers the smallest size classes and is held constant as a background resource spectrum for zooplankton. Size-dependent processes of growth and mortality drive the zooplankton and fish components.  
 
-Heneghan RF, Everett JD, Blanchard JL and Richardson AJ (2016) Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs. *Front. Mar. Sci.* **3**:201. doi: 10.3389/fmars.2016.00201
-https://doi.org/10.3389/fmars.2016.00201  
+## Overview of ZooMSS
+
+The Zooplankton Model of Size Spectra (ZooMSS) is a functional size-spectrum model of the marine ecosystem (following Heneghan et al. 2016, ZooMSSv1) to resolve phytoplankton, nine zooplankton functional groups (heterotrophic flagellates and ciliates, omnivorous and carnivorous copepods, larvaceans, euphausiids, salps, chaetognaths and jellyfish) and three size-based fish groups. Zooplankton functional groups are resolved using their body-size ranges, size- based feeding characteristics and carbon content, and the zooplankton community emerges from the model across global environmental gradients, depending on the functional traits of the different groups. 
+
+We developed the Zooplankton Model of Size Spectra version 2 (ZooMSSv2) based on the prototype of Heneghan et al. (2016). ZooMSSv2 uses the functional size-spectrum framework (Blanchard et al., 2017) to resolve the body size ranges, size-based feeding characteristics and carbon content of nine zooplankton groups and three fish groups. The model is run on 5° grid squares across the global ocean. For each region, the model is forced with the long- term mean satellite sea surface temperature and chlorophyll a from MODIS-Aqua.
+
+ZooMSSv2 represents the marine ecosystem as three communities: phytoplankton, zooplankton and fish. The zooplankton community consists of nine of the most abundant zooplankton groups, and the fish community was made up of a small, medium and large group. Dynamics of the phytoplankton are not explicitly resolved in the model, rather the mean size structure of the phytoplankton community is estimated directly from satellite chlorophyll a observations (Brewin et al., 2010; Barnes et al., 2011; Hirata et al., 2011). Abundances of the zooplankton and fish communities are driven by size-dependent processes of growth and mortality, with the temporal dynamics of each functional group governed by separate second- order McKendrick-von Foerster equations.
+
+
+## How to run ZooMSS using `R`
+
+### Model Structure
+
+The model is split into 5 main files.
+
+1. Setup_RUN_NAME.R - Set up your own specific simulation. This is the only file you usually need to edit.
+
+2. fZooMSS_Model.R - The top-level model function which seuqentially runs the following 3 functions, and saves the output.
+
+    + fZooMSS_Params.R - Set all the parameter values for the model.
+  
+    + fZooMSS_Setup.R - Do all calculations which can be done in advance of the time loop.
+  
+    + fZooMSS_Run.R - Contains the time-dependent loop.
+
+### Input files
+
+ZooMSS requires two input files:   
+
+1. TestGroups.csv which contains all the taxa-specific parameter values, including size ranges.   
+
+2. And RDS file which contains all environmental data for each grid cell. This includes chlorophyll a (chlo), sea surface temperature (sst), bathymetry (bathy), phytoplankton slope and intercept, and latitude/longitude which is useful for plotting.   
+
+
+### Running
+To run the model, you simply need to provide the relevent information in Setup_RUN_NAME.R and run the script.
+
+## Publications
+
+1. Heneghan, R.F., Everett, J.D., Blanchard, J.L., Richardson, A.J., 2016. Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs. Front. Mar. Sci. 3, 1–15. https://doi.org/10.3389/fmars.2016.00201
+
+2. Heneghan, R.F., Everett, J.D., Sykes, P., Batten, S.D., Edwards, M., Takahashi, K., Suthers, I.M.,  Blanchard, J.L., Richardson, A.J., in review, A global size-spectrum model of the marine ecosystem that resolves zooplankton composition. Ecological Modelling
+
+
+## Getting Help
+
+If you are running the model come across a problem, raise an issue on GitHub: https://github.com/MathMarEcol/ZoopSizeSpectraModel/issues
+
+If you find some errors, or just want to improve the model, we'd love you to go ahead and make the changes and then submit a pull request for us to check and approve.
+
 
-## SCRIPTS 
-ZooMSS_Shell.R: You run the model from this script.    
-ZooMSS_Parameters.R: Parameters for the model.   
-ZooMSS_Stability: Code for stability analysis via Newton-Raphson method.    
-ZooMSS_PDE: The model code.    

TestGroups.csv

@@ -0,0 +1,13 @@
+Species,Type,FeedType,Prop,W0,Wmax,Wmat,SearchCoef,SearchExp,PPMRscale,PPMR,FeedWidth,GrossGEscale,Carbon,Repro,PlotColour
+Flagellates,Zooplankton,Heterotroph,0.1,-10.7,-6.8,-8.8,640,0.8,1.5,NA,0.36,2.5,0.15,0,seagreen3
+Ciliates,Zooplankton,Heterotroph,0.1,-9.3,-6.3,-8.3,640,0.8,0.04,NA,0.47,2.5,0.15,0,seagreen4
+Larvaceans,Zooplankton,FilterFeeder,0.1,-6.4,-3.2,-5.2,640,0.8,-3,NA,0.7,2.5,0.02,0,lightgoldenrod4
+OmniCopepods,Zooplankton,Omnivore,0.04,-7.5,-3.5,-5.5,640,0.8,-0.5,NA,0.57,2.5,0.12,0,cornflowerblue
+CarnCopepods,Zooplankton,Carnivore,0.06,-7.5,-2.5,-4.5,640,0.8,1.5,NA,0.36,2.5,0.12,0,darkred
+Euphausiids,Zooplankton,Omnivore,0.1,-4.2,0.2,-1.8,640,0.8,-2,NA,0.7,2.5,0.12,0,darkblue
+Chaetognaths,Zooplankton,Carnivore,0.1,-5.9,-0.9,-2.9,640,0.8,1,NA,0.46,2.5,0.04,0,firebrick1
+Salps,Zooplankton,FilterFeeder,0.01,-4.7,1.4,-0.6,640,0.8,-2.68,NA,0.7,2.5,0.02,0,lightgoldenrod3
+Jellyfish,Zooplankton,Carnivore,0.01,-3,2,0,640,0.8,0.73,NA,0.52,2.5,0.005,0,orange
+Fish_Small,Fish,Carnivore,NA,-3,2,0,640,0.8,NA,100,1.3,2.5,0.1,0,lightcoral
+Fish_Med,Fish,Carnivore,NA,-3,4,2,640,0.8,NA,100,1.3,2.5,0.1,0,coral3
+Fish_Large,Fish,Carnivore,NA,-3,7,4,640,0.8,NA,100,1.3,2.5,0.1,0,coral4

ZooMSS.Rproj

@@ -11,3 +11,7 @@ Encoding: UTF-8
 
 RnwWeave: Sweave
 LaTeX: pdfLaTeX
+
+StripTrailingWhitespace: Yes
+
+QuitChildProcessesOnExit: Yes

ZooMSS_Plot_GlobalMap.R

@@ -0,0 +1,66 @@
+library(tidyverse)
+library(rnaturalearth)
+library(rnaturalearthdata)
+library(sf)
+library(raster)
+
+source("fZooMSS_Xtras.R")
+
+w <- 10^(seq(from = -10.7, to =  7, 0.1)) # Size bins of whole model
+
+res <- read_rds("Output/res_20200429_Control.RDS")
+groups <- read_csv("TestGroups.csv")
+enviro <- read_rds("envirofull_20200317.RDS")
+
+# Calculate biomass of each species for each cell
+bioms <- fZooMSS_SpeciesBiomass(res,w)
+bioms <- as.data.frame(t(as.data.frame(bioms)))
+colnames(bioms) <- groups$Species # Add names
+rownames(bioms) <- c() # Remove row names
+
+bioms <- bind_cols(bioms, enviro) # Bind enviro data
+
+# Which variable to plot
+var = "Fish_Small"
+
+# What to plot
+df <- bioms %>%
+  dplyr::select("Lon", "Lat", all_of(var))
+
+# Now do the plotting
+world <- ne_countries(scale = "medium", returnclass = "sf")
+world_sf <- st_transform(world, crs = 54009) # Convert to Mollweide
+
+df_raster <- rasterFromXYZ(df)  #Convert first two columns as lon-lat and third as value
+df_poly <- rasterToPolygons(df_raster, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=FALSE) # Convert to polygon which is better for plotting
+df_sf <- st_as_sf(df_poly, crs = 4326) # Convert to sf
+
+st_crs(df_sf) <- 4326 # Set base projection
+df_sf_moll <- st_transform(df_sf, crs = 54009) # Alter the CRS to mollweide
+
+gg_world <- ggplot() +
+  geom_sf(data = df_sf_moll, aes_string(fill = var), color = NA) +
+  geom_sf(data = world_sf, size = 0.05, fill = "grey50") +
+  scale_fill_gradient(low = "white",
+                       high = "lightcoral",
+                       space = "Lab",
+                       na.value = "grey50",
+                       aesthetics = "fill",
+                       oob = scales::squish,
+                       guide = guide_colourbar(title = paste0(var, " Biomass"),
+                                               title.position = "top",
+                                               title.hjust = 0.5)) +
+  theme(plot.background = element_rect(fill = NA),
+        plot.margin = unit(c(0,0,0,0), "mm"),
+        panel.grid.major = element_line(colour = "grey50", size = 0.2),
+        panel.background = element_blank(),
+        panel.border = element_blank(),
+        legend.position = "bottom",
+        legend.key.width = unit(20,"mm"),
+        legend.key.height = unit(4,"mm"))
+
+
+graphics.off()
+x11(width = 8, height = 6)
+gg_world
+ggsave(paste0("Figures/Biomass_",var,".png"), dpi = 300)

ZooMSS_Plot_RunTime.R

@@ -0,0 +1,40 @@
+library(tidyverse)
+library(patchwork)
+
+# Load all files
+fi <- list.files("RawOutput/", full.names = TRUE)
+
+# Create empty df to store data
+df <- tibble("No" = rep(NA, length(fi)),
+             "SST" = rep(NA, length(fi)),
+             "Chl" = rep(NA, length(fi)),
+             "RunTime" = rep(NA, length(fi)))
+
+for (i in 1:length(fi)){
+
+  dat <- read_rds(fi[i])
+
+  df$No[i] <- i
+  df$SST[i] <- dat$model$param$sst
+  df$Chl[i] <- dat$model$param$chlo
+  df$RunTime[i] <- as.numeric(dat$model$model_runtime[3])/60
+  rm(dat)
+
+}
+# Plot
+gg1 <- ggplot(data = df, aes(x = log10(Chl), y = RunTime)) +
+  geom_point() +
+  geom_smooth(method = "loess") +
+  theme_bw() +
+  ylab("Model Run Time (mins)")
+
+gg2 <- ggplot(data = df, aes(x = SST, y = RunTime)) +
+  geom_point() +
+  geom_smooth(method = "loess") +
+  theme_bw() +
+  ylab("Model Run Time (mins)")
+
+x11(width = 12, height = 6)
+gg1 / gg2
+ggsave("Figures/Model_RunTime.png", dpi = 300)
+

ZooMSS_Plot_SpeciesAccum.R

@@ -0,0 +1,75 @@
+
+## Now lets have a look at how the output changes in 50 year increments
+
+library(tidyverse)
+library(patchwork)
+fi <- list.files("RawOutput", full.names = TRUE)
+
+int <- 50
+tm <- seq(int,1000,int)
+
+# Create empty df to store data
+out <- tibble("SST" = numeric(),
+              "Chl" = numeric(),
+              "Run" = numeric(),
+              "Time" = numeric(),
+              "Abund" = numeric(),
+              "Growth" = numeric(),
+              "Mort" = numeric())
+
+for (i in 1:length(fi)){
+  dat <- read_rds(fi[i])
+  for (j in 1:length(tm)){
+
+    sp <- dat$model$param$Groups$Species
+
+    df <- tibble("SST" = rep(NA, length(sp)),
+                 "Chl" = rep(NA, length(sp)),
+                 "Run" = rep(NA, length(sp)),
+                 "Time" = rep(NA, length(sp)),
+                 "Abund" = rep(NA, length(sp)),
+                 "Growth" = rep(NA, length(sp)),
+                 "Mort" = rep(NA, length(sp)))
+
+    for (k in 1:length(sp)){
+
+      df$Run[k] <- i
+      df$Time[k] <- tm[j]
+      df$SST[k] <- dat$model$param$sst
+      df$Chl[k] <- dat$model$param$chlo
+      df$Species[k] <- sp[k]
+      df$Abund[k] <- sum(colMeans(dat$model$N[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+      df$Growth[k] <- sum(colMeans(dat$model$gg[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+      df$Mort[k] <- sum(colMeans(dat$model$Z[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+
+
+      # df$Abund[k] <- sum(colMeans(dat$model$N[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+      # df$Growth[k] <- sum(colMeans(dat$model$gg[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+      # df$Mort[k] <- sum(colMeans(dat$model$Z[(ceiling(0.5*tm[j]):tm[j]),k,], dims = 1))
+    }
+
+    out <- bind_rows(out, df)
+    rm(df)
+
+  }
+}
+
+
+x11(width = 12, height = 6)
+out %>%
+  # filter(Species == "OmniCopepods") %>%
+  mutate(Run = as.factor(Run),
+         Species = factor(Species, levels=sp),
+         ChlSST = factor(paste0(Chl, "Chl-", SST, "SST"))) %>%
+  ggplot(aes(x = Time, y = log10(Abund))) +
+  # geom_point(size  = 0.5) +
+  geom_line(aes(colour = Run), size = 0.1) +
+  facet_wrap(~ Species, ncol = 4, scales = "free") +
+  theme_bw() +
+  theme(legend.position = "none") +
+  scale_x_continuous(expand = c(0, 0)) +
+  scale_y_continuous(expand = c(0, 0))
+
+ggsave("Figures/SpeciesAccumulation.png", dpi = 200)
+
+

ZooMSS_RUN_NAME.pbs

@@ -0,0 +1,17 @@
+#!/bin/bash
+ 
+#PBS -l nodes=1:ppn=1,mem=1gb
+#PBS -l walltime=11:59:59
+#PBS -M [email protected]
+#PBS -m ae
+#PBS -o JobOutput/
+#PBS -e JobOutput/
+#PBS -J 1-1638
+#PBS -j oe
+
+module purge
+module add R/3.6.1
+
+cd $PBS_O_WORKDIR
+
+Rscript setup_RUNS_NAME.R

ZooMSS_User_Manual.Rmd

@@ -0,0 +1,95 @@
+---
+title: "ZooMSS User Manual"
+author: "The UQ ZooMSS Team"
+date: "18/03/2020"
+output: pdf_document
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+
+Note: This is still very much a work in progress.
+
+## Overview of ZooMSS
+
+The Zooplankton Model of Size Spectra (ZooMSS) is a functional size-spectrum model of the marine ecosystem (following Heneghan et al. 2016, ZooMSSv1) to resolve phytoplankton, nine zooplankton functional groups (heterotrophic flagellates and ciliates, omnivorous and carnivorous copepods, larvaceans, euphausiids, salps, chaetognaths and jellyfish) and three size-based fish groups. Zooplankton functional groups are resolved using their body-size ranges, size- based feeding characteristics and carbon content, and the zooplankton community emerges from the model across global environmental gradients, depending on the functional traits of the different groups. 
+
+We developed the Zooplankton Model of Size Spectra version 2 (ZooMSSv2) based on the prototype of Heneghan et al. (2016). ZooMSSv2 uses the functional size-spectrum framework (Blanchard et al., 2017) to resolve the body size ranges, size-based feeding characteristics and carbon content of nine zooplankton groups and three fish groups. The model is run on 5° grid squares across the global ocean. For each region, the model is forced with the long- term mean satellite sea surface temperature and chlorophyll a from MODIS-Aqua.
+
+ZooMSSv2 represents the marine ecosystem as three communities: phytoplankton, zooplankton and fish. The zooplankton community consists of nine of the most abundant zooplankton groups, and the fish community was made up of a small, medium and large group. Dynamics of the phytoplankton are not explicitly resolved in the model, rather the mean size structure of the phytoplankton community is estimated directly from satellite chlorophyll a observations (Brewin et al., 2010; Barnes et al., 2011; Hirata et al., 2011). Abundances of the zooplankton and fish communities are driven by size-dependent processes of growth and mortality, with the temporal dynamics of each functional group governed by separate second- order McKendrick-von Foerster equations.
+
+
+## How to run ZooMSS using `R`
+
+### Model Structure
+
+The model is split into 5 main files.
+
+1. Setup_RUN_NAME.R - Set up your own specific simulation. This is the only file you usually need to edit.
+
+2. fZooMSS_Model.R - The top-level model function which seuqentially runs the following 3 functions, and saves the output.
+
+    + fZooMSS_Params.R - Set all the parameter values for the model.
+  
+    + fZooMSS_Setup.R - Do all calculations which can be done in advance of the time loop.
+  
+    + fZooMSS_Run.R - Contains the time-dependent loop.
+
+### Input files
+
+ZooMSS requires two input files:   
+
+1. TestGroups.csv which contains all the taxa-specific parameter values, including size ranges.   
+
+2. And RDS file which contains all environmental data for each grid cell. This includes chlorophyll a (chlo), sea surface temperature (sst), bathymetry (bathy), phytoplankton slope and intercept, and latitude/longitude which is useful for plotting.   
+
+
+### Parameter Values
+
+Insert table from Heneghan et al (in review)
+
+### Output
+
+TBC
+
+### Running
+To run the model, you simply need to provide the relevent information in Setup_RUN_NAME.R and run the script.
+
+### Plotting
+
+TBC
+
+## Useful reading
+
+Datta Paper
+The other thing to look at would be some of the numerical methods we use 
+  - Simpson’s method for the growth, mortality and diffusion integrals and 
+  - a semi-implicit upward finite difference scheme to solve the PDE at each time step. 
+  For the PDE numerical method, check out Appendix 1 of my thesis (attached)
+  Info on the difference scheme in Press’ 2007 book: Numerical Recipes, the Art of Scientific Computing (Available at UQ Library)
+
+
+On 1 Apr 2020, at 06:50, Richardson, Anthony (O&A, St. Lucia) <[email protected]> wrote:
+
+<Datta et al. - 2010 - A Jump-Growth M odel for Predator–Prey Dynamics De.pdf>
+
+## Publications
+
+1. Heneghan, R.F., Everett, J.D., Blanchard, J.L., Richardson, A.J., 2016. Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs. Front. Mar. Sci. 3, 1–15. https://doi.org/10.3389/fmars.2016.00201
+
+2. Heneghan, R.F., Everett, J.D., Sykes, P., Batten, S.D., Edwards, M., Takahashi, K., Suthers, I.M.,  Blanchard, J.L., Richardson, A.J., in review, A global size-spectrum model of the marine ecosystem that resolves zooplankton composition. Ecological Modelling
+
+## References
+
+TBC
+
+## Getting Help
+
+If you are running the model come across a problem, raise an issue on GitHub: https://github.com/MathMarEcol/ZooMSS/issues
+
+If you find some errors, or just want to improve the model, we'd love you to go ahead and make the changes and then submit a pull request for us to check and approve.
+
+Same for this manual. If there is information lacking, errors etc, please fix it rather than just reporting it. This should be a team effort.
+
+