# ==============================================================================
# Didactic Example: Simple Lake Phytoplankton Model
# ==============================================================================



# Load class definitions:
# =======================

source("ecosim.r")


# ==============================================================================
# Model with constant external influence factors
# ==============================================================================


# Definition of parameters:
# =========================


param    <- list(k.gro.ALG   = 1.0,      # 1/d   # average growth rate
                 k.death.ALG = 0.8,      # 1/d   # death rate includes grazing by zooplankton
                 K.HPO4      = 0.002,    # gP/m3
                 alpha.P.ALG = 0.002,    # gP/gDM
                 A           = 8.5e+006, # m2
                 h.epi       = 4,        # m
                 Q.in        = 4,        # m3/s
                 C.ALG.ini   = 0.05,     # gDM/m3
                 C.HPO4.ini  = 0.02,     # gP/m3
                 C.HPO4.in   = 0.04,     # gP/m3             
                 J.col.ALG   = 1000)     # gDM/d             


# Definition of transformation processes:
# =======================================

# Growth of algae:
# ----------------

gro.ALG   <- new(Class  = "process",
                 name   = "Growth of algae",
                 rate   = expression(k.gro.ALG
                                     *C.HPO4/(K.HPO4+C.HPO4)
                                     *C.ALG),
                 stoich = list(C.ALG  = expression(1),              # gDM/gDM
                               C.HPO4 = expression(-alpha.P.ALG)))  # gP/gDM

# Death of algae:
# ---------------

death.ALG <- new(Class = "process",
                 name   = "Death of algae",
                 rate   = expression(k.death.ALG*C.ALG),
                 stoich = list(C.ALG  = expression(-1)))            # gDM/gDM


# Definition of reactor:
# ======================

# Epilimnion:
# -----------

epilimnion <- 
   new(Class            = "reactor",
       name             = "Epilimnion",
       volume.ini       = expression(A*h.epi),
       conc.pervol.ini  = list(C.HPO4 = expression(C.HPO4.ini),     # gP/m3
                               C.ALG  = expression(C.ALG.ini)),     # gDM/m3
       input            = list(C.ALG  = expression(J.col.ALG)),     # gDM/d
       inflow           = expression(Q.in*86400),                   # m3/d
       inflow.conc      = list(C.HPO4 = expression(C.HPO4.in),
                               C.ALG  = 0),
       outflow          = expression(Q.in*86400),
       processes        = list(gro.ALG,death.ALG))


# Definition of system:
# =====================

# Lake system:
# ------------

lake <- new(Class    = "system",
            name     = "Lake",
            reactors = list(epilimnion),
            param    = param,
            t.out    = seq(0,365,by=1),
            dt.max   = 0.2)


# Perform simulation:
# ===================

res1 <- calc.res(lake)


# Plot results:
# =============
# the function plot.res has the following arguments: res, colnames, file, 
# if only res is given, each column of res is plotted in a seperate plot
# if colnames is given, each element of colnames is plotted in a seperate plot
# if a filename is specified with file="filename.pdf" the plot is written
# to a pdf file.

                 
plot.res(res1)               # plot to screen

plot.res(res=res1, colnames=list("C.HPO4","C.ALG"))

plot.res(res      = res1,    # plot to pdf file
         colnames = list("C.HPO4","C.ALG"),
         file     = "didacticmodel_lakephyto_1.pdf",
         width    = 8,
         height   = 4)


# ==============================================================================
# Change of Parameter Values
# ==============================================================================


#lake@param$k.gro.ALG <- 2
#res1 <- calc.res(lake)                  
#plot.res(res=res1, colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
#lake@param <- param


#lake@param$k.gro.ALG <- 0.5
#res1 <- calc.res(lake)                
#plot.res(res=res1, colnames=list("C.HPO4",c("C.ALG", "C.ZOO")))
#lake@param <- param



# ==============================================================================
# Model with seasonally varying external influence factors
# ==============================================================================


# Growth of algae extended process:
# ---------------------------------

gro.ALG.ext <-
   new(Class  = "process",
       name   = "Growth of algae extended",
       rate   = expression(k.gro.ALG
                           *exp(beta.ALG*(T-T0))
                           *C.HPO4/(K.HPO4+C.HPO4)
                           *log((K.I+I0)
                                /(K.I+I0*exp(-(lambda.1+lambda.2*C.ALG)*h.epi)))
                            /((lambda.1+lambda.2*C.ALG)*h.epi)
                           *C.ALG),
       stoich = list(C.ALG  = 1,                            # gDM/gDM
                     C.HPO4 = expression(-alpha.P.ALG)))    # gP/gDM

# Change processes in the reactor "epilimnion"
epilimnion@processes <- list(gro.ALG.ext,death.ALG)


# Changed parameter values and additional parameters:
# ===================================================


param    <- list(k.gro.ALG   = 1.5,      # 1/d   # maximum growth rate
                 k.death.ALG = 0.8,      # 1/d   # death rate includes grazing by zooplankton
                 K.HPO4      = 0.002,    # gP/m3
                 alpha.P.ALG = 0.002,    # gP/gDM
                 A           = 8.5e+006, # m2
                 h.epi       = 4,        # m
                 Q.in        = 4,        # m3/s
                 C.ALG.ini   = 0.05,     # gDM/m3
                 C.HPO4.ini  = 0.02,     # gP/m3
                 C.HPO4.in   = 0.04,     # gP/m3
                 J.col.ALG   = 1000,     # gDM/d             
                 beta.ALG    = 0.046,    # 1/degC
                 T0          = 20,       # degC
                 K.I         = 30,       # W/m2
                 lambda.1    = 0.05,     # 1/m
                 lambda.2    = 0.03)     # m2/gDM 
                 

# Add additional parameters for environmental conditions:
# =======================================================

param <- c(param,
           list(t.max  = 200,
                I0.min = 50,
                I0.max = 250,
                T.min  = 5,
                T.max  = 20))

lake@param       <- param
  

epilimnion@cond <- list(I0 = expression(0.5*(I0.min+I0.max)+
                                        0.5*(I0.max-I0.min)*
                                        cos(2*pi/365.25*(t-t.max))),   # W/m2
                        T  = expression(0.5*(T.min+T.max)+
                                        0.5*(T.max-T.min)*
                                         cos(2*pi/365.25*(t-t.max))))  # degC

epilimnion@input$C.ALG <- param$J.col.ALG

lake@reactors    <- list(epilimnion)


t.out <- seq(0,1461,by=1)

lake@t.out       <- t.out


# Perform simulation:
# ===================

res2 <- calc.res(lake)
       

# Plot results:
# =============
       
plot.res(res2)                       # plot to screen

plot.res(res=res2, colnames=list("C.HPO4","C.ALG"))

plot.res(res      = res2,    # plot to pdf file
         colnames = list("C.HPO4","C.ALG"),
         file     = "didacticmodel_lakephyto_2.pdf",
         width    = 8,
         height   = 4)