
# R program to make frames for raised peat-bog growth simulation

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#  You can run the program with 'R CMD BATCH ???.R' where '???' is the stem name  #
#  of the program e.g. 'Bog_sim_V6.R'. The program should run for several minutes #
#  without any indication of progress. It will produce a file ???.Rout (e.g.      #
#  'Bog_sim_V6.Rout' that lists messages. If it fails, an explanation should be   #
#  in this ???.Rout file                                                          #
#                                                                                 #
#  Or you can run the program with 'Rscript ???.R' (e.g. 'Rscript bog_sim_V6.R'). #
#  This mode will produce several bits of information, and will record progress   #
#  on the terminal.                                                               #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# R is meant for vector processing, and has good graphics abilities.
# I scarcely use vector processing here, but need the graphics
# Much of the program uses 'for' loops as it deals with complicated
# different actions for each of the cells in a vector.

# The program is unlikely to be run often, so I have made no attempt
# to make it efficient, or to remove redundancies.

# All the important settings can be changed: find 'GLOBAL' and read on

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#  Copyright (C) 2017 by R. S. Clymo                                #
#  r.clymo@QMUL.ac.uk                                               #
#                                                                   #
#  This program is free software; you can redistribute it and/or    #
#  modify it under the terms of version 2 of the GNU General        #
#  Public License as published by the Free Software Foundation.     #
#                                                                   #
#  This program is distributed in the hope that it will be useful,  #
#  but WITHOUT ANY WARRANTY; without even the implied warranty of   #
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    #
#  GNU General Public License for more details.                     #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# current settings, all changeable
# 25 fps for 20 s = 500 frames
# + intro (125) + coda (125) = 750 frames total;

# grow for 10000 year, so each frame is 10000 / 500 = 20 yr;

# for linear depth, not mass (but easily adapted to mass if dry bulk
# density profiles are available);

# choice of decay functions: Boudreau & Ruddick, exponential, constant

# TERMINOLOGY
# 'wide' means half the full width of the generally symmetrical peat mass
# 'high' means the vertical height, relative to the base.



# FUNCTION LIST:

# SERVICE
# NL, SP, trk, png.dev, time.name, tweak,
# STRUCTURES
#   decay, peat.sim, sim.try, plot.ellipse, sym.tree,
#   plot.growth, plot.ghosts,
#   calc.slopes, plot.pools, hyperb.tree.tall, plot.trees, plot.lagg, plot.woods,
#
# CONTROL
# run.intro, run.sim, ORCHESTRATE

# CONVENIENCE
NL            <- "\n"                       # for cat
SP            <- ' '

# SECTIONS are GLOBAL, FUNCTIONS, MAIN


# GLOBAL


# USER MUST CHOOSE

# for EX)panding across a plain or ST)ationary (everywhere at once)
# 'topog' must be "EXos" or "STos"
topog         <- "EXos"
topog         <- "STos"


# Decay models

# 1. Constant rate of decay (unrealistic)
# 2. Exponential decline in decay rate with time
# 3. BandR: Boudreau & Ruddick (1991) On a reactive continuum representation
#    of organicmatter diagenesis. American Journal of Science, 291, 507–538.

# Decay model: choose one of "cnst", "xpl", "BandR"
dk.model <- "xpl"
dk.model <- "cnst"
dk.model <- "BandR"

# USER may wish to change one or more of these settings

# tracking controls (for debugging)
from.frame    <- 9999 # large e.g. 9999 SUPPRESSES tracking
fcn.trkd      <- 0    # counts track calls
max.trkd      <- 200  # arbitrary track calls print limit

# Bog growth parameters

# p.1 productivity;    vp.2 vary p.1;
# a.3 decay;           va.4 vary decay
# tlim.5 years to run; dt.6 evaluate every years

p.1             <- NA # calculated later to end at specified depth and bog width
vp.2            <- 0  # '0' inactivates exponential decline in 'p.1'
                      # -0.002 -> convex age vs depth

if        (dk.model == "cnst") {
	a.3           <- 0.98 * 1.0e-4 # avoid error at t=10000
	va.4          <- 0

} else if	(dk.model == "xpl") {
	a.3           <- 0.98 * 1.0e-4
	va.4          <- 0.0001

} else if (dk.model == "BandR") {
	a.3           <- 100
	va.4          <- 0.5

} else {
  cat ("decay model '", dk.model, "not one of 'cnst', 'xpl', 'BandR'", NL)
  quit ()
}

tlim.5          <- 10000
dt.6            <- 20;


# other bog parameters


# 'bog.wide.mx', 'bog.high.mx' are (vertical) ellipse semi radii 'a'(x) and 'b' (y)
bog.wide.mx     <- 500   # / m, half width
bog.high.mx     <- 7     # / m, depth at tlim.5
asp.x.ov.y      <- bog.wide.mx / bog.high.mx # aspect quotient ('ratio')

x.step          <- 2         # / m, x-axis steps, bog, lagg, woods
half.step       <- x.step / 2 # for plotting about a cell centre

# acrotelm plotting
acro.dpth       <- 0.35  # / m, acrotelm above peat surface
acro.lwd        <- 1 + x.step # plotting line width, equivalent to acro

# specify scre  en dimensins: wider is generally better
fnt.size        <- 10 # font size
scrn.high       <-  5 # inches
scrn.wide       <- 18 # inches

# three phases
# the 'coda' section simply repeats the last frame of 'sim' for several s;
# it is managed inside 'run.sim' by setting 'change <- F'
do.intro        <- T
do.sim          <- T
do.coda         <- T

# simulate features
do.cato         <- T # basic peat surface
do.acro         <- T # basic acrotelm
do.ghosts       <- T # show old ('ghost') surfaces
do.pools        <- T # add pools
do.trees        <- T # pine/larch trees on bog
do.lagg         <- T # deciduous trees in lagg fen
do.woods        <- T # spruce, some pine, trees around bog

do.grth.grph    <- T # show growth at left side as it happens


# seconds for intro, simulation, still
intro.sec       <-  5
sim.sec         <- 20
coda.sec        <-  5

# Note. Some functions (e.g. plot.pools. plot.lagg, plot.woods) need
# thousands of preferably reproducible pseudo-random numbers. One or more
# of these functions may be switched off. A single PRNG then produces different
# sets of PRNs in the functions that are switched on. But there seems to
# be no simple way in R of running independent PRNGs for each function

prng.seed       <- 2567   # arbitrary, OK to change

# video controls

# set fps (frames per second)
fps             <- 25





# OTHER SETTINGS (should not need, or must not be, changing)

intro.frames    <- fps * intro.sec
sim.frames      <- fps * sim.sec
coda.frames     <- fps * coda.sec

tot.frames      <- NA

change          <- NA # later set to 'T' at first, but changed to 'F' for coda

cur.frame       <- 0  # important base setting for all 3 sections
bog.frm.ct      <- 0  # bog frame count
bog.frm.start   <- NA
frm.ct          <- 0  # running count
t.yr            <- NA # to be set later

