#-----------------------#
  # the model parameters: #
  #-----------------------#
  
  parameters<-c(maxPhotoSynt   =0.125,      #molC/molC/hr
                rMortPHY       =0.001,      #/hr
                alpha          =-0.125/150, #uEinst/m2/s/hr
                pExudation     =0.0,        #-
                maxProteinSynt =0.136,      #molC/molC/hr
                ksDIN          =1.0,        #mmolN/m3
                minpLMW        =0.05,       #molC/molC
                maxpLMW        =0.15,       #molC/molC
                minQuotum      =0.075,      #molC/molC
                maxStorage     =0.23,       #/h
                respirationRate=0.0001,     #/h
                pResp          =0.4,        #-
                catabolismRate =0.06,       #/h
                dilutionRate   =0.01,       #/h
                rNCProtein     =0.2,        #molN/molC
                inputDIN       =10.0,       #mmolN/m3
                rChlN          =1,          #gChl/molN
                parMean        =250.,       #umolPhot/m2/s
                dayLength      =15.         #hours
                )
  
  #-------------------------#
  # the initial conditions: #
  #-------------------------#
  
  state     <-c(DIN     =6.,     #mmolN/m3
                PROTEIN =20.0,   #mmolC/m3
                RESERVE =5.0,    #mmolC/m3
                LMW     =1.0)    #mmolC/m3
  
  #----------------------#
  # RUNNING the model:   #
  #----------------------#
  
  times <-seq(0,24*20,1)
  
  out <- as.data.frame(aquaphy(times,state,parameters))
  
  #------------------------#
  # PLOTTING model output: #
  #------------------------#
  
  par(mfrow=c(2,2), oma=c(0,0,3,0))    
  col <- grey(0.9)
  ii <- 1:length(out$PAR)              
  
  plot (times[ii],out$Chlorophyll[ii],type="l",
        main="Chlorophyll",xlab="time, hours",ylab="ug/l")
  polygon(times[ii],out$PAR[ii]-10,col=col,border=NA);box()
  lines (times[ii],out$Chlorophyll[ii]  ,lwd=2 )
  
  
  plot (times[ii],out$DIN[ii]        ,type="l",main="DIN",
        xlab="time, hours",ylab="mmolN/m3")
  polygon(times[ii],out$PAR[ii]-10,col=col,border=NA);box()
  lines (times[ii],out$DIN[ii]  ,lwd=2 )
  
  
  plot (times[ii],out$NCratio[ii]    ,type="n",main="NCratio",
        xlab="time, hours",ylab="molN/molC")
  polygon(times[ii],out$PAR[ii]-10,col=col,border=NA);box()
  lines (times[ii],out$NCratio[ii]  ,lwd=2 )
  
  
  plot (times[ii],out$PhotoSynthesis[ii],type="l",
         main="PhotoSynthesis",xlab="time, hours",
         ylab="mmolC/m3/hr")
  polygon(times[ii],out$PAR[ii]-10,col=col,border=NA);box()
  lines (times[ii],out$PhotoSynthesis[ii]  ,lwd=2 )
  
  mtext(outer=TRUE,side=3,"AQUAPHY",cex=1.5)
  
  #------------------------#
  # SUMMARY  model output: #
  #------------------------#
  t(summary(out))Run the code above in your browser using DataLab