  > library(tmdl)
  > data(pol)
  > data(pol)
  > flow.ranked <- rankflow(flow, pol, names = list(date = "Date", Q = "cfs"))
  > fd.fig <- fdplot(flow.ranked, names = list(Q = "cfs"), values = TRUE)
  > WQS <- 5
  > WLA.construction = function(exc, LC, MOS, LC.fun, MOS.fun) {
  +     LA = ifelse(exc > 75, LC.fun(75) - MOS.fun(75), LC - MOS)
  +     0.01 * LA
  + }
  > MS4.Duluth = function(exc, LC, MOS, LC.fun, MOS.fun) {
  +     0.16 * (LC - MOS - WLA.construction(exc, LC, MOS, LC.fun, MOS.fun))
  + }
  > WLA = list(WLA = WLA.construction, MS4 = MS4.Duluth)
  > kr.tmdl <- tmdl(flow.ranked, names = list(Q = "cfs", pol = "TSS"), WQS = WQS, WLA = WLA)
  > summary(kr.tmdl)

              High flows Moist conditions Mid-range flows Dry conditions  Low flows
  LC          5211.20688       819.602960      269.123360    117.4356480 39.1452160
  MOS         2546.88562       428.150800       73.397280     57.4945360 19.1648453
  90th      165781.70236      4190.984688      158.538125    112.0531808 11.8658936
  WLA           26.64321         3.914522        1.957261      0.5994111  0.5994111
  MS4          422.02849        62.006022       31.003011      9.4946721  3.1009535
  LA          2664.32126       391.452160      195.726080     59.9411120 19.9803707
  Reduction     98.39287        90.659661              NA     46.5065502         NA

  > ld.fig <- plot(kr.tmdl, ylab = expression(paste("Load , ", over(kg, day))), plot = list(WLA = TRUE))
  > mv.fig <- mvplot(flow.ranked)
  > ld2.fig <- update(ld.fig, par.settings = theEconomist.theme())

> plot(fd.fig)

Figure 1:

Flow duration based on USGS

> plot(ld.fig)

Figure 2:

Load duration

> plot(mv.fig)

Figure 3:

Monthly variation

> plot(ld2.fig)

Figure 4:

Load duration with custom theme