# years increment at each step
t.frm.inc     <- tlim.5 / sim.frames


# corners for main and subsidiary (height vs time) plot position and size
#                   xlo xhi ylo yhi
main.crnrs      <- c (0.15, 1, 0, 1)
sbsid.crnrs     <- c (0,  0.2, 0, 1)

# stores for subsidiary plot (in each frame) of height vs time
height.at       <- rep (NA, sim.frames)       # height.at ...
time.n          <- rep (NA, sim.frames)       # age / yr


# separate storage for old heights
# NA is dummy to set store size; reset later
n.ghosts        <- 5 # equal times in tlim.5 e.g.10000 / 5 -> 2000
ghost.intvl     <- ceiling (tlim.5 / n.ghosts) # interval to save surface
# n.b. grey12 is dark, grey60 is light !
ghost.clr       <- c ("grey12", "grey24", "grey36", "grey48", "grey60")
ghost.yr        <- rep (-1, n.ghosts)
ghost.wide.0    <- rep (-1, n.ghosts) # width  when formed
ghost.high.0    <- rep (-1, n.ghosts) # height when formed
ghost.high.t    <- rep (-1, n.ghosts) # height at end
ghost.slb       <- rep (-1, n.ghosts)
ghost.ct        <- 0
ghost.nxt       <- ghost.intvl

# tree scaling factor: alters the shape of all trees on bog, lagg, woods
tree.high.fct   <- 1.0
tree.wide.fct   <- 1.0

tree.grth.fct   <- 0.3

# Common spaces, cells and distances  / m

# the maximum; actual for topog == "EXos" is less until end
bog.wide.cells     <- ceiling (bog.wide.mx / x.step)

if (do.lagg)  {
  stream.cells     <- 2
  lagg.tree.cells  <- 3
  lagg.wide.cells  <- stream.cells + lagg.tree.cells
  lagg.wide.mx     <- lagg.wide.cells * x.step
} else {
  stream.cells     <- 0
  lagg.tree.cells  <- 0
  lagg.wide.cells  <- 0
  lagg.wide.mx     <- 0
}
# lagg depth below mineral ground zero
lagg.stream   <- -0.3  # / m deep


if (do.woods) { # woods around the bog
  woods.wide.cells <- 15
  woods.wide.mx  <- woods.wide.cells * x.step
} else {
  woods.wide.cells <- 0
  woods.wide.mx    <- 0
}
woods.st          <- NA
woods.nd          <- NA

all.wide.cells <- bog.wide.cells + lagg.wide.cells + woods.wide.cells
all.wide.mx    <- bog.wide.mx  + lagg.wide.mx + woods.wide.mx # / m

# global common distance and bog cells (only) for variables v?: with uses
# depending on structure; all or some of these cells carry history
# from one frame to the next

# use as e.g.v3[1,?] and v3[2,?] for left -ve and right +ve distances
# from centre by  'for (rl in 1:2) { v3[rl,?] <- ??? }'


# store values of different meanings in different places
all.bog.cells <- 2 * bog.wide.cells
all.cells     <- 2 * all.wide.cells
vx            <- rep (0, all.cells)       # x-coords
vs            <- rep (99, bog.wide.cells) # slopes (right and left same)
vp            <- rep (0, all.cells)       # pools
vt            <- rep (0, all.cells)       # trees
v1            <- rep (0, all.cells)
v2            <- rep (0, all.cells)
v3            <- rep (0, all.cells)
v4            <- rep (0, all.cells)
v5            <- rep (0, all.cells)
v6            <- rep (0, all.cells)
# specify as a 2-row array
dim (vx)      <- c(2, all.wide.cells)
#    vs       is a single dimension filled with 99
dim (vp)      <- c(2, all.wide.cells)
dim (vt)      <- c(2, all.wide.cells)
dim (v1)      <- c(2, all.wide.cells)
dim (v2)      <- c(2, all.wide.cells)
dim (v3)      <- c(2, all.wide.cells)
dim (v4)      <- c(2, all.wide.cells)
dim (v5)      <- c(2, all.wide.cells)
dim (v6)      <- c(2, all.wide.cells)

# set x-coordinates
up.to         <- all.wide.mx
# xxx is distance from centre to MIDDLE of cell
xxx           <- seq (from=half.step, to=up.to, by=x.step)

vx[1,]        <-  xxx # right side
vx[2,]        <- -xxx # left  side

rm (xxx)

cat (6 * all.bog.cells, " cells needed; distance", all.wide.mx, " m/n")

# these are set in 'plot.pools' and 'plot.trees'
pool.cell.lim <- -1
cell.pool.nd  <- -1
cell.tree.nd  <- -1
cell.tree.st  <- -1

# skip.b4.frame <-  5 # no pools until 5th frame

# current width (x), height (y), and their squares; bog + lagg + woods
cur.high      <- NA
cur.wide      <- NA       # to be set later
cur.wide.sq   <- NA       # to be set later

# pools: space and critical slope (hydraulic gradient)
is.pool       <-  999     #is a pool
h.p.crit      <-  0.004   # hydraulic gradient (slope) less than this for pools
h.d.crit      <-  1.5 * h.p.crit # slope for pool to dry out
pool.penalty  <-  0.1     # factor reducing pool potential non-pools either side of pool
pool.risk     <-  300     # larger decreases chance of a pool forming
pool.lwd      <-  6       # pool width, equivalent to cell width

pool.dpth.st  <-  0.1     # proportion of depth at start
pool.dpth.nd  <-  2.0     # proportion of full depth at end
pool.dpth.inc <-  1.5 * acro.dpth * (pool.dpth.nd - pool.dpth.st) / sim.frames



# trees on bog perimeter
h.t.crit          <-   0.0222 # 1 in 45 slope
tree.mx.death.yr  <- 200      # tree max age at death
tree.high.at.death <-  12.0   # / m
tree.start.prb    <-   0.5    # probability that a new tree will start
tree.cell.overrun <-   5      # search beyond h.t.crit
ix.trees          <-   0      # saved index
new.trees.ct      <-   0


x.lim         <- c (-all.wide.mx, +all.wide.mx)
y.mx          <- high.mx + 1.0              # max y value for plot
y.mn          <- lagg.stream                 # lagg stream depth
y.lim         <- c (y.mn, y.mx)              # a bit more than high.mx

nodk.mx.yr    <- NA

#Plotting libraries
# N.B.: needs 'cairoDevice' and 'plotrix' libraries
library (cairoDevice)
library (plotrix)


# FUNCTIONS

trk <- function (txt="none") {
# lists functions called; for debugging
  if ((cur.frame >= from.frame) && (fcn.trkd < max.trkd)) {
    cat ("Frame", frm.ct, ":")
    if (txt == "none") {
      cat ("No function name in 'trk'. QUIT\n")
      quit (status=59)
    }
    fcn.trkd <<- fcn.trkd + 1
    cat ("In", txt, ": call #", fcn.trkd, NL)
  }
} # trk


