tests/testthat/test-fish-fluxes.R

In StrathE2E2: End-to-End Marine Food Web Model

test_that("check fish feeding flux calculations", {

#-----------------------------------------------------------------------------------------------------------------
#
# This test uses a 'Testbed' version of the model in which all of the external drivers are set to constant values.
# We expect the model outputs to converge to a steady state under these conditions. Then we take the masses of the
# food web variables from the steady state model output, and manually derive the uptake fluxes between them.
# These manually derived fluxes are then compared with the fluxes generated in the 'flux_matrix' outputs from the model.
# So the test is checking that the flux calculations in the C-code of the package are correctly evaluated.
#
# The Testbed model is stored in /StrathE2E2/tests/testdata/models/
#
# The model has been run for >1000 years to ensure a steady state, and the end-state saved as initial conditions.
# So the model should be at steady state from time 0. In fact, the input csv file containing the initial conditions
# does not hold sufficient precision (number of decimal places) to ensure perfect steady state from the onset of a run.
# This is manifest as some tiny oscillations in the first year or so of a run before everything settles to its steady state.
# Hence we run the Testbed model for 10 years and then base the test on the final year of the run.
# 
# For convenience, the Testbed model uses the fitted parameters from the North Sea implementation of StrathE2E2.
# 
# Two important points to note:
# 1) The temperature driving data for all three spatial compartmets of the Testbed
# model are set to 10C, wich matches the Q10 reference temperature in the fixed parameter inputs. This means that when manually
# calculating the uptake fluxes there is no need to compute any temperature corrections for the maximum uptake rate
# parameters. However, the C-code in the model IS doing the Q10 calculation, so this is an important aspect of the testing.
# As a consequence, it is very important not to chnage the temperature driving values in the Testbed model - they must be
# equal to the Q10 reference temperature value.
# 2) The sediment habitats are all set to the same area (0.125). This affects the calculation of uptake of corpses by demersal fish.
# By setting the habitats to equal area, we can simply sum the mass of corpses across each zone and use that as input to the
# manual calculation of uptake by demersal fish. The C code performs this calculation separately for each habitat and then adds the
# results. So this is an additional test of the model code. But it means that the habitat areas should not be edited in the
# input csv files.
#
# This test applies the above approach to uptake by planktivorous, demersal and migratory fish classes.
#
#-----------------------------------------------------------------------------------------------------------------
#
# The following non-exported functions in the package are needed for the test:
# source("internal.R")
# source("read_fitted_parameters.R")
# source("read_fixed_parameters.R")
# source("read_physics_drivers.R")
# source("read_physical_parameters.R")
#
#-----------------------------------------------------------------------------------------------------------------

#Some additional functions required for assembling the tests

#######################
#Function to pull the final year annual average mass of a variable out of the results
pull_final_year_mean<-function(results,vtoget){
vgot<-(elt(results$output, vtoget))
vmean<-mean(vgot[(length(vgot)-360):(length(vgot))])
vmean
}
######################
#Function to pull the final year annual flux of a variable out of the results
pull_final_year_flux<-function(results,vtoget){
vgot<-(elt(results$output, vtoget))
vflux<- vgot[(length(vgot))]   -  vgot[(length(vgot)-360)]
vflux
}
#######################
#Function to pull flux data out of the flow matrix generated by the model
extract_model_flux<-function(matrixname,preyname,predname){
predi<-which(colnames(matrixname)==predname)
preyi<-which(colnames(matrixname)==preyname)
model_flux<-(matrixname[preyi,predi]/360)
model_flux
}
#######################
#Heterotrophic uptake function
f1<-function(a,b,k1,k2)
#          prey,pred,umax,hs
{
  return ((b*k1*a)/(a+k2))
}
#######################




#-----------------------------------------------------------------------------------------------------------------

#Run the Testbed model and pull out data from the results object

#Run the testbed model - assuming R home dir is "/GitLab/StrathE2E2/tests"
model<-e2e_read(model.name="Testbed",
                  model.variant="Const",
                  model.ident="base",
                  models.path="../testdata/models")
results <- e2e_run(model,nyears=10,csv.output=FALSE)
#plot_full_length_timeseries(model,results)
#extract_new_initial_cond_from_end_of_run(model, results)

#Pull out the inshore and offshore final year data for phytoplankton and zoopkankton groups

#names(results$output[1:78]) # Lists the names of state varibales in the header

corpse_s0<-pull_final_year_mean(results,"corpse_s0")
corpse_s1<-pull_final_year_mean(results,"corpse_s1")
corpse_s2<-pull_final_year_mean(results,"corpse_s2")
corpse_s3<-pull_final_year_mean(results,"corpse_s3")
discard_i<-pull_final_year_mean(results,"discard_i")
benthslar_i<- pull_final_year_mean(results,"benthslar_i")
benthclar_i<- pull_final_year_mean(results,"benthclar_i")
fishplar_i<- pull_final_year_mean(results,"fishplar_i")
fishdlar_i<- pull_final_year_mean(results,"fishdlar_i")
omni_i<-pull_final_year_mean(results,"omni_i")
carn_i<-pull_final_year_mean(results,"carn_i")
benths_i<-pull_final_year_mean(results,"benths_i")
benthc_i<-pull_final_year_mean(results,"benthc_i")
fishp_i<-pull_final_year_mean(results,"fishp_i")
fishd_i<-pull_final_year_mean(results,"fishd_i")
fishm_i<-pull_final_year_mean(results,"fishm_i")

corpse_i <- corpse_s0 + corpse_s1 + corpse_s2 + corpse_s3

corpse_d0<-pull_final_year_mean(results,"corpse_d0")
corpse_d1<-pull_final_year_mean(results,"corpse_d1")
corpse_d2<-pull_final_year_mean(results,"corpse_d2")
corpse_d3<-pull_final_year_mean(results,"corpse_d3")
discard_o<-pull_final_year_mean(results,"discard_o")
benthslar_o<- pull_final_year_mean(results,"benthslar_o")
benthclar_o<- pull_final_year_mean(results,"benthclar_o")
fishplar_o<- pull_final_year_mean(results,"fishplar_o")
fishdlar_o<- pull_final_year_mean(results,"fishdlar_o")
omni_o<- pull_final_year_mean(results,"omni_o")
carn_o<- pull_final_year_mean(results,"carn_o")
benths_o<-pull_final_year_mean(results,"benths_o")
benthc_o<-pull_final_year_mean(results,"benthc_o")
fishp_o<-pull_final_year_mean(results,"fishp_o")
fishd_o<-pull_final_year_mean(results,"fishd_o")
fishm_o<-pull_final_year_mean(results,"fishm_o")

corpse_o <- corpse_d0 + corpse_d1 + corpse_d2 + corpse_d3

Pfish_recruit<-(pull_final_year_flux(results,"Pfish_recruit"))/360
Dfish_recruit<-(pull_final_year_flux(results,"Dfish_recruit"))/360

#Pull out the flow matrix from the results
flow_matrix<-elt(results$final.year.outputs,"flow_matrix_all_fluxes")
flow_matrix_ns<-elt(results$final.year.outputs,"flow_matrix_excl_spawn_recruit")

#colnames(flow_matrix)   #Lists the column names in the flow matrix
#rownames(flow_matrix)   #Lists the row names in the flow matrix

#-----------------------------------------------------------------------------------------------------------------




#Pull data out of the model object, ie data that have been assembled from the input csv files

#Pull out the model configuration to get areas and volumes

model.path<-model$setup$model.path
	read.model.setup(model.path)			# Models/Model/Variant/MODEL_SETUP.csv
	physical.parameters	<- read_physical_parameters(model.path)
	so_depth<-elt(physical.parameters,"so_depth")
	si_depth<-elt(physical.parameters,"si_depth")
	d_depth<-elt(physical.parameters,"d_depth")
        area_s0<-elt(physical.parameters,"x_area_s0")
        area_s1<-elt(physical.parameters,"x_area_s1")
        area_s2<-elt(physical.parameters,"x_area_s2")
        area_s3<-elt(physical.parameters,"x_area_s3")
        area_d0<-elt(physical.parameters,"x_area_d0")
        area_d1<-elt(physical.parameters,"x_area_d1")
        area_d2<-elt(physical.parameters,"x_area_d2")
        area_d3<-elt(physical.parameters,"x_area_d3")
shallowprop<-(area_s0+area_s1+area_s2+area_s3)
volume_so<-so_depth*(1-shallowprop)
volume_d<-d_depth*(1-shallowprop)
volume_si<-si_depth*(shallowprop)



# Pull out the uptake function parameters needed to manually calculate the fluxes

model.path<-model$setup$model.path
	read.model.setup(model.path)			# Models/Model/Variant/MODEL_SETUP.csv
	fitted.parameters	<- read_fitted_parameters(model.path)
	fixed.parameters	<- read_fixed_parameters(model.path)

CZ_inedible <- elt(fixed.parameters,"CZ_inedible")

#Calculate edible component of carn zoo
CZ_inedible_i <- CZ_inedible * shallowprop
CZ_inedible_o <- CZ_inedible * (1-shallowprop)
 if(carn_o>0) CZ_edible_o  <- max(0,carn_o-(CZ_inedible_o/carn_o))   else   CZ_edible_o  <- 0
 if(carn_i>0) CZ_edible_i  <- max(0,carn_i-(CZ_inedible_i/carn_i))   else   CZ_edible_i  <- 0


u_fishp<-elt(fitted.parameters,"u_fishp")
h_fishp<-elt(fitted.parameters,"h_fishp")

PREF_omni_fishp<-elt(fitted.parameters,"PREF_omni_fishp")
PREF_carn_fishp<-elt(fitted.parameters,"PREF_carn_fishp")
PREF_benthslar_fishp<-elt(fitted.parameters,"PREF_benthslar_fishp")
PREF_benthclar_fishp<-elt(fitted.parameters,"PREF_benthclar_fishp")
PREF_fishplar_fishp<-elt(fitted.parameters,"PREF_fishplar_fishp")
PREF_fishdlar_fishp<-elt(fitted.parameters,"PREF_fishdlar_fishp")

#....

u_fishd<-elt(fitted.parameters,"u_fishd")
h_fishd<-elt(fitted.parameters,"h_fishd")

PREF_corpse_fishd<-elt(fitted.parameters,"PREF_corp_fishd")
PREF_discard_fishd<-elt(fitted.parameters,"PREF_disc_fishd")
PREF_carn_fishd<-elt(fitted.parameters,"PREF_carn_fishd")
PREF_benths_fishd<-elt(fitted.parameters,"PREF_benths_fishd")
PREF_benthc_fishd<-elt(fitted.parameters,"PREF_benthc_fishd")
PREF_fishplar_fishd<-elt(fitted.parameters,"PREF_fishplar_fishd")
PREF_fishdlar_fishd<-elt(fitted.parameters,"PREF_fishdlar_fishd")
PREF_fishp_fishd<-elt(fitted.parameters,"PREF_fishp_fishd")
PREF_fishd_fishd<-elt(fitted.parameters,"PREF_fishd_fishd")
PREF_fishm_fishd<-elt(fitted.parameters,"PREF_fishm_fishd")

#....

u_fishm<-elt(fitted.parameters,"u_fishm")
h_fishm<-elt(fitted.parameters,"h_fishm")

PREF_omni_fishm<-elt(fitted.parameters,"PREF_omni_fishm")
PREF_carn_fishm<-elt(fitted.parameters,"PREF_carn_fishm")
PREF_benthslar_fishm<-elt(fitted.parameters,"PREF_benthslar_fishm")
PREF_benthclar_fishm<-elt(fitted.parameters,"PREF_benthclar_fishm")
PREF_fishplar_fishm<-elt(fitted.parameters,"PREF_fishplar_fishm")
PREF_fishdlar_fishm<-elt(fitted.parameters,"PREF_fishdlar_fishm")


#-----------------------------------------------------------------------------------------------------------------

#Now calculate the fluxes


#............................

calc_flux_omni_fishp_o<- f1(omni_o,(fishp_o),u_fishp*PREF_omni_fishp,h_fishp*(volume_so+volume_d))
calc_flux_omni_fishp_i<- f1(omni_i,(fishp_i),u_fishp*PREF_omni_fishp,h_fishp*(volume_si))

calc_flux_omni_fishp_w <- calc_flux_omni_fishp_o + calc_flux_omni_fishp_i

model_flux_omni_fishp<-extract_model_flux(flow_matrix,"omnivzoo","pfish")

#calc_flux_omni_fishp_w
#model_flux_omni_fishp

#............................

calc_flux_carn_fishp_o<- f1(CZ_edible_o,(fishp_o),u_fishp*PREF_carn_fishp,h_fishp*(volume_so+volume_d))
calc_flux_carn_fishp_i<- f1(CZ_edible_i,(fishp_i),u_fishp*PREF_carn_fishp,h_fishp*(volume_si))

calc_flux_carn_fishp_w <- calc_flux_carn_fishp_o + calc_flux_carn_fishp_i

model_flux_carn_fishp<-extract_model_flux(flow_matrix,"carnzoo","pfish")

#calc_flux_carn_fishp_w
#model_flux_carn_fishp

#............................

calc_flux_benthslar_fishp_o<- f1(benthslar_o,(fishp_o),u_fishp*PREF_benthslar_fishp,h_fishp*(volume_so+volume_d))
calc_flux_benthslar_fishp_i<- f1(benthslar_i,(fishp_i),u_fishp*PREF_benthslar_fishp,h_fishp*(volume_si))

calc_flux_benthslar_fishp_w <- calc_flux_benthslar_fishp_o + calc_flux_benthslar_fishp_i

model_flux_benthslar_fishp<-extract_model_flux(flow_matrix,"benthslar","pfish")

#calc_flux_benthslar_fishp_w
#model_flux_benthslar_fishp

#............................

calc_flux_benthclar_fishp_o<- f1(benthclar_o,(fishp_o),u_fishp*PREF_benthclar_fishp,h_fishp*(volume_so+volume_d))
calc_flux_benthclar_fishp_i<- f1(benthclar_i,(fishp_i),u_fishp*PREF_benthclar_fishp,h_fishp*(volume_si))

calc_flux_benthclar_fishp_w <- calc_flux_benthclar_fishp_o + calc_flux_benthclar_fishp_i

model_flux_benthclar_fishp<-extract_model_flux(flow_matrix,"benthclar","pfish")

#calc_flux_benthclar_fishp_w
#model_flux_benthclar_fishp

#............................

calc_flux_fishplar_fishp_o<- f1(fishplar_o,(fishp_o),u_fishp*PREF_fishplar_fishp,h_fishp*(volume_so+volume_d))
calc_flux_fishplar_fishp_i<- f1(fishplar_i,(fishp_i),u_fishp*PREF_fishplar_fishp,h_fishp*(volume_si))

calc_flux_fishplar_fishp_w <- calc_flux_fishplar_fishp_o + calc_flux_fishplar_fishp_i

model_flux_fishplar_fishp<-(extract_model_flux(flow_matrix,"pfishlar","pfish")) - Pfish_recruit  # Flow matrix includes the recruitment flux as well

#calc_flux_fishplar_fishp_w
#model_flux_fishplar_fishp

#............................

calc_flux_fishdlar_fishp_o<- f1(fishdlar_o,(fishp_o),u_fishp*PREF_fishdlar_fishp,h_fishp*(volume_so+volume_d))
calc_flux_fishdlar_fishp_i<- f1(fishdlar_i,(fishp_i),u_fishp*PREF_fishdlar_fishp,h_fishp*(volume_si))

calc_flux_fishdlar_fishp_w <- calc_flux_fishdlar_fishp_o + calc_flux_fishdlar_fishp_i

model_flux_fishdlar_fishp<-extract_model_flux(flow_matrix,"dfishlar","pfish")

#calc_flux_fishdlar_fishp_w
#model_flux_fishdlar_fishp

#............................
#............................

calc_flux_omni_fishm_o<- f1(omni_o,(fishm_o),u_fishm*PREF_omni_fishm,h_fishm*(volume_so+volume_d))
calc_flux_omni_fishm_i<- f1(omni_i,(fishm_i),u_fishm*PREF_omni_fishm,h_fishm*(volume_si))

calc_flux_omni_fishm_w <- calc_flux_omni_fishm_o + calc_flux_omni_fishm_i

model_flux_omni_fishm<-extract_model_flux(flow_matrix,"omnivzoo","mfish")

#calc_flux_omni_fishm_w
#model_flux_omni_fishm

#............................

calc_flux_carn_fishm_o<- f1(CZ_edible_o,(fishm_o),u_fishm*PREF_carn_fishm,h_fishm*(volume_so+volume_d))
calc_flux_carn_fishm_i<- f1(CZ_edible_i,(fishm_i),u_fishm*PREF_carn_fishm,h_fishm*(volume_si))

calc_flux_carn_fishm_w <- calc_flux_carn_fishm_o + calc_flux_carn_fishm_i

model_flux_carn_fishm<-extract_model_flux(flow_matrix,"carnzoo","mfish")

#calc_flux_carn_fishm_w
#model_flux_carn_fishm

#............................

calc_flux_benthslar_fishm_o<- f1(benthslar_o,(fishm_o),u_fishm*PREF_benthslar_fishm,h_fishm*(volume_so+volume_d))
calc_flux_benthslar_fishm_i<- f1(benthslar_i,(fishm_i),u_fishm*PREF_benthslar_fishm,h_fishm*(volume_si))

calc_flux_benthslar_fishm_w <- calc_flux_benthslar_fishm_o + calc_flux_benthslar_fishm_i

model_flux_benthslar_fishm<-extract_model_flux(flow_matrix,"benthslar","mfish")

#calc_flux_benthslar_fishm_w
#model_flux_benthslar_fishm

#............................

calc_flux_benthclar_fishm_o<- f1(benthclar_o,(fishm_o),u_fishm*PREF_benthclar_fishm,h_fishm*(volume_so+volume_d))
calc_flux_benthclar_fishm_i<- f1(benthclar_i,(fishm_i),u_fishm*PREF_benthclar_fishm,h_fishm*(volume_si))

calc_flux_benthclar_fishm_w <- calc_flux_benthclar_fishm_o + calc_flux_benthclar_fishm_i

model_flux_benthclar_fishm<-extract_model_flux(flow_matrix,"benthclar","mfish")

#calc_flux_benthclar_fishm_w
#model_flux_benthclar_fishm

#............................

calc_flux_fishplar_fishm_o<- f1(fishplar_o,(fishm_o),u_fishm*PREF_fishplar_fishm,h_fishm*(volume_so+volume_d))
calc_flux_fishplar_fishm_i<- f1(fishplar_i,(fishm_i),u_fishm*PREF_fishplar_fishm,h_fishm*(volume_si))

calc_flux_fishplar_fishm_w <- calc_flux_fishplar_fishm_o + calc_flux_fishplar_fishm_i

model_flux_fishplar_fishm<-extract_model_flux(flow_matrix,"pfishlar","mfish")

#calc_flux_fishplar_fishm_w
#model_flux_fishplar_fishm

#............................

calc_flux_fishdlar_fishm_o<- f1(fishdlar_o,(fishm_o),u_fishm*PREF_fishdlar_fishm,h_fishm*(volume_so+volume_d))
calc_flux_fishdlar_fishm_i<- f1(fishdlar_i,(fishm_i),u_fishm*PREF_fishdlar_fishm,h_fishm*(volume_si))

calc_flux_fishdlar_fishm_w <- calc_flux_fishdlar_fishm_o + calc_flux_fishdlar_fishm_i

model_flux_fishdlar_fishm<-extract_model_flux(flow_matrix,"dfishlar","mfish")

#calc_flux_fishdlar_fishm_w
#model_flux_fishdlar_fishm

#............................
#............................

calc_flux_corpse_fishd_o<- f1(corpse_o,(fishd_o),u_fishd*PREF_corpse_fishd,h_fishd*(volume_so+volume_d))
calc_flux_corpse_fishd_i<- f1(corpse_i,(fishd_i),u_fishd*PREF_corpse_fishd,h_fishd*(volume_si))

calc_flux_corpse_fishd_w <- calc_flux_corpse_fishd_o + calc_flux_corpse_fishd_i

model_flux_corpse_fishd<-extract_model_flux(flow_matrix,"corpses","dfish")

#calc_flux_corpse_fishd_w
#model_flux_corpse_fishd

#............................

calc_flux_discard_fishd_o<- f1(discard_o,(fishd_o),u_fishd*PREF_discard_fishd,h_fishd*(volume_so+volume_d))
calc_flux_discard_fishd_i<- f1(discard_i,(fishd_i),u_fishd*PREF_discard_fishd,h_fishd*(volume_si))

calc_flux_discard_fishd_w <- calc_flux_discard_fishd_o + calc_flux_discard_fishd_i

model_flux_discard_fishd<-extract_model_flux(flow_matrix,"discards","dfish")

#calc_flux_discard_fishd_w
#model_flux_discard_fishd

#............................

calc_flux_carn_fishd_o<- f1(CZ_edible_o,(fishd_o),u_fishd*PREF_carn_fishd,h_fishd*(volume_so+volume_d))
calc_flux_carn_fishd_i<- f1(CZ_edible_i,(fishd_i),u_fishd*PREF_carn_fishd,h_fishd*(volume_si))

calc_flux_carn_fishd_w <- calc_flux_carn_fishd_o + calc_flux_carn_fishd_i

model_flux_carn_fishd<-extract_model_flux(flow_matrix,"carnzoo","dfish")

#calc_flux_carn_fishd_w
#model_flux_carn_fishd

#............................

calc_flux_benths_fishd_o<- f1(benths_o,(fishd_o),u_fishd*PREF_benths_fishd,h_fishd*(volume_so+volume_d))
calc_flux_benths_fishd_i<- f1(benths_i,(fishd_i),u_fishd*PREF_benths_fishd,h_fishd*(volume_si))

calc_flux_benths_fishd_w <- calc_flux_benths_fishd_o + calc_flux_benths_fishd_i

model_flux_benths_fishd<-extract_model_flux(flow_matrix,"benths","dfish")

#calc_flux_benths_fishd_w
#model_flux_benths_fishd

#............................

calc_flux_benthc_fishd_o<- f1(benthc_o,(fishd_o),u_fishd*PREF_benthc_fishd,h_fishd*(volume_so+volume_d))
calc_flux_benthc_fishd_i<- f1(benthc_i,(fishd_i),u_fishd*PREF_benthc_fishd,h_fishd*(volume_si))

calc_flux_benthc_fishd_w <- calc_flux_benthc_fishd_o + calc_flux_benthc_fishd_i

model_flux_benthc_fishd<-extract_model_flux(flow_matrix,"benthc","dfish")

#calc_flux_benthc_fishd_w
#model_flux_benthc_fishd

#............................

calc_flux_fishp_fishd_o<- f1(fishp_o,(fishd_o),u_fishd*PREF_fishp_fishd,h_fishd*(volume_so+volume_d))
calc_flux_fishp_fishd_i<- f1(fishp_i,(fishd_i),u_fishd*PREF_fishp_fishd,h_fishd*(volume_si))

calc_flux_fishp_fishd_w <- calc_flux_fishp_fishd_o + calc_flux_fishp_fishd_i

model_flux_fishp_fishd<-extract_model_flux(flow_matrix,"pfish","dfish")

#calc_flux_fishp_fishd_w
#model_flux_fishp_fishd

#............................

calc_flux_fishd_fishd_o<- f1(fishd_o,(fishd_o),u_fishd*PREF_fishd_fishd,h_fishd*(volume_so+volume_d))
calc_flux_fishd_fishd_i<- f1(fishd_i,(fishd_i),u_fishd*PREF_fishd_fishd,h_fishd*(volume_si))

calc_flux_fishd_fishd_w <- calc_flux_fishd_fishd_o + calc_flux_fishd_fishd_i

model_flux_fishd_fishd<-extract_model_flux(flow_matrix,"dfish","dfish")

#calc_flux_fishd_fishd_w
#model_flux_fishd_fishd

#............................

calc_flux_fishplar_fishd_o<- f1(fishplar_o,(fishd_o),u_fishd*PREF_fishplar_fishd,h_fishd*(volume_so+volume_d))
calc_flux_fishplar_fishd_i<- f1(fishplar_i,(fishd_i),u_fishd*PREF_fishplar_fishd,h_fishd*(volume_si))

calc_flux_fishplar_fishd_w <- calc_flux_fishplar_fishd_o + calc_flux_fishplar_fishd_i

model_flux_fishplar_fishd<-extract_model_flux(flow_matrix,"pfishlar","dfish")

#calc_flux_fishplar_fishd_w
#model_flux_fishplar_fishd

#............................

calc_flux_fishdlar_fishd_o<- f1(fishdlar_o,(fishd_o),u_fishd*PREF_fishdlar_fishd,h_fishd*(volume_so+volume_d))
calc_flux_fishdlar_fishd_i<- f1(fishdlar_i,(fishd_i),u_fishd*PREF_fishdlar_fishd,h_fishd*(volume_si))

calc_flux_fishdlar_fishd_w <- calc_flux_fishdlar_fishd_o + calc_flux_fishdlar_fishd_i

model_flux_fishdlar_fishd<-(extract_model_flux(flow_matrix,"dfishlar","dfish")) - Dfish_recruit  # Flow matrix includes recruitment flux as well

#calc_flux_fishdlar_fishd_w
#model_flux_fishdlar_fishd


#............................

calc_flux_fishm_fishd_o<- f1(fishm_o,(fishd_o),u_fishd*PREF_fishm_fishd,h_fishd*(volume_so+volume_d))
calc_flux_fishm_fishd_i<- f1(fishm_i,(fishd_i),u_fishd*PREF_fishm_fishd,h_fishd*(volume_si))

calc_flux_fishm_fishd_w <- calc_flux_fishm_fishd_o + calc_flux_fishm_fishd_i

model_flux_fishm_fishd<-extract_model_flux(flow_matrix,"mfish","dfish")

#calc_flux_fishm_fishd_w
#model_flux_fishm_fishd

#............................
#............................



#Implement testthat checks...

expect_equal(model_flux_omni_fishp,calc_flux_omni_fishp_w , 1e-7)
expect_equal(model_flux_carn_fishp,calc_flux_carn_fishp_w , 1e-7)
expect_equal(model_flux_benthslar_fishp,calc_flux_benthslar_fishp_w , 1e-7)
expect_equal(model_flux_benthclar_fishp,calc_flux_benthclar_fishp_w , 1e-7)
expect_equal(model_flux_fishplar_fishp,calc_flux_fishplar_fishp_w , 1e-7)
expect_equal(model_flux_fishdlar_fishp,calc_flux_fishdlar_fishp_w , 1e-7)

expect_equal(model_flux_omni_fishm,calc_flux_omni_fishm_w , 1e-7)
expect_equal(model_flux_carn_fishm,calc_flux_carn_fishm_w , 1e-7)
expect_equal(model_flux_benthslar_fishm,calc_flux_benthslar_fishm_w , 1e-7)
expect_equal(model_flux_benthclar_fishm,calc_flux_benthclar_fishm_w , 1e-7)
expect_equal(model_flux_fishplar_fishm,calc_flux_fishplar_fishm_w , 1e-7)
expect_equal(model_flux_fishdlar_fishm,calc_flux_fishdlar_fishm_w , 1e-7)

expect_equal(model_flux_corpse_fishd,calc_flux_corpse_fishd_w , 1e-7)
expect_equal(model_flux_discard_fishd,calc_flux_discard_fishd_w , 1e-7)
expect_equal(model_flux_carn_fishd,calc_flux_carn_fishd_w , 1e-7)
expect_equal(model_flux_benths_fishd,calc_flux_benths_fishd_w , 1e-7)
expect_equal(model_flux_benthc_fishd,calc_flux_benthc_fishd_w , 1e-7)
expect_equal(model_flux_fishplar_fishd,calc_flux_fishplar_fishd_w , 1e-7)
expect_equal(model_flux_fishdlar_fishd,calc_flux_fishdlar_fishd_w , 1e-7)
expect_equal(model_flux_fishp_fishd,calc_flux_fishp_fishd_w , 1e-7)
expect_equal(model_flux_fishm_fishd,calc_flux_fishm_fishd_w , 1e-7)
expect_equal(model_flux_fishd_fishd,calc_flux_fishd_fishd_w , 1e-7)

})