Brandon's Productivity.R Homework (09/4/2012)

Code

Set some phytosysnthesis physiology parameters


Ik = 15  ###unit is W/m^2)
Pmax = 1e-04  ###unit is mg*C/(hr*cell)

#### respiration term per cell
R <- 1e-05  ####unit is mg*C/(hr*cell)

##### Since Pmax is per cell, then we need to have the number of cells at
##### each depth and time to compute productivity
cell.field <- matrix(10000, nrow(uI.field), ncol(uI.field))  ##We will hold the number of cells per grid cell constant for this exercise

prod.field <- cell.field * Pmax * tanh(uI.field/Ik) - R * cell.field

Plotting the light field vs the producivity field

par(mar = c(2.8, 2.8, 1, 1), mgp = c(2, 0.5, 0))
plot(uI.field, prod.field, xlim = c(1, 2), ylim = c(-0.1, 0.1), las = 1)
abline(h = 0, lty = 2)
abline(v = 1.5, col = "red")
text(1.6, 0.05, "~1.5 w/m^2")

plot of chunk Plot1

Graphically, this appears to be ~ 1.5 W/m2

Making a contour plot of the productivity field

par(mar = c(2, 2, 1, 1), mgp = c(1, 0.5, 0))
contour(sun_time, z, prod.field, las = 1, ylim = c(-40, 0))
abline(h = -30, lty = 2, col = "red")
text(1104628141, -37.94193, "~30 m")

plot of chunk Plot2

What is the approximate compensation depth of this system?

Integrating the total production

#### note, this only works using the sum() function because the
#### producivity is in units of per hour, and the time step of the model
#### is in hours.
sum(prod.field)
## [1] -501

Questions

1. What is the total production with the above attenuation, light levels, and physiology?

2. Is this a net autotrophic or net heterotrphic system?

3. At what Ik level is the whole system producing as much as it is respiring?


Productivity <- function(cell.field, Pmax, uI.field, Ik, R) {
    cell.field * Pmax * tanh(uI.field/Ik) - R * cell.field
}
net.auto <- {
}
for (Ik in seq(3, 4, by = 0.01)) {
    net.auto <- rbind(net.auto, c(Ik, sum(Productivity(cell.field, Pmax, uI.field, 
        Ik, R))))
}
net.auto[3:8, ]
##      [,1]    [,2]
## [1,] 3.02  4.0622
## [2,] 3.03  2.9731
## [3,] 3.04  1.8878
## [4,] 3.05  0.8061
## [5,] 3.06 -0.2720
## [6,] 3.07 -1.3465
net.auto <- {
}
for (Ik in seq(3.05, 3.07, by = 0.001)) {
    net.auto <- rbind(net.auto, c(Ik, sum(Productivity(cell.field, Pmax, uI.field, 
        Ik, R))))
}
net.auto[7:9, ]
##       [,1]     [,2]
## [1,] 3.056  0.15879
## [2,] 3.057  0.05104
## [3,] 3.058 -0.05668