png.dev <- function (out.f=NULL, wide=6, high=9, size=12) {
# set graphic device for 'png'; needs library 'cairoDevice'
# 'out.f' is stem name of output file
  graphics.off ()
  fname <- paste (out.f, ".png", sep="")
  Cairo (filename=fname, surface="png",
         width=wide, height=high, pointsize=size)
} # png.dev

time.name <- function (fr.num) {
# for ffmpeg names: allows frames 0 to 9999
# preceding '0's needed for compatibility with ffmpeg
  if      (fr.num >= 1000) fr.nam <- fr.num
  else if (fr.num >= 100)  fr.nam <- paste ("0",   fr.num, sep="")
  else if (fr.num >= 10)   fr.nam <- paste ("00",  fr.num, sep="")
  else                     fr.nam <- paste ("000", fr.num, sep="")
  return (fr.nam)
} # time.name

tweak <- function (vv=NA, prp=0.1) {
# tweak value 'vv' randomly by up to +/- half proportion 'prp'
  twk.vv <- vv * (1 + runif (1, max=prp) - prp / 2)
  return (twk.vv)
} # tweak


decay <- function (kind=NULL, par.1=NA, par.2=NA, tt=NA) {
  if        (kind == "BandR") {
    # Boudreau & Riddick
	  w <- (par.1 / (par.1 + tt)) ^ par.2
  } else if (kind == "xpl") {
    # exponential
    w <- exp (-par.2 * tt) * (1 - par.1 * tt)
  } else if (kind == "cnst") {
    # constant
    w <- 1 - par.1 * tt
  } else { # should never reach this: kind already checked
    cat ("\ndecay kind already checked: should be impossible", NL)
    quit ()
  }

  if ((w <= 0) || (w >= 1)) {
    cat ("\nERROR: decay factor outside 0+ to 1-:", w, NL)
    cat ("par.1, par.2, tt: ", par.1, SP, par.2, SP, tt, NL)
    quit (status=99)
  }

  return (w)

} # decay


peat.sim <- function (p=NULL, a=NULL, vp=NULL, va=NULL, t.mx=NULL,
                      dt=NULL, sav.ghost=T) {
# 4 parameter simulation for t.mx yr, each step dt yr, save ghost of
# surface every ghost.intvl yr

  trk ("peat.sim")

  # allocate space
  ns.mx <- ceiling (t.mx / dt) + 1 # number of slabs

  yr.made      <- rep (0, ns.mx) # year made (base = 0)
  high.0       <- rep (0, ns.mx) # height.at start
  high.t       <- rep (0, ns.mx) # height now
  dk.fct       <- rep (0, ns.mx) # decay factor
  cum.bog.wide <- rep (0, ns.mx) # cumulative bog width at each slab
  cum.bog.high <- rep (0, ns.mx) # cumulative height at each slab
  cum.nodk     <- rep (0, ns.mx) # cumulative height if no decay

  # set basal slab
  yr           <- dt / 2         # slab ages in centre of slab
  yr.made[1]   <- yr
  w            <- p * yr         # p * (dt / 2)
  high.0[1]    <- w;             # time at half dt

  # main loop: yr, time at half slab from base;
  for (ns in 2:ns.mx) {
    yr <- yr + dt                # next slab, yr from start
    yr.made[ns] <- yr
    v          <- exp (-vp * yr) # exponential diminishing p (if vp != 0)
    high.0[ns] <- v * p * dt;    # new top slab, no decay yet

    # decay factor; each slab follows the track of the (smaller) one
    # below, so hoist factors up one slab, and calculate new factor for
    # basal slab ONLY

    # (inefficient but is there a better method?)
    # dk.fct is decay factor
    for (j in ns:2) dk.fct[j] <- dk.fct[j-1] # hoist up, from top dowmwards

    # set new values in basal slab
    w          <- decay (dk.model, a.3, va.4, yr)
    dk.fct[1]  <- w
    high.t[1]  <- high.0[1] * w # height at current slab


    cum       <- 0.0                # height of top slab; will be returned
    cum.ndk   <- 0.0                # height if no decay
    for (j in 1:ns) {               # for slabs up to slab to date (in ns.mx)
      u <- high.0[j]                # initial height of this slab
      w <- dk.fct[j] * u            # current height of this slab
      high.t[j]    <- w             # save current height of this slab
      cum          <- cum + w       # cumulative height so far
      cum.bog.high[j]  <- cum       # save cumulative height so far (for later use)
      cum.ndk      <- cum.ndk + u   # cumulative height if no decay
      cum.nodk[j]  <- cum.ndk       # save cumulative height if no decay
      if (cum.ndk < high.mx) nodk.mx.yr <- yr
    } # for j

    # save ghost?
    if (yr >= ghost.nxt) {
      if (sav.ghost) {
        ghost.ct <<- ghost.ct + 1
        # is this is the first time this height has been reached?
        if (ghost.high.0[ghost.ct] < 0) { # it was set to -1 in global

          # set fixed values, once only
          # width when formed
          if (topog == "STos") { # full width, stationary
            ghost.wide.0[ghost.ct]  <<- bog.wide.mx
          } else {               # topog == "EXos"; expanding
            ghost.wide.0[ghost.ct]  <<- cum * asp.x.ov.y      # width when formed
          }
          ghost.high.0[ghost.ct]  <<- cum                     # height when formed
          ghost.slb[ghost.ct]     <<- ns  # slab number when formed
        }

        ghost.yr[ghost.ct]        <<- ghost.nxt               # idealised yr
        ghost.nxt                 <<- ghost.nxt + ghost.intvl # next ghost
      }
    }

  } # for ns

  # height at end of growth
  if (ghost.ct > 0) {
    for (j in 1:ghost.ct) {
      k <- ghost.slb[j]
      ghost.high.t[j] <<- cum.bog.high[k]
    }
  }

  rm (yr.made, high.0, high.t, dk.fct, cum.bog.high)

  return (cum)

} # peat.sim

sim.try <- function (pp) {
# prepare for uniroot (zero of a function)
# calc difference of final depth from that set in high.mx
  trk ("sim.try")

  sims <- peat.sim (pp, a.3, vp.2, va.4, tlim.5, dt.6, sav.ghost=F)
  res  <- sims - high.mx # the value to be minimised by uniroot
  return (res)
} # sim.try


plot.ellipse <- function (height=NULL, offs=0.0, width=NULL, colr=NULL,
                          thick=NULL) {
# plot an ellipse
  trk ("    plot.ellipse")
  width.sq <- width * width
  # points on ellipse: / m
  xx <- seq (from=0.0, to=width, by=x.step)
  ww <- xx^2 / width.sq
  if ((is.na (ww)) || (ww > 1.0)) return ()  # skip
  yy <- height * sqrt ( 1 - ww) + offs

  # negative half (to left)
  points (x=-xx, y=yy, type="l", lwd=thick, col=colr)
  # positive half (to right)
  points (x=xx,  y=yy, type="l", lwd=thick, col=colr)

  # tidy up
  rm (xx, yy)

} # plot.ellipse



