box_model.R

library(deSolve)
rho <- 1000                            # kg / m^3
Cp <- 4100                             # J / kg / degC
H <- 1                                 # water depth
B <- 100                               # noontime flux
t <- seq(0, 2*86400, 3600/4)           # 2 days, half-hourly


## No need to edit below this line
day <- t / 86400                       # days, for plots
flux <- function(t)
{
    f0 <- B * sin(2 * pi * t / 86400)
    f <- ifelse(f0 < 0, -B * 2 / pi, f0)
    f
}
dTdt <- function(t, y, parms) list(flux(t) / rho / Cp / parms$H)
par(mfrow=c(3,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
plot(day, flux(t), ylab="Flux [W/m^2]", type='l')
abline(h=0, col='blue')
grid()

sol <- lsoda(10, t, dTdt, parms=list(H=1))
plot(day, sol[,2], ylab="T [degC]", type='l')
grid()
mtext("H=1 m", line=-2)

sol <- lsoda(10, t, dTdt, parms=list(H=10))
plot(day, sol[,2], ylab="T [degC]", type='l')
grid()
mtext("H=10 m", line=-2)