#' Clam bioenergetic population model (alternative version) preprocessor
#' @param userpath the path where folder containing model inputs and outputs is located
#' @param forcings a list containing model forcings
#' @return a list containing the time series in the odd positions and realted forcings in the even positions. Forcings returned are: Water temperature [Celsius degrees], Chlorophyll a concentration [mgChl-a/m^3]
#' @import matrixStats plotrix rstudioapi
#' @import grDevices graphics utils stats


  rm(list=ls())       # Clean workspace

  cat("Data preprocessing")

  # Extracts forcings values from the list

# Read forcings and parameters from .csv files
Param_matrix=read.csv(paste0(userpath,"/ClamF_population/Inputs/Parameters//Parameters.csv"),sep=",")                      # Reading the matrix containing parameters and their description

  # Extract parameters and forcing values from parameters matrix and convert to type 'double' the vector contents
Param=as.matrix(Param_matrix[1:25,3])     # Vector containing all parameters
Dates=Param_matrix[27:28,3]                          # Vector containing the starting and ending date of teh simulation
IC=as.double(as.matrix(Param_matrix[26,3]))          # Initial weight condition
CS=as.double(as.matrix(Param_matrix[29,3]))                # Commercial size

# Prepare data for ODE solution
t0=min(as.numeric(as.Date(timeT[1], "%d/%m/%Y")), as.numeric(as.Date(timeChl[1], "%d/%m/%Y")),  as.numeric(as.Date(Dates[1], "%d/%m/%Y"))) # starting minimum starting date for forcings and observations
timestep=1                                           # Time step for integration [day]
ti=as.numeric(as.Date(Dates[1], "%d/%m/%Y"))-t0      # Start of integration [day]
tf=as.numeric(as.Date(Dates[2], "%d/%m/%Y"))-t0      # End of integration [day]
Ww=as.vector(matrix(0,nrow=ti))                      # Initialize vector weight
Ww[ti]=IC                                            # Wet weight initial value [g]
times<-cbind(ti, tf, timestep)                       # Vector with integration data

# Initial conditions
a=Param[8]               # [cm] Wet weight - length conversion coefficient
k=Param[6]               # [-] Wet weight - length conversion exponent
b=Param[9]               # [gWW] Dry weight - wet weight conversion coefficient
p=Param[5]               # [-] Dry weight - wet weight conversion exponent

Wd=as.vector(matrix(0,nrow=ti))   # Initialize vector dry weight
Wd[ti]=b*Ww[ti]^p                  # Dry weight initial value [g]
L=as.vector(matrix(0,nrow=ti))    # Initialize vector length
L[ti]=(Ww[ti]/a)^k                # Length of the mussel initial value [cm]

# Read files with population parameters and management strategies
Pop_matrix=read.csv(paste0(userpath,"/ClamF_population/Inputs/Parameters//Population.csv"),sep=",")   # Reading the matrix containing population parameters and their description
Management=read.csv(paste0(userpath,"/ClamF_population/Inputs/Population_management//Management.csv"),sep=",")   # Reading the matrix containing seeding and harvesting management

# Extract population parameters
meanWw=as.double(as.matrix(Pop_matrix[1,3]))    # [g] Wet weight average
deltaWw=as.double(as.matrix(Pop_matrix[2,3]))   # [g] Wet weight standard deviation
Wwlb=as.double(as.matrix(Pop_matrix[3,3]))      # [g] Wet weight lower bound
meanGdmax=as.double(as.matrix(Pop_matrix[4,3])) # [l/d gDW] Maximum growth rate on a dry weight average
deltaGdmax=as.double(as.matrix(Pop_matrix[5,3]))# [l/d gDW] Maximum growth rate on a dry weight standard deviation
Nseed=as.double(as.matrix(Pop_matrix[6,3]))     # [-] number of seeded individuals
mort=as.double(as.matrix(Pop_matrix[7,3]))      # [1/d] natural mortality rate
nruns=as.double(as.matrix(Pop_matrix[8,3]))     # [-] number of runs for population simulation

# Prepare management values
for (i in 1:length(Management[,1])) {
     manag[i,1]=as.numeric(as.Date(Management[i,1], "%d/%m/%Y"))-t0
  if ((Management[i,2])=="h") {
    } else {

# Population differential equation solution
N<-Pop_fun(Nseed, mort, manag, times)

# Print to screen inserted parameters

cat(" \n")
cat('The model will be executed with the following parameters:\n');
cat(" \n")
for (i in 1:25){
cat(paste0(toString(Param_matrix[i,2]), ": ", toString(Param_matrix[i,3]), " " ,toString(Param_matrix[i,4])),"\n")

cat(" \n")
cat("Integration is performed between ", toString(Dates[1]), " and ", toString(Dates[2]),"\n")
cat(" \n")
cat('Commercial size is ', toString(CS)," mm")
cat(" \n")

# Print to screen population characteristics

cat(" \n")
cat('The population is simulated by assuming that initial wet weight and maximum growth rate on a dry weight basis are normally distributed:\n');
cat(" \n")
for (i in 1:5){
  cat(paste0(toString(Pop_matrix[i,2]), ": ", toString(Pop_matrix[i,3]), " " ,toString(Pop_matrix[i,4])),"\n")

cat(" \n")
cat("The population is initially composed by ", toString(Pop_matrix[6,3]), " Individuals\n")
cat(" \n")
cat("The mortality rate is:", toString(Pop_matrix[7,3]),'1/d\n' )

# Print to screen management actions
cat(" \n")
cat('The population is managed according with following list (h:harvesting s:seeding):\n');
cat(" \n")
for (i in 1:length(Management[,1])){
  cat(paste0(toString(Management[i,1])," ", toString(Management[i,2]), " " ,toString(Management[i,3])),"individuals\n")

cat(" \n")
cat("The individual model will be executed ", toString(nruns), " times in order to simulate a population\n")
cat(" \n")

# Plot to file inserted forcing functions

cat(" \n")
cat("Forcings are represented in graphs available at the following folder:\n")

# plot Temperature forcing
days <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), by = "days", length = tf-ti+1)
plot(days, Tintsave, ylab="Water temperature (Celsius degrees)", xlab="",xaxt = "n",type="l",cex.lab=1.4)
labDates <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), tail(days, 1), by = "months")
axis.Date(side = 1, days, at = labDates, format = "%d %b %y", las = 2)

# Plot Chlorophyll a forcing
days <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), by = "days", length = tf-ti+1)
plot(days, Chlaintsave, ylab="Chlorophyll_a (mg/m^3)", xlab="", xaxt = "n",type="l",cex.lab=1.4)
labDates <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), tail(days, 1), by = "months")
axis.Date(side = 1, days, at = labDates, format = "%d %b %y", las = 2)

# plot population dynamics
days <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), by = "days", length = tf-ti)
plot(days, Nsave, ylab="Individuals", xlab="", xaxt = "n",type="l",cex.lab=1.4)
labDates <- seq(as.Date(Dates[1], format = "%d/%m/%Y"), tail(days, 1), by = "months")
axis.Date(side = 1, days, at = labDates, format = "%d %b %y", las = 2)

output=list(Param, times, Dates, IC, Tint, Chlint,N, CS)

Try the RAC package in your browser

Any scripts or data that you put into this service are public.

RAC documentation built on May 2, 2019, 3:26 a.m.