sym.tree <- function (x=NA, y=NA, kind="T", high=1.0, y.size=1.0, x.size=1.0,
                      trk.top=0.4, trk.wyd=NA, trk.wyd.ov.high=0.04,
                      trk.prpy=-2,
                      cpy.bot=0.2, cpy.wyd=NA, cpy.wyd.ov.high=0.2,
                      cpy.prpy=2, from=NULL) {
# draw 1 of 3 kinds of stylised (symbolic) tree

# GENERAL PROPERTIES of all types of tree
# x, y coordinates of tree trunk base
# kind = one of "T"riangle, "E"llipse, "Z"ig-zag (spruce, broadleaf, pine)
# high           = overall height / y-units
# y.size         = factor to multiply height by
# x.size         = factor to multiply widths by

# TRUNK properties
# trk.top        = proportion of height (range 0 to <1)
# trk.wyd        = trunk width / x-units
# trk.wyd.ov.high = trunk width / height
# trk.prpy       = 1 - 4 fill with "brown[1-4]";
#                  5 - 8 fill with "grey[20, 40, 60, 80]";
#                  0, 9 no fill;
#                  +ve rectangular; -ve taper towards point at top of tree

# CANOPY properties
# cpy.bot        = canopy lower wide as proportion of height (0 to <1))
# cpy.wyd        = canopy width
# cpy.wyd.ov.high = canopy width / height
# cpy.prpy       = 1 - 4 fill with "green[1-4]";
#                  5 - 8 fill with "grey[20, 40, 60, 80]";
#                  0, 9 no fill;
#                  only +ve accepted

# trk.wyd, cpy.wyd default to NA. If supplied they have priority over
# trk.wyd.ov.high, cpy.wyd.ov.high (which specify shape as constant
# proportions)

# NB few checks! Caveat emptor!
  trk ("    sym.tree")

  if ((is.na (x)) || (length (x) != 1)) {
    cat ("x:", x, NL)
    cat ("kind, high, cpy.wyd:", kind, high, cpy.wyd, NL)
    cat ("'x' must be single value\n")
    cat ("Called 1 from:", from, NL)
    quit (status=80)
  }

  # COMMON to canopy and trunk
  high         <- y.size  * high   # overall height / y-units


  # CANOPY shape depends on kind
  cpy.lo      <- cpy.bot * high    # base of canopy above y / y-units
  if (!is.na (cpy.wyd)) w <- cpy.wyd
  else                  w <- cpy.wyd.ov.high * high
  cpy.wyd     <- x.size * w       # canopy width / x-units
  cpy.wide.wyd <- 0.5 * cpy.wyd    # canopy half-width / x-units

  # canopy colour green/grey and density (for all 3 kinds)
  if ((cpy.prpy == 9) || (cpy.prpy == 0)) cpy.clr <- NULL
  else {
    if (cpy.prpy < 5) {
      colr <- "green"
      num  <- cpy.prpy
    } else if (cpy.prpy < 9) {
      colr <- "grey"
      num  <- 20 * (9 - cpy.prpy)
    } else {
      cat ("cpy.prpy", cpy.prpy, "must be in range 0 to +9 \n")
      cat ("Called 2 from:", from, NL)
      quit (status=98)
    }
    cpy.clr <- paste (colr, as.character (num), sep="")
  }

  if (kind == "T") { # a triangle, cf. spruce
    # SW, N, SE, SW
    xx  <- x + c (-cpy.wide.wyd, 0, +cpy.wide.wyd, -cpy.wide.wyd)
    w   <- y + cpy.lo
    yy  <- w + c (0, high, 0, 0)
    polygon (x=xx, y=yy, col=cpy.clr, border=NA)

  } else if (kind == "E") { # an ellipses cf. broadleaf
    a.x <- cpy.wide.wyd
    b.y <- 0.5 * (high - cpy.lo)
    xx  <- x
    yy  <- y + cpy.lo + b.y # ellipse centre
    draw.ellipse (x=xx, y=yy, a=a.x, b=b.y, nv=30, arc.only=T,
                  col=cpy.clr)

  } else if (kind == "Z") { # zig-zag line cf. pine
    if (cpy.prpy == 2) cpy.prpy <- 5 # override default
    pts <- 2 * cpy.prpy + 1 # override default (= 5)
    if ((cpy.prpy < 4) || (cpy.prpy > 9)) {
      cat ("cpy.prpy", cpy.prpy, " is outside 4 to 9 for kind = 'Z'\n")
      cat ("Called 3 from:", from, NL)
      quit (status=97)

    }
    xx  <- rep (NA, pts)
    yy  <- xx # yy also NA
    del.y <- (high - cpy.lo) / (pts - 1)
    for (j in 1:pts) yy[j] <- (j-1) * del.y
    yy <- y + cpy.lo + yy
    del.x <- cpy.wide.wyd / (pts-1)
    dir <- -1
    x.cur <- cpy.wide.wyd
    for (j in 1:pts) {
      xx[j] <- dir * x.cur
      x.cur <- x.cur - del.x
      dir   <- -dir
    }
    xx      <- xx + x
    xx[1]   <- x # 1st point and ...
    xx[pts] <- x # ... last point are central
    points (x=xx, y=yy, type="l", col="black")

  } else {
    cat ("Tree kind '", kind, "' must be one in 'TEZ'\n" )
    cat ("Called 4 from:", from, NL)
    quit (status=96)
  }

  # trunk, being later, over-plots canopy

  # TRUNK  is common to all 3 kinds of tree
  trk.hi <- trk.top       * high      # top of trunk above y / y-units
  if (!is.na (trk.wyd)) w <- trk.wyd
  else                  w <- trk.wyd.ov.high * high
  trk.wyd <- x.size       * w        # trunk width / x-units

  wide.wyd.lo      <- 0.5 * trk.wyd
  if (trk.prpy > 0) wide.wyd.hi <- wide.wyd.lo # +ve = rectangular
  else              wide.wyd.hi <- wide.wyd.lo * (1 - (trk.hi / high))
  # NW, SW, SE, NE, SW
  xx <- x + c (-wide.wyd.hi, -wide.wyd.lo, +wide.wyd.lo, +wide.wyd.hi)
  yy <- y + c (trk.hi, 0, 0, trk.hi)
  # trunk colour brown/grey and density
  fill.code <- abs (trk.prpy)
  if ((fill.code == 9) || (fill.code == 0)) clr <- NULL
  else {
    if (fill.code < 5) {
      colr <- "brown"
      num  <- fill.code
    } else if (fill.code < 9) {
      colr <- "grey"
      num  <- 20 * (9 - fill.code)
    } else {
      cat ("trk.prpy", trk.prpy, "must be in range -9 to +9 \n")
      cat ("Called 5 from:", from, NL)
      quit (status=95)
    }
    clr <- paste (colr, as.character (num), sep="")
  }
  # draw the trunk
  polygon (x=xx, y=yy, border="black", col=clr)

} # sym.tree


plot.growth <- function () {
# supplementary plot of growth vs age
  trk ("  plot.growth")
  par(fig=sbsid.crnrs, new=T)
  plot (x=time.n, y=height.at, type="l",
        xlim=c(0,1.1*tlim.5), ylim=y.lim,
        yaxp=c(0,y.mx,4),
        xlab="", ylab="",
        col="red", bty="n", lwd=2, las=1)
  mtext ("Time / yr", side=1, line=+0.2)
  mtext ("Height / m",   side=2, line=+2)

  mlt <- 0.3 # omit lower 30 %
  at.t <- c (mlt * nodk.mx.yr, nodk.mx.yr)
  height <- c (mlt * high.mx, high.mx)
  points (x=at.t, y=height, type="l", lwd=1.5,  lty="dotted", col="red")

} # plot.growth


plot.ghosts <- function () {
# plot current position of old surfaces

  # place to put symbol for original position
  w.co <- c (-40, -22, +40, +20) # reversed 3,4 is OK

  ix <- 1
  ic <- ghost.ct

  while ((ix <= ghost.ct) && (ghost.yr[ix] < tlim.5)) { # exclude very last
    clr <- ghost.clr[ic] # colours from low to high index
    ic  <- ic - 1

    high.0 <- ghost.high.0[ix]
    high.t <- ghost.high.t[ix]
    wide.0 <- ghost.wide.0[ix]
    lin.wt <- 1.5

    # symbol for original surface top
    x.w <- w.co[ix]
    points (x=x.w, y=high.0, col=clr, type="p", pch=16, cex=1)


    # line to centre of current position
    xx <- c (x.w, 0)
    yy <- c (high.0, high.t)
    points (x=xx, y=yy, type="l", lwd=lin.wt, col=clr, lty=3)
    # lty=3 for dotted line

    # text age
    x.p <- 0.35 * cur.wide + 10 * ix # shift right
    y.p <- high.t - 0.4
    text (x=x.p,  y=y.p, col=clr,
          labels=paste (as.character(ghost.yr[ix]), " yr\n"), pos=2)

    # semi-ellipse of old surface as it is now
    plot.ellipse (height=high.t, offs=-0.1, width=wide.0, colr=clr, thick=lin.wt)

    ix <- ix + 1

  } # while

} # plot.ghosts


calc.slopes <- function () {
# needed for plot.pools and plot.trees
# vs values are 99 (unset) or slope

  if (change) {
    ix.lim <- ceiling (cur.wide / x.step)
	  ix <- 0
		repeat {
			ix <- ix + 1
			if ((ix > ix.lim) || (ix > bog.wide.cells)) break

			x  <- vx[1,ix]
			w  <- x^2 / cur.wide^2
			if ((is.na (w)) || (w > 1.0)) return ()      # skip
			y  <- cur.high * sqrt (1 - w)

			slope   <- (cur.high^2 * x) / (cur.wide^2 * y) # ellipse has 2 +ve and 2 -ve slopes
      vs[ix] <<- slope

      # is slope so steep that pools dry out?
      if (slope >= h.d.crit) {
        for (rl in 1:2) vp[rl,ix] <<- 0
        pool.cell.lim <<- ix
      }
		}

	} # if change

  return ()

} # calc.slopes


plot.pools <- function () {
# generate and plot pools: new version
  trk ("plot.pools")

  if (bog.frm.ct < skip.b4.frame) return () # skip early frames to avoid possible 0 indices

  # vx are cell distances from centre
  # vp 'is.pool' = 1; '!is.pool' = 0; saved
  # v3 pool potential
  # v1 pool depth


	if (change) {
		v3[1,] <- rep (0, all.wide.cells) # pool potential as 1 - slope / h.p.crit
	  v3[2,] <- rep (0, all.wide.cells)

	  # pool potential inversely proportional to slope, range 1:0 flat:h.p.crit
	  ix <- 0
		repeat {
			ix <- ix + 1
			if (ix > bog.wide.cells) return () # skip pools (should never happen)

			slope <- vs[ix]

	    # too steep?
			if (slope > h.p.crit) {
				ix     <- ix - 1                 # has overshot
			  if (ix < 1) return ()            # skip pools
			  cell.pool.nd <<- ix
				break
			}
			# save pool potential in v3, range 0:1
	    v <- 1 - slope / h.p.crit
			v3[1,ix] <- v # left  side
			v3[2,ix] <- v # right side
		}

		# factor by which to reduce potential of those cells adjacent to a pool
		# pool.penalty 0.1
		pen.ct <- 0        # penalty count; for info during development
		for (rl in 1:2) {  # right, then left side

	  	if (cell.pool.nd <= 3) break   # skip early frames only
			for (ix in 2:cell.pool.nd-1) { # not end cells
				if (vp[rl,ix] == is.pool) { # saved from frame to frame

					# to right of current pool
					ix.p1 <- ix + 1
					if (vp[rl,ix.p1] != is.pool) { # penalise potential of non-pool in v3
						v3[rl,ix.p1] <<- pool.penalty * v3[rl,ix.p1]
						pen.ct <- pen.ct + 1
					}

					# to left of current pool
					ix.m1 <- ix - 1
					if ((ix.m1 > 0) && (vp[rl,ix.m1] != is.pool)) { # penalise potential in v3
						v3[rl,ix.m1] <<- pool.penalty * v3[rl,ix.m1]
						pen.ct <- pen.ct + 1
					}
				}
			}

		} # for rl

		# generate pools, using potential in v3
		for (rl in 1:2) { # right, then left side
	    if (cell.pool.nd < 1) break # safety
			new.pool.ct <- 0
			old.pool.ct <- 0

			for (ix in 1:cell.pool.nd) {
				if (vp[rl,ix] != is.pool) {   # not (yet) a pool
					# should it become one?

					v <- runif (1, min=0, max=pool.risk)
					if (v < v3[rl,ix]) {  # compare with (possibly penalised) value
						vp[rl,ix]  <<- is.pool    # is.pool
						v1[rl,ix]  <<- pool.dpth.st # starting depth 0.1
						new.pool.ct <- new.pool.ct + 1
					}
				} else old.pool.ct <- old.pool.ct + 1
			}
		}

		# smooth base of contiguous pools
	  first.nd <- c (0,0)
	  for (rl in 1:2) { # right then left sides

	    contig.sm <- 0.0
	    contig.ct <- 0
	    ix        <- 0
	    repeat {
	      ix <- ix + 1
	      # seek contiguous pools
	      if (vp[rl,ix] == is.pool) { # start or continue contiguous
	        contig.ct <- contig.ct + 1
	        contig.sm <- contig.sm + v1[rl, ix]
	      } else {                    # not a pool: smooth current block
	        if (first.nd[rl] == 0) first.nd[rl] <- ix - 1
	        if (contig.ct > 1) {
		        contig.dp <- contig.sm / contig.ct
		        for (j in ix-contig.ct:ix) v1[rl,j] <- contig.dp # install smoothed depth
		      }
	        contig.sm <- 0.0
	        contig.ct <- 0
	      }
	      if (ix > pool.cell.lim) break;
	    }
	  }
	  # one pool at centre?
	  dp <- c (0,0)
	  if ((vp[1,1] == is.pool) && (vp[2,1] == is.pool)) { # harmonise centre
	    # extra weight to deeper pool
	    for (rl in 1:2) dp[rl] <- v1[rl,1] # depths of right and left pools
	    wt <- c (0,0)
	    if (dp[1] > dp[2]) wt <- c (1.9, 0.1) else wt <- c (0.1, 1.9)
	    sm <- wt[1] * first.nd[1] * v1[1,1] + wt[2] * first.nd[2] * v1[2,1]
	    mn <- sm / (first.nd[1] + first.nd[2])
	    for (rl in 1:2) {
	      for (n in 1:first.nd[rl]) v1[rl,n] <- mn
			}
	  }
	} # if change

	# now, plot pools
	for (rl in 1:2) {  # right, then left
	  for (ix in 1:cell.pool.nd) {
	    if (cell.pool.nd < 1) break # safety

	    if (vp[rl,ix] == is.pool) { # is.pool
	      x   <- vx[rl,ix]      # centre of pool
	      x.m <- x - half.step
	      x.p <- x + half.step

	      # CCW from SW
	      xx  <- c (x.m, x.p, x.p, x.m, x.m)

	      adj  <- 0.2
	      y.dp <- v1[rl,ix] + adj
	      w    <- x^2 / cur.wide^2
	      if (w <= 1) { # avoid sqrt (-ve)
		      y.hi <- cur.high * sqrt (1 - w)
		      y.lo <- y.hi - y.dp

		      # CCW from SW
		      yy   <- c (y.hi, y.hi, y.lo, y.lo, y.hi) + acro.dpth
		      polygon (x=xx, y=yy, border="blue4", col="blue4", lend=2)
        }

        # deepen pools
	      if (change) v1[rl,ix] <<- v1[rl,ix] + v3[rl,ix] * pool.dpth.inc
	    }
	  }
	} # for rl

  return ()
} # plot.pools



hyperb.tree.tall <- function (age) {
# hyperbola for tree height, mx / 2 at x (age) = tree.mx.death.yr
  b   <- tree.mx.death.yr
  a   <- tree.high.at.death
  high <- (a * age) / (age + b)
  return (high)
} # hyperb.tree.tall

plot.trees <- function () {
# plot trees on periphery of raised bog

  # vt tree age
  # v2 slope.p
  # v3 y when tree started (save)



  if (change) {
    # find limits to tree area (same on right and left sides)
#    b.a.sq   <- (cur.high * cur.high) / cur.wide.sq
    cell.tree.nd <<- floor (cur.wide / x.step)   # current last bog cell
    if (cell.tree.nd <= 3) return () # too few cells
    cell.tree.st <<- -1 # unset
    ix.st         <- cell.pool.nd + 1
    if (ix.st < 0) ix.st <- 0
    ix            <- cell.tree.nd # bog wide
    # now work back towards centre
    repeat {
      if ((ix < ix.st) || (ix < 1))  break
#      x   <- vx[1,ix]
#      w   <- x^2 / cur.wide^2
#      if (is.na(w) || (w > 1.0)) return ()    # skip
#      y.high         <- cur.high * sqrt (1 - (w))
#      slope          <- b.a.sq * x / y.high
      slope          <- vs[ix]
      slope.p        <- h.p.crit / slope # < 1 moving towards tiny
      v2[1,ix]      <<- slope.p # tree age??
      v2[2,ix]      <<- slope.p
      if ((slope < h.t.crit) && (cell.tree.st < 0)) {
        cell.tree.st <- ix
        ix.st        <- ix - tree.cell.overrun
        if (ix.st < cell.pool.nd) ix.st <- cell.pool.nd + 1
        if (ix.st < 0) ix.st <- 0
      }
      ix             <- ix - 1         # next toward bog centre
    } # repeat

    # safety check
    if ((ix.st < 1) || (cell.tree.nd < 1) || (cell.tree.nd < ix.st)) return ()

    # remove trees?
    for (rl in 1:2) {
      for (ix in ix.st:cell.tree.nd) {
        if (vt[rl,ix] > 0) { # is.tree
          r <- runif (1, min=0, max=tree.mx.death.yr)
          # tree dies? (age weighted by shallower slope
          if ((vt[rl,ix] / v2[rl,ix]) > r) {
            vt[rl,ix] <<- 0  # is.tree
            v3[rl,ix] <<- 0  # missing
          }
        }
      }
    }

    # grow, or start a new tree?
    for (rl in 1:2) {
      for (ix in ix.st:cell.tree.nd) {
        if (vt[rl,ix] > 0.0) { # is.tree
          # tree ages
          vt[rl,ix] <<- vt[rl,ix] + tree.grth.fct * t.frm.inc
          # to do (perhaps) modify to allow growth = f (steepness of slope)
        } else { # NOT a tree; start a new one?
          if (runif (1, min=0, max=1) > tree.start.prb) {
            vt[rl,ix]   <<- t.frm.inc # age
            x            <- vx[rl,ix]
            w            <- x^2 / cur.wide.sq
            y.high       <- cur.high * sqrt (1 - w)
            v3[rl,ix]   <<- y.high
            new.trees.ct <<- new.trees.ct + 1
          }
        }
      }
    }

    ix.trees <<- ix.st # carry ouside 'if change'

  } # if change


  # draw trees
  ix.st <- ix.trees
  for (rl in 1:2) {
    for (ix in ix.st:cell.tree.nd) {
      if (vt[rl,ix] > 0) {
        xx <- vx[rl,ix]
        yy <- v3[rl,ix]
        points (x=xx, y=yy, pch=21, cex=0.5)
        width <- 1.8 * x.step
        high   <- hyperb.tree.tall (vt[rl,ix])
        sym.tree (x=xx, y=yy, kind="Z", high=high,
                  trk.prpy=3,
                  cpy.prpy=2, cpy.wyd=width,
                  from="plot.trees")
      }
    }
  }

return ()

} # plot.trees



plot.lagg <- function () {
# trees in lagg stream

  trk ("  plot.lagg")

  # lagg stream, as open topped rectangle
  u  <- cur.wide + x.step
  uu <- cur.wide + stream.cells * x.step
  v  <-  0.0
  vv <- -0.25
  # topleft in U to top right
  xx <- c (u,  u, uu, uu)
  yy <- c (v, vv, vv,  v)
  dir <- 1
  for (i in 1:2) {
    points (x=dir*xx,  y=yy, type="l", col="blue", lwd=5)
    dir <- -dir
  }

  # lagg fen trees
  # settings for 3 sorts of lagg trees of high, cpy.wyd, trk.top
  lagg.tree.par       <- c (1.8, 6.0, 0.4,  1.0, 8.0, 0.2,  2.1, 5.0, 0.4)
  dim (lagg.tree.par) <- c (3, 3)

  xx <- uu
  ix <- floor (xx / x.step)
  for (loc in 1:3) { # 3 deciduous trees, differing shapes
    dir <- -1
    if (change) {
      ht <- tweak (lagg.tree.par[1, loc], 0.5)
      cw <- tweak (lagg.tree.par[2, loc], 0.6)
      tt <- tweak (lagg.tree.par[3, loc], 0.2)
      # save
      v2[1,ix] <<- ht
      v3[1,ix] <<- cw
      v4[1,ix] <<- tt
      v5[1,ix] <<- xx
    } else {
      # restore
      ht <- v2[1,ix]
      cw <- v3[1,ix]
      tt <- v4[1,ix]
      xx <- v5[1,ix]
    }
    for (j in 1:2) { # left and right sides
      sym.tree (x=dir*xx, y=0, kind="E", high=ht,
                y.size=tree.high.fct, x.size=tree.wide.fct,
                trk.top=tt, trk.wyd=0.2, trk.prpy=8,
                cpy.bot=0.2, cpy.wyd=cw, cpy.prpy=2,
                from="plot.lagg")
      dir <- -dir
    } # j in 1:2
    xx <- xx + tweak (2, 0.2) * x.step
    ix <- ix + 1
  } # loc in 1:3
  cur.wide <<- xx
} # plot.lagg


plot.woods <- function () {
# settings for 'woods' around raised bog

  trk ("  plot.woods")

  woods.par     <- c (2.8, 7.0, 0.2)

  if (change) {
    xx      <- cur.wide + x.step
    xx.mx   <- all.wide.mx
    ix.st   <- floor (xx / x.step)
    ix.nd   <- floor (xx.mx / x.step)
    woods.st <<- ix.st
    woods.nd <<- ix.nd
  } else {
    ix.st <- woods.st
    ix.nd <- woods.nd
  }
  # woods around
  for (rl in 1:2) {
    for (ix in ix.st:ix.nd) {
      if (change) {
        if (runif (1, max=1.0) > 0.2) tree.sort <- "T" # 'Triangle' (spruce)
        else tree.sort <- "Z"                          # 'Zig-zag' (pine)
      }
      # tree sizes
      if (change) {
        xx <- vx[rl,ix]
        ht <- tweak (woods.par[1], 0.2)
        cw <- tweak (woods.par[2], 0.2)
        tt <- tweak (woods.par[3], 0.2)
        # store
        v2[rl,ix] <<- ht
        v3[rl,ix] <<- cw
        v4[rl,ix] <<- tt
        v5[rl,ix] <<- xx
        v6[rl,ix] <<- tree.sort
      } else {
        # retrieve
        ht <- v2[rl,ix]
        cw <- v3[rl,ix]
        tt <- v4[rl,ix]
        xx <- v5[rl,ix]
        tree.sort <- v6[rl,ix]
      }
      if (tree.sort == "T") cp <- 4
      else                  cp <- 8
      sym.tree (x=xx, y=0, kind=tree.sort, high=ht,
                y.size=tree.high.fct, x.size=tree.wide.fct,
                trk.top=tt, trk.wyd=0.2, trk.prpy=3,
                cpy.bot=0.1, cpy.wyd=cw, cpy.prpy=cp,
                from="plot.woods")
    } # for ix
  } # for rl

} # plot.woods

run.intro <- function () {
# Introduction explaining how to pause, and colour meanings

  trk ("run.intro")
  # first (or only) part of labels 1:6, 8:11
  sim.sec.txt <- as.character (sim.sec)
  duration    <- paste ("10,000 years in ", sim.sec.txt, " seconds", sep="")
  tx1     <- c (
                "Try <space> or <ctrl-P> to PAUSE",       # 1
                "VERTICAL SECTION OF GROWING RAISED BOG", # 2
                duration,                                 # 3
                "Green:",                                 # 4
                "Brown:",                                 # 5
                "Grey: ",                                 # 6
                "Blue  "                                  # 7
               )
  # second (uncoloured) part of coloured labels 4:7
  tx2     <- c (" ", " ", " ",                            # 1:3
                "living vegetation and acrotelm",         # 4
                "top of catotelm peat",                   # 5
                "history of older surfaces",              # 6
                "pools"                                   # 7
               )
  clr     <- c ("black", "black", "black",
               "green4", "brown4", "grey50", "blue4")
  #               4         5         6         7

  n.lin       <- NROW (tx1)

  # line spacing
  # lin down from top
  lin.tall    <- c (14, 14, 14, 14, 11, 11, 11) # individual line spaces

  lin.below   <- 15 # for trees

  y.val       <- rep (NA, n.lin)
  y.cum       <- 0
  for (i in 1:n.lin) {
    y.cum     <- y.cum + lin.tall[i]
    y.val[i]  <- y.cum
  }

  y.top <- y.cum + lin.below + 1
  y.val <- y.top - y.val # y coordinate
  y.val <- y.val + lin.below

  show.pause <- T
  in.sec     <- intro.sec

  flicker <- 3 # smaller for faster

  # produce frames
  for (fr in 1:intro.frames) {

  # file name and device
	  fr.name <- paste (topog, "_", dk.model, "_", time.name(fr), sep="");
	  png.dev (out.f=fr.name,
	           wide=scrn.wide, high=scrn.high, size=fnt.size)

    # first line flicker control
    if ((fr %% flicker) == 0) show.pause <- !show.pause

    # seconds countdown
    if ((fr %% fps) == 0) in.sec <- in.sec - 1

    # blank plot
    x.lim <- c (0, 10) # 10 is arbitrary, but other things then depend on it
    y.lim <- c (0, y.top)
    plot(1, type="n", axes=F, xlab="", ylab="",
         xlim=x.lim, ylim=y.lim)

    # first line, flickers
    t.s <- as.character (in.sec)
    txp <- paste (t.s, "        ", tx1[1], sep="")
    colr <- "gray20"
    if (show.pause) colr <- "gray40"
    text (x=1, y=y.val[1], label=txp, pos=4,
          lwd=2, col=colr, font=2) # font = bold

    # other lines
    for (j in 2:n.lin) {
      n.y <- y.val[j]
      if (j < 4) {  # lines 2, 3 simple text
        text   (x=1,   y=n.y, label=tx1[j], pos=4, col=clr[j])
      }
      if (j >= 4) { # lines in 2 parts,coloured, plain
        points (x=1,   y=n.y, type="p", pch=15, col=clr[j], cex=2)
        text   (x=1.3, y=n.y, label=tx1[j], pos=4, col=clr[j])
        # first part was coloured; now do 2nd part plain
        text   (x=2.5, y=n.y, label=tx2[j], pos=4, col="black")
      }
    }


    # trees, at bottom
    xx   <- 1
    x.on <- 1.2
    yy   <- 0
    y.hi <- 20.0
    x.wd <- 0.5
    x.sz <- 1 / y.hi
    if (do.trees) {
       sym.tree (x=xx, y=yy, x.size=x.sz, kind="Z", high=y.hi,
                 from="intro/do.trees")
       text (x=xx+x.wd, y=yy, labels="Pine")
       # advance to next
       xx <- xx + x.on
    }
    if (do.lagg) {
       sym.tree (x=xx, y=yy, x.size=0.5*x.sz, kind="E", high=y.hi,
                 from="intro/do.lagg")
       text (x=xx+x.wd, y=yy, labels="Deciduous")
       # advance to next
       xx <- xx + x.on
    }
    if (do.woods) {
       sym.tree (x=xx, y=yy, x.size=x.sz, kind="T", high=y.hi,
                 from="intro/do.woods")
       text (x=xx+x.wd, y=yy, labels="Spruce")
       # advance to next
       xx <- xx + x.on
    }

    dev.off () # done with current device

  } # for fr

  cur.frame <<- fr

} # run.intro


run.sim <- function () {
# simulate and plot

  trk ("run.sim")

  change <<- T

  tot.frames    <<- sim.frames
  if (do.intro) {
    tot.frames <<- tot.frames + intro.frames
    bog.frm.start <<- intro.frames
  } else {bog.frm.start <<- 0}

  freeze         <- tot.frames + 1

  if (do.coda)  tot.frames <<- tot.frames + coda.frames + 1

  prog.ch <- '@' # progress monitor (in '.Rout')

  offs        <- 0
  # the main loop
  for (fr in 1:tot.frames) {

    frm.ct   <<- fr + cur.frame

    if (frm.ct == freeze) {
      change <<- F # freeze state at end of main run
      cat ("\nStarting coda. No more changes after frame ", freeze, NL)
    }

    if ((frm.ct > bog.frm.start) && change) bog.frm.ct <<- bog.frm.ct + 1

    if (change) t.yr <<- fr * t.frm.inc

    # copy last sim.frame to '0000'
    if (frm.ct == freeze + 1) {
      t.name <- time.name (0)
      offs   <- 1
    } else  t.name <- time.name(frm.ct - offs)

    # set graphic device
    fr.name <- paste (topog, "_", dk.model, "_", t.name, sep="");
    png.dev (out.f=fr.name,
            wide=scrn.wide, high=scrn.high, size=fnt.size)

    # skeleton plot
    par(fig=main.crnrs) # set plot space and position
    plot (-1, type="l", xlim=x.lim, ylim=y.lim,
          bty="n", las=1,
          xlab="", ylab="", axes=F)
    x.tck <- c (-800, -600, -400, -200, 0, 200, 400, 600, +800)
    axis (1, at=x.tck, labels=NULL, pos=0.0, xaxs="i")
    y.tck <- c (0, 2, 4, 6, 8)
    axis (2, at=y.tck, labels=NULL)
    mtext (text="Distance / m", side=1, line=1.5, las=1)
    mtext (text="Height / m",   side=2, line=2, las=0)
    # baseline
    points (x=x.lim, y=c(0.0,0.0), type="l",
            lwd=1.5, col="black")

    # calculate height, and save ghosts
    cur.high <<- peat.sim (p.1, a.3, vp.2, va.4, t.yr, dt.6, sav.ghost=T)

    # store current height and time for growth vs time plot
    height.at[fr] <<- cur.high
    time.n[fr]    <<- t.yr

    # set global current width, height for ellipse plotting
    if (topog == "EXos") {
      cur.wide <<- cur.high * asp.x.ov.y # expanding
    } else if (topog == "STos") {
      cur.wide <<- wide.mx               # stationary, full width
    } else {
      cat ("'topog' must be either 'STos' or 'EXos', but is", topog, NL)
      quit (status=77)
    }
    cur.wide.sq <<- cur.wide * cur.wide


    # catotelm top
    if (do.cato)
      plot.ellipse (height=cur.high, offs=0.0, width=cur.wide,
                    colr="brown4", thick=6)

    # acrotelm
    if (do.acro)
      plot.ellipse (height=cur.high, offs=acro.dpth, width=cur.wide,
                    colr="green", thick=acro.lwd)

    # add old surfaces
    if (do.ghosts) plot.ghosts ()

    # slopes, needed by plot.pools and plot.trees
		if (do.pools || do.trees) calc.slopes ()

    # add pools
    if (do.pools) plot.pools ()

    # add trees to bog
    if (do.trees) plot.trees ()

    # add trees to lagg fen
    if (do.lagg) plot.lagg ()

    # add surrounding forest
    if (do.woods) plot.woods ()

    # plot size vs time
    if (do.grth.grph) plot.growth ()

    # frame done
    dev.off () # close this device

    # monitor progress
    cat (prog.ch)
    if (prog.ch == '@') prog.ch <- '='
    else prog.ch <- '@'

  } # for fr

  cat (NL)
} # run.sim

# MAIN (ORCHESTRATE whole process)
ORCHESTRATE <- function () {

  trk ("ORCHESTRATE")

  # start pseudo-random number generator
  set.seed (prng.seed)

  # check settings
  if (do.trees && (x.step >= 10)) {
    cat ("do.trees=T, but x.step", x.step, " must be < 10 m\n")
    quit (status=90)
  }

  fr.offs <- 0

  if (do.intro) run.intro ()
  cur.frame # check

  if (do.sim) {
    # find productivity (p.1) to give high.mx after tlim.5 years
    w           <- uniroot (f=sim.try, lower=1e-7, upper=1e-1, maxiter=500)
    p.1        <<- w$root
    nodk.mx.yr <<- high.mx / p.1              # yr for top if no decay
    cat ("'Productivity' (p.1) calculated as", 1000 * p.1, " mm / yr\n")
    run.sim ()
  }

} # ORCHESTRATE


# MAIN


ORCHESTRATE () # controls what to do



warnings () # sometimes useful/needed during development

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#  You can run the program with 'R CMD BATCH ???.R' where '???' is the stem name  #
#  of the program e.g. 'Bog_sim_V6.R'. The program should run for several minutes #
#  without any indication of progress. It will produce a file ???.Rout (e.g.      #
#  'Bog_sim_V6.Rout' that lists messages. If it fails, an explanation should be   #
#  in this ???.Rout file                                                          #
#                                                                                 #
#  Or you can run the program with 'Rscript ???.R' (e.g. 'Rscript bog_sim_V6.R'). #
#  This mode will produce several bits of information, and will record progress   #
#  on the terminal.                                                               #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# You can run the program with 'R CMD BATCH ???.R' where '???' is the stem name   #
# of the program e.g. 'Bog_sim_V6.R'. The program should run for several minutes  #
# without any indication of progress. It will produce a file ???.Rout (e.g.       #
# 'Bog_sim_V6.Rout' that lists messages. If it fails, an explanation should be    #
# in this ???.Rout file                                                           #
#                                                                                 #
# Or you can run the program with 'Rscript ???.R' (e.g. 'Rscript bog_sim_V6.R').  #
# This mode will produce several bits of information, and will record progress    #
# on the terminal.                                                                #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# To make the video: in a terminal use the command:                               #
# ffmpeg -f image2 -i ??%4d.png -r 25 ##.mp4                                      #
# where '??' is stem name of files, and '##' is chosen output name                #
# (the '%4d' shows that each file name is followed by 4 digits).                  #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

quit (status=0)
