#  **May 15, 2007**

# **Documentation for SPBLOB: A Computer Simulation** 

# **Model of the Joint Population Dynamics for** 

#  **Loblolly Pine and the Southern Pine Beetle** 

# 

# 

by

> John Bishir

James H. Roberds

Brian L. Strom

Xiaohai Wan

# 

# 

#  Publication No.

####### USDA Forest Service

#  Southern Research Station

#### Preface

This publication describes and documents **SPBLOB**, a computer
simulation model for initiation, expansion, and decline of Southern Pine
Beetle (SPB: *Dendroctonus frontalis* Zimm.) infestations in Loblolly
Pine (Pinus taeda L.) plantations from time of stand establishment until
harvest. It consists of four components:

-   **Documentation** for **SPBLOB**, which provides information about
    the code structure and the values for the parameters contained in
    the **SPBLOB** code

```{=html}
<!-- -->
```
-   Two **User Guides** that explain the choices available in setting up
    stand simulations, together with interpretation of simulation output

-   A **Glossary** that lists and defines the constants, variables and
    parameters used in the simulation code

-   **References** upon which **SPBLOB** model structure and the major
    code parameter values are based

Bishir et al. (2007) give a description and explanation of the
**SPBLOB** model. The **Source Code** for **SPBLOB,** and two
**Executable Versions** - one for specifying and running single
simulations, the other for multiple independent simulations - can be
downloaded from **www.srs.fs.usda.gov/4501**.

While this code has been extensively tested and checked against results
of field studies reported in the literature, undetected errors have no
doubt crept in. We will be most grateful to users who inform us of
these.

**Legalities**

The compilers, random number generators and command files used in
connection with the **SPBLOB** code have been downloaded from sites in
the public domain. Citations and credits appear as appropriate in
**Section V** of this Documentation. More detailed information,
including Disclaimers and Legal Requirements associated with these uses,
is contained in the **Appendix**.

**Disclaimers Associated with SPBLOB Itself**

We consider the **SPBLOB** model and its accompanying simulation code
primarily as research tools, rather than accurate predictors of stand
timber volume or degree of beetle activity in *particular* stand
scenarios. The reason is that fundamental beetle activities - egg
production, development in juvenile stages, movement, and mortality -
are controlled largely by temperature (Payne, 1980; Coulson, 1980;
Wagner et al. 1983/1984; Wagner et al. 1984; Stephen and Lih 1985;
Ungerer et al.1999) which, in turn, is subject to unpredictable chance
fluctuations. In addition, field observations indicate that in some
stands beetles can appear in considerable numbers (e.g., via migration
from outside) by age 10 or before, while others remain free of beetle
infestations throughout their existence. Thus, model comparisons of
management strategies (e.g., clonal plantings vs those generated from
seedlings or, in clonally generated stands, block plantings vs random
mix of ramets) should be based, not on single stand simulations, but
rather on multiple, independent replications of any potential
strategies.

No model is perfect or complete, in the sense that no model can include
all the subtle details and interactions that exist in a real system. Our
model focuses on the joint dynamics of southern pine beetles and
loblolly pines, in particular on tree growth and oleoresin flow rates,
the former because it allows prediction of timber volume (Burkhart et
al. 1987) and the latter because \"the oleoresin system in pines is the
primary defensive mechanism against attack by bark beetles\" (Blanche et
al. 1983). In doing so, it does not include beetle competitors or
predators, nor does it include effects on trees of important pathogens
introduced by attacking beetles. And it omits a variety of mechanisms,
such as hypersensitive response (Berryman, 1972) and
growth/differentiation balance (Lorio, 1986; Lorio et al. 1990) known to
affect beetle populations.

In addition, estimates of important model parameters vary in the
literature, and we have had to choose among them. Further, to our
knowledge some processes and parameters included in the model and its
simulation code have not been studied quantitatively. For these, we have
assigned parameter values that result in code output that approximates
observed stand outcomes, in particular, the qualitative behavior of
beetles and the tree mortality they cause.

Finally, while this code has been extensively tested and checked against
results of field studies reported in the literature, neither the USDA
Forest Service nor the authors accept responsibility for errors or
inaccuracies that arise. In particular, we provide absolutely no
warranty of any kind, either expressed or implied, including but not
limited to, warranties of merchantability or fitness for a particular
purpose.

##### Authors 

John Bishir is Professor Emeritus in the Department of Mathematics,
North Carolina State University, Raleigh, NC 27695-8205, USA
([bishir@math.ncsu.edu]{.underline}). Xiaohai (Robert) Wan received his
PhD in Biomathematics at the same university
([xwan@unity.ncsu.edu)]{.underline} and is currently located at Capital
One Services, 15000 Capital One Drive, Richmond, VA 23238
([robertwan94@yahoo.com)]{.underline}. James H. Roberds is Research
Geneticist, Southern Institute of Forest Genetics, Southern Research
Station, USDA Forest Service, Saucier, MS 39574, USA
([jroberds@fs.fed.us]{.underline}). Brian Strom is Research
Entomologist, Research Unit on Insects, Diseases and Invasive Plants,
Southern Research Station, USDA Forest Service, Pineville, LA 71360, USA
([brianstrom@fs.fed.us)]{.underline}

##### Acknowledgments

Development of this model was supported by grants provided by the
Southern Institute of Forest Genetics, Southern Research Station, USDA
Forest Service, Saucier, Mississippi, and by the Southern Research
Station, IDIP, USDA Forest Service, Pineville, Louisiana.

Our special thanks go to Professor **Harold Burkhart**, Department of
Forestry, Virginia Polytechnic Institute and State University (VPI), and
his colleagues at VPI and the Westvaco Corporation, Summerville, South
Carolina, who generously allowed us to include parts of **PTAEDA2**,
their simulation code for competitive tree growth and mortality
(**Daniels and Burkhart, 1975; Burkhart et al. 1987**).

##### Abstract

# 

**SPBLOB** is a computer simulation model for the interaction between
loblolly pine (Pinus taeda L.), the economically most important forest
crop in the United States, and the southern pine beetle (SPB:
*Dendroctonus frontalis* Zimm.), the major insect pest for this species.
The model simulates loblolly pine stands from time of planting until
harvest. It mimics day-to-day changes in SPB populations, and the
associated tree mortality caused by these bark beetles. In addition, it
provides yearly updates of tree mortality due to competition and of
growth of the surviving trees.

Model output relates to both trees and beetles. Simulations provide
daily and yearly fluctuations in the size, composition and spread of
beetle populations within a stand. Stand types -- e.g. plantations
derived from seedlings or from clonally generated materials -- can be
compared in terms of average total merchantable volume of timber per
hectare at time of harvest and/or by average proportions of surviving
trees. Average stand yields in the presence of beetles can also be
compared with those when beetles are absent, thus providing a basis for
estimating the real cost of beetle damage in loblolly stands.

In this **Documentation** we briefly describe the model, then outline
the concepts and structure upon which the **SPBLOB** simulation source
code is based. The included **User Guides** specify the input
information regarding a loblolly stand and its environment required to
run the Executable **Simulation Codes**. These **Guides** also describe
the extensive variety of information about simulated tree and beetle
populations that can be retrieved during simulations.

> **Contents**

##### I. INTRODUCTION 1

II\. THE MODEL AND ITS SIMULATION CODE 2

A. Introduction 2

B. Overview 3

C. Plantation Establishment 5

C.1 Setting the Stand Parameters -- Subroutine MAIN 5

C.2 Choosing Tree Materials -- Subroutines NOCLONE, CLONE

and DEPLOY 6

C.2.1 Non-clonal Stands - Subroutine NOCLONE 6

> C.2.2 Stands Using Clonal Materials - Subroutines DEPLOY 6
>
> and CLONE

#####  D. Simulating the Planted Stand 7

D.1 The Tree Juvenile Period -- Subroutine JUVEN 7

#####  D.2 Setting Beetle Population Parameters - Subroutine PARBTL 7

#####  E. Age 9 until Harvest -- SIMUL and its Subroutines 8

#####  E.1 Tree Growth and Mortality Related to Intrastand Competition 8 E.1.1 Merchantable Timber Volume -- Subroutine OUTPUT 8

> E.1.2 Tree Competition Indices -- Subroutine COMP2 9
>
> E.1.3 Tree Growth and Mortality -- Subroutine GROW2 9
>
> E.2 Setting Beetle-related Stand and Tree Parameters --
>
> Subroutine BEGYR 9
>
> E.3 Daily Simulation of Lightning Strikes -- Subroutine
>
> LIGHTNING 9
>
> E.4 Daily Updates of Beetle Populations on Each Tree and of
>
> the Tree Mortality they Cause -- Subroutines TEMP
>
> through DISPERSE 10
>
> E.4.1 Generating Daily Temperatures - Subroutine TEMP 12 E.4.2 Today's
> Oviposition Rate - Subroutine EGG 12
>
> E.5 (Continuation of Section E.4) Command Central for the
>
> Simulation Code -- Subroutine UPDATE 12

#####  E.5.1 Imposing Winter Mortality - Subroutine COLD 13 E.5.2 The Spread of SPB Pheromones - Subroutine DIFFUSE 13

#####  E.5.3 Beetle Penetration of Tree Bark - Subroutine ATTACK 13

E.5.4 Tree Death from Beetle Penetrations 13

#####  E.5.5 Summing up the Day's Activities -- Subroutine ENDOFDAY 14

#####  E.5.6 Moving the Departing Adult Beetles - Subroutine DISPERSE 14 

III\. USER GUIDE TO IMPLEMENTATION OF THE INTERACTIVE

WINDOWS-BASED SPBLOB CODE 16

A. Introduction and Overview 16

B. Downloading and Installing SPBLOB on Your Machine 16

C. The Interactive, Input Phase of SPBLOB 20

C.1. The Input Phase 20

C1.1 Number of Clones When Blocks Are Not Used 23

C.1.2 More Information Related to Blocks 23 C.1.3 Clonal Parameter
Values 24

C.1.4 Wrapup When Clones Are Being Used 25

C.1.5 No Clones 26

C.2. Simulation Output 26

C.2.1 What Output Information Is Saved to Files 28

C.2.1.1. The tree-information file 28

C.2.1.2. The stand-information file 30

IV\. GUIDE TO USING THE EXECUTABLE MULTIPLE-SIMULATION

> VERSION OF SPBLOB 32

#####  A. How to Run Batch Simulation Jobs 32

> A.1. The file configs.txt 32

A.2. The file input.txt 35

A.3. The file run.txt 36

> A.4. The file tree.exe 37

B. The Simulations Themselves - What Output Information Is Saved

> to Files 40

B.1. The tree-information files 40 B.2. The stand-information files 42

##### V. THE SPBLOB Model SOURCE CODE AND ITS Parameters, WITH 

#####  EMPHASIS ON PARAMETER Values and THEIR Effects 

##### ON MODEL OUTPUT 44

A. Introduction 44

A.1. The SPBLOB Source Code 44

A.1.1 Modifying and Recompiling the Source Code 44

#####  A.2. SPBLOB Parameters: Values, Sources, and Effects on Model 

#####  Output 46

##### B. Parameters Associated With Stand Planting 46

#####  B.1. The Parameter g0 47

#####  B.2. The Parameter RMAX 47

##### B.3. Assignment of Phenotypic Values to Trees 47

B.4. Edge Effects 48

##### C. Parameters Associated with Beetle Population Dynamics 48

##### C.1. Computation of SPB Vital Rates - Subroutine PARBTL 49

> C.2. Beetle Immigration and Emigration 49
>
> C.2.1. Immigration 49
>
> C.2.2. Emigration 50

D. Simulation of Temperatures 50

##### E. Frequency of Lightning Strikes 51

F. Beetle Mortality Due to Extremely Cold Temperatures 52

G. Juvenile Development and Mortality over the Course of a Day 52

H. The Process of Beetle Attack on Trees - Subroutine ATTACK 53

##### J. Tree Response to Beetle Attack - Subroutine ENDOFDAY 53

##### J.1. Resin Loss from Flow Through Attack Holes 54

> J.2. Tree Mortality Due to Beetle Penetrations 54

##### K. Activities of Parent Beetles -Subroutine ENDOFDAY 55

##### K.1. Proportion of Her Total Eggs That a Female Parent Deposits 

##### in a Day 55

> K.2. Total Number (newegg) of Eggs Laid in Each Tree 55
>
> K.3. Daily Parent Survival and Reemergence 56

##### 

##### L. Pheromone Production, Concentrations and Diffusion 57

> L.1. Production - Subroutines UPDATE and ENDOFDAY 58
>
> L.2. Diffusion and Resulting Concentrations - Subroutine DIFFUSE 58

M. Beetle Dispersal - Subroutine DISPERSE 59

N. Tree Decay and Fall - Subroutine SIMUL 59

P. Merchantable Volume - Subroutine OUTPUT 59

##### VI. Glossary for the SPBLOB Source Code 60

VII\. References 76

VIII\. Appendix - Legalities 81

A. Introduction 81

B. Web Sources 81

C. Legalities 82

C.1 f2kcli 82

C.2 SYBASE, INC. 82

C.3 RANDLIB 82

C.3.1 From the RANDLIB90 README FILE 82

C.3.2 The RANDLIB90 LEGALITIES FILE 83

##### I. Introduction

This report documents and provides User Guides for the computer
simulation model **SPBLOB** (Bishir et al. 2007), a model for the joint
population dynamics between loblolly pine and the southern pine beetle
over the entire life of a stand, from planting until harvest. Loblolly
pine (*Pinus taeda* L.) is the leading commercial timber species in the
United States, growing on more than 13.4 million hectares (Schultz
1997). Southern pine beetle (SPB: *Dendroctonus frontalis* Zimm.) is the
most destructive insect pest of pine forests in the southeastern United
States. Although SPB attacks all pine species within its range (Payne
1980), *Pinus taeda* is its major host**.**

The **SPBLOB** model tracks day to day changes in beetle populations on
and in each tree in a stand, as well as the tree mortality caused by
beetle attacks. In addition, it provides yearly updates of tree growth
and mortality associated with intra-stand competition. Other models and
simulation packages for projecting SPB dynamics exist (Stephen and Lih
1985; Lih et al. 1995; Coulson et al. 1989). However, these are
primarily 'single season' models, focusing on \'spot growth\', the rapid
increase and spread of SPB that periodically occurs when beetle
populations reach critical mass. Beetle populations wax and wane in our
model as well, and spots occur, these being driven by the fundamental,
underlying time-varying interactions of temperature, numbers of beetles
at different levels of maturity, and individual tree locations and their
susceptibility to beetle attack.

The **Source Code** for model simulation can be obtained by going to
**www.srs.fs.usda.gov/4501**. Two **Executable Versions** of this code
can also be downloaded from this site. The first is an **interactive**
code that enables a user to specify a desired stand setup, and to run
individual simulations of that stand. The second runs **multiple
independent simulations** of any of 23 different stand scenarios, by
entering a single command. Sections III and IV below serve as respective
**User Guides** for these executable codes.

Section II gives a brief description of the overall structure of the
model code, including the purpose and output of each of its subroutines.
Following the **User Guides**, Sections V, VI and VII contain,
respectively, detailed information concerning the values and literature
sources for the most important **parameters** in the model; a
**Glossary** of principal variables and parameters that appear in the
source code; and a list of **References** cited in this
**Documentation**.

**II. The Model and Its Simulation Code**

**A. Introduction**

This Section provides an overview of the **SPBLOB** source code,
including a descriptive summary of each of the major subroutines.
Further details appear in Section V, which explains the principal
computations in the code, along with the major parameters, their values
and their sources in the literature.

A broad description of the purpose, structure and outputs of the
**SPBLOB** model appears in Bishir et al. (2007). Briefly, **SPBLOB** is
a computer simulation model of the joint population dynamics between
loblolly pine (*Pinus taeda* L.) and the southern pine beetle (SPB:
*Dendroctonus frontalis* Zimm.). Our original impetus for developing
this model was to evaluate deployment patterns of clonally generated
tree materials. Our goal was to determine the pattern -- individual
clonal blocks or a random mixture of trees from all clones -- that
produces, on average, the greatest amount of merchantable timber from a
stand (see, e.g., Muhs 1993; Roberds et al. 1990; Roberds and Bishir
1997; Bishir and Roberds 1999). The nature of this question requires
knowledge of individual tree attributes and their effects on beetle
populations. Thus, we settled on a model that is discrete in

-   **space** -- Each tree is a separate entity, its interactions with
    its neighbors being affected by distances to them and by its size
    relative to theirs

-   **population size** -- Trees and beetles can be counted, and
    individuals of each species react to their own immediate
    surroundings

-   **time** -- There are yearly updates of the growth of surviving
    stems and of tree mortality due to competition. During each year,
    the model tracks day-to-day changes in the beetle numbers on each
    tree, and the resulting tree mortality they cause.

The flexibility of this structure allows the model code to simulate
joint loblolly/SPB population dynamics over the entire life of a stand,
from planting until harvest. In this generality, the model is capable of
application to any loblolly stand - ranging from natural stands to
commercial plantings - together with any level of beetle populations.

This having been said, we consider the model and its accompanying
simulation code primarily as research tools, rather than accurate
predictors of stand timber volume or the degree of beetle activity in
*particular* stand scenarios. The reason is that primary beetle
activities - egg production, development in juvenile stages, movement,
and mortality - are controlled largely by temperature (Payne 1980;
Coulson 1980; Wagner et al. 1983/1984; Wagner et al. 1984; Stephen and
Lih 1985; Ungerer et al. 1999) which, in turn, is subject to
unpredictable chance fluctuations. In addition, field observations
indicate that in some stands beetles can appear in considerable numbers
(e.g., via migration from outside) by age 10 or before, while others
remain free of beetle infestations throughout their existence. Thus,
model comparisons of management strategies (e.g., clonal plantings vs
those generated from seedlings or, in clonally generated stands, block
plantings vs random mix of ramets) should be based, not on single stand
simulations, but rather on multiple, independent replications of any
proposed strategy.

Model computations relating to SPB focus almost entirely on females
because they alone oviposit and they, not males, select host trees. In
the model code, males play only one significant role, the production of
repelling pheromone that leads to termination of attack (Section E.5.5
below). Otherwise, all beetle population categories and beetle
activities described below pertain solely to the female portion of the
population.

For convenience, the description that follows concerns commercial
loblolly pine plantations; minor modifications are required for natural
stands. These simulated stands are subject to beetle infestation. If
desired, however, a user can specify a 'no-beetle' option. In this case
the code is a simplified version of the tree growth simulator
**PTAEDA2** developed by Daniels and Burkhart (1975) and Burkhart et al.
(1987).

**B. Overview**

The source code for **SPBLOB** is written in Fortran, a Windows-based
version of which can be downloaded from **www.srs.fs.usda.gov/4501**.
Executable versions of this code are also available at this site.

Each stand simulation consists of five stages:

-   Choosing tree materials and planting the stand at the desired
    density

-   Tree growth and survival through an 8 year juvenile period in which
    no beetles are present

-   Subsequent daily updates of the numbers of beetles on each tree, and
    the tree mortality they cause

-   Yearly updates of tree growth and mortality related to tree
    competition

-   Timber harvest

**Figure 1** shows the relationships among the source code subroutines.
Initially, program flow is governed by subroutine **MAIN** which, in
subroutines **NOCLONE** through **CORRELATE** (read left to right in the
top level of the figure), translates user specifications into a
particular stand to be simulated. Subroutine **JUVEN** then projects
stand growth and survival through an 8-year juvenile period during which
no beetles are present, after which subroutine **PARBTL** assigns the
parameters that govern beetle vital rates. At this point, control passes
to subroutine **SIMUL**, which directs traffic for the remainder of the
simulation (2^nd^ and 3^rd^ levels of the figure).

**Figure 1.** Flow diagram for source code subroutines.

Under **SIMUL**, there are:

-   yearly updates by routines **OUTPUT, COMP2** and **GROW2** of tree
    growth and survival related to intrastand competition

-   yearly reinitialization of tree resin flow rates in subroutine
    **BEGYR**

```{=html}
<!-- -->
```
-   lightning strikes simulated daily by subroutine **LIGHTNING**

-   daily changes in beetle populations on individual trees due to
    beetle mortality, oviposition,and juvenile maturation, computed in
    subroutines **TEMP** through **DISPERSE**

-   yearly end-of-season summaries of stand status contained in
    subroutine **STANDINFO**

```{=html}
<!-- -->
```
-   merchantable timber at harvest approximated by subroutine **OUTPUT**

-   summary information about each planted tree recorded at harvest time
    by subroutine **TREEINFO**

In the remainder of this section II we focus on the *structure* of the
code, including the purposes and outputs of the individual subroutines.
Section V below can be consulted for information about the principal
code parameters, their numerical values and their literature sources.

**C. Plantation Establishment**

Planting occurs in stages that correspond to major management decisions
underlying stand establishment

-   assignment of overall stand parameters

-   choice of materials to be planted and the pattern in which to plant
    them

-   planting of the chosen materials

**C.1 Setting the Stand Parameters -- Subroutine MAIN**

This subroutine allows users to specify (Section III below) values of
individual parameters or to choose an entire parameter group from a
prescribed list of parameters (Section IV) that among them determine
stand geometry and harvest age, beetle vital rates, and relevant
environmental quantities. Major parameters include

-   TMEAN = yearly (long term) mean average temperature (degrees C)

```{=html}
<!-- -->
```
-   nyears = harvest age

-   nc = number of columns in the stand

-   nr = number of rows in the stand

-   px = spacing (meters) between columns

-   py = spacing (meters) between rows

-   g0MEAN = average total height of the trees in the stand at age 25

-   numcln = number of clones, if any, to be used (numcln = 0 indicates
    a

> non-clonal stand established from seedlings using selected tree

materials)

In addition, when clonal blocks are used:

-   rspb = number of rows of trees per block

-   cspb = number of columns of trees per block

-   bspr = number of rows of blocks in the stand

-   bspc = number of columns of blocks in the stand

**C.2 Choosing Tree Materials - Subroutines NOCLONE, CLONE and DEPLOY**

In the model, each 'tree' is a pair of fundamental parameter values
assigned by the user or chosen randomly from default distributions
specified in the code. These parameters are

-   g0 = \'expected\' height at age 25 (roughly, an individual tree site
    index)

-   RMAX = the maximum rate at which resin flows from holes in the tree
    bark

that are produced by attacking beetles; RMAX decreases as

greater numbers of beetle penetrations occur

Unless a user specifies otherwise, RMAX follows a gamma distribution, as
described in Roberds et al. (2003), while g0 is normally distributed. In
the case of clonally generated materials, specifically selected for
rapid growth, g0 values are chosen from the upper tail of a distribution
having a higher mean than the distributions used for non-clonal stands,
but RMAX values are randomly generated from the same gamma distribution,
whether or not clones are used. Regardless of how these \'genetic\'
growth and resin rates are assigned, \'site variations\' are added when
trees are planted. Details are in Section V.B.

\[**Note** - Although laboratory studies suggest that resin compounds
are toxic to southern pine beetles (Coyne and Lott 1976)**,** most field
studies show little or no effective resin toxicity (F. Hain, personal
communication). Hence, in the current version of **SPBLOB** we assign a
default value of 0.0 to the toxicity parameter included in the code.
However, unrestricted user choices, including non-zero toxicity values,
can be introduced via changes in the source code.\]

**C.2.1 Non-clonal Stands - Subroutine NOCLONE**

If a non-clonal stand is indicated (numcln = 0) subroutine **NOCLONE**
is called to set the numbers of rows and columns in the stand, after
which each individual tree's (g0,RMAX) pair is randomly chosen from the
distributions described above**.**

**C.2.2 Stands Using Clonal Materials - Subroutines DEPLOY and CLONE**

In clonal plantings all trees of any particular clone have the **same**
pair of g0 and RMAX values. Thus, if numcln \> 0, so clones are to be
used, we generate only a limited number (numcln) of pairs, the default
procedure being to chose randomly and independently from distributions
appropriate to clonal materials. However, users may specify their own
pairs if desired.

Once clones are selected the question becomes one of deployment. If
clones are to be planted in blocks, each parameter pair is assigned to
trees in the same number of blocks. (There are some restrictions
inherent in such plantings, as indicated in Section III.) The general
patterns included in the code are

-   **RM** Randomly mixed clonal ramets - There are clones, and the

> *previously chosen* parameter pairs are sampled randomly, with
>
> replacement, until parameters are assigned to all trees in the stand

-   **CB** Clonal Blocks - There are clones, and trees are planted in a
    mosaic

> of clonal blocks, each block containing trees of a single clone

##### D. Simulating the Planted Stand

Simulation of a stand consists of four stages:

-   tree survival and growth through an 8 year juvenile period

-   from age 9 until harvest, yearly updates of tree growth and
    mortality related to

> competition

-   also beginning in year 9 - daily updates of numbers of SPB on each
    tree and of

> the tree mortality they cause (as indicated in Figure 1, some
> particularly
>
> frequent beetle activities are updated 4 times per day)

-   timber harvest

##### D.1 The Tree Juvenile Period -- Subroutine JUVEN

Because beetles rarely injure young stands, the first 8 years are
treated as a single juvenile period in which there is no beetle activity
and through which tree growth and mortality are projected in one step
using subroutine **JUVEN**, (adapted from **PTAEDA2** with permission:
see Daniels and Burkhart 1975; Burkhart et al. 1987). This routine
randomly assigns appropriate tree mortality during the period, then,
based on the previously assigned growth rates, and competition via
**PTAEDA2**, computes each surviving tree's height, crown ratio and dbh.

##### 

**D.2 Setting Beetle Population Parameters -- Subroutine PARBTL**

If a simulation involves beetles, migration into and out of the stand
occurs during randomly selected years from year 9 until stand harvest.
In addition, there are beetle population changes related to

-   survival

```{=html}
<!-- -->
```
-   movement

```{=html}
<!-- -->
```
-   oviposition

-   juvenile maturation

that occur within the stand.

Rates of these SPB processes depend principally on temperature (Payne,
1980; Coulson 1980; Wagner et al. 1983/1984; Wagner et al. 1984; Stephen
and Lih 1985; Ungerer et al.1999). At the beginning of stand year 9,
subroutine **PARBTL** applies temperature dependent formulas for SPB
vital rates drawn from Ungerer et al. (1999) and Wagner et al.
(1983/1984) to compute daily rates of oviposition and juvenile
maturation for each temperature between 1^o^C and 45^o^C. These rates
are stored in an array and used throughout the remainder of the
simulation.

In addition to temperature, rates of beetle survival and movement depend
on resin flows and pheromone levels associated with individual trees.
Hence, these rates are recomputed daily as population updates proceed
(See Section E.5).

**E. Age 9 Until Harvest -- SIMUL and its Subroutines**

#####  Control now passes to subroutine SIMUL (Figure 1). The first task is to call subroutines OUTPUT and STANDINFO to summarize overall stand status (number of surviving trees, total merchantable timber volume) as of the end of the 8-year juvenile period. In each stand year thereafter, we begin by assessing tree growth and mortality related to intrastand tree competition, then turn to the complex processes that take place daily within the beetle populations on each tree.

##### E.1 Tree Growth and Mortality Related to Intrastand Competition

At the beginning of each year from stand age 9 until harvest, total
merchantable timber in the stand is computed using **PTAEDA2**
subroutine **OUTPUT** (Daniels and Burkhar, 1975); Burkhart et al.
1987). Subroutines **COMP2** and **GROW2**, also adapted from
**PTAEDA2**, then project tree growth and mortality due to competition
in that year (Figure 1).

**E.1.1 Merchantable Timber Volume -- Subroutine OUTPUT**

This subroutine computes an estimate of the total merchantable timber
volume of all live trees in the stand, based on each tree's height and
dbh. The original code (see Daniels and Burkhart (1975); Burkhart et al.
1987) offers a variety of options for computing volume. In **OUTPUT**,
the volume computed is total 'inside bark' volume of trees whose
diameters exceed 4.5 inches (in **PTAEDA** notation, dbh \> 5.0). For
this computation, \'volume\' includes only that part of the tree bole
between stump height and the upper point at which the diameter = 4".

**E.1.2 Tree Competition Indices -- Subroutine COMP2**

This subroutine computes the 'competition index' of each tree, T, as

cip(T) = Σ \[diam(A)/diam(T)\] /dist(A,T),

where diam(A) and diam(T) are the respective diameters (in inches) of
trees A and T, and dist(A,T) is the distance (in feet) between the two
trees. The sum extends over all trees A for which dist(A,T), in feet, is
numerically less than 2.75\*diam(A), in inches.

**E.1.3 Tree Growth and Mortality -- Subroutine GROW2**

Here, the code computes individual tree mortality and growth associated
with intrastand tree competition in each year. For each particular tree,
T, the *probability* that T survives tree competition during a year is

plive = \[1.028\*(cr\*\*0.038)\]\*exp\[-0.0023\*(cip(T)\*\*2.65)\]

where cr is the crown ratio of tree T (vertical length of the crown as a
fraction of the total height of the tree), cip(T) is as computed above,
and exp denotes the exponential function. In the simulation, tree T
survives the year if a uniform random deviate drawn from the interval
0.0 to 1.0 is smaller than plive.

Tree growth is computed as a function of height and diameter increments
(see the **SPBLOB** source code or Burkhart et al. 1987 for details).
These increments are then used together with stand age to determine new
values for tree height, diameter, crown ratio, crown length, and bole
area.

**E.2 Setting Beetle-related Stand and Tree Parameters -- Subroutine
BEGYR**

If the simulated stand contains beetles, at the beginning of each year
subroutine **BEGYR** reinitializes each tree's maximum resin flow rate,
RMAX, to account for tree age. Since the current (2006) literature on
resin flow is inconclusive, the present version of the code holds RMAX
constant throughout a tree's life. However, any desired sequence of
yearly changes can easily be introduced through changes in the source
code.

**E.3 Daily Simulation of Lightning Strikes - Subroutine LIGHTNING**

As many as 30 to 45 percent of all SPB caused 'spots' of tree death in
loblolly stands originate in lightning-struck trees (Hodges and Pickard
1971; Coulson et al. 1983). Thus, simulation of each day's activity
begins with random generation of a Poisson distributed number and
location of lightning strikes occurring that day. The daily mean number
of strikes depends on both the day of the year and the geographical
location of the stand, and can be user-assigned by changing one line in
the source code. Section V.E gives the details.

**E.4 Daily Updates of Beetle Populations on Each Tree and of the Tree**

**Mortality They Cause -- Subroutines TEMP Through DISPERSE**

As noted above, in the model, SPB can enter the stand as early as stand
age 9. Because beetles arriving on a tree decide quickly whether to
attack or move on (Bunt et al., 1980), beetle populations on individual
trees are updated as many as four times per day.

Model recursions track only the female portion of the population, the
male population in each tree being assumed equal to the number of female
parents already in the tree at the beginning of a day (see Bishir et al.
2007 and references therein). Repelling male pheromone levels are
computed from this number, while attracting female pheromone production
depends on the number of successful bark penetrations by females that
day.

A model update begins with all beetles in the stand located on
particular trees, those on each tree being grouped into three
categories: Landed, Parent, and Juvenile. *Landed* beetles are on the
tree exterior and engage in such activities as walking more or less
aimlessly, active searching, investigation of existing entrance holes,
and fighting with other SPB (Bunt et al. 1980; Bishir et al. 2004).
*Parents* are adults that entered the tree previously and are engaged in
gallery construction and oviposition. *Juveniles* are immature beetles,
including eggs and teneral adults, inside the tree bark.

Beetle mortality is assessed in every update period, while other
activities, such as oviposition and juvenile maturation, occur only when
simulated temperatures exceed 7^o^C (Thatcher 1967; Gagne et al. 1982;
Wagner et al 1983/1984; Ungerer et al. 1999).

Figure 2 shows how category populations in and on a particular tree
change during an update. All beetle categories suffer natural mortality,
for which no arrows are shown. In the figure, landed beetles leave or
attack. Those that do not attack enter a temporary group that will
disperse from the tree at the end of the update period. Attackers either
die, from resin flow or predation, or successfully enter the tree to
become parents. Parent beetles oviposit and, when their egg production
in a tree is completed, re-emerge to join the dispersing group.
Juveniles gain maturity and those that reach full adulthood emerge to
join the group of departing beetles.

Other arrows correspond to beetles that leave or enter the stand, or
those that disperse within the stand at the end of the update. For each
individual tree, the new beetles that arrive, either via dispersal from
other trees or by immigration from outside the stand, constitute the
*landed* beetles with which the next update begins.

Because the timing of the attack process varies from tree to tree, on
any particular day, beetles at all stages of development usually exist
in the stand. Thus we cannot speak meaningfully of beetle
\'generations,\' and must treat each day\'s new eggs as a

Figure 2. Flow of beetles among behavioral categories during daily
updates in **SPBLOB**. Solid arrows indicate movement of individual
beetles and the

dashed arrow represents contribution of parents to progeny.

> ![](media/image1.wmf)

separate juvenile cohort (Coulson et al. 1977; Sharpe et al. 1977; Curry
et al. 1978; Coulson, 1980). Members of each cohort have their own level
of maturity, determined by temperatures experienced since they entered
the population. On days with sufficiently high temperatures, each
cohort's maturity level increases. Those cohorts attaining full maturity
emerge to join the dispersers, while new eggs are placed in a cohort of
their own.

#####  

##### E.4.1 Generating Daily Temperatures - Subroutine TEMP

Because of the fundamental role of temperature, each simulated day
begins with subroutine **TEMP**, which generates the day\'s high, low
and average temperatures. We assume long-term daily average temperatures
follow a sine curve, with yearly average and amplitude appropriate to
the geographical location of the stand. Successive daily temperatures
are correlated by setting the expected average temperature for a day
equal to the average of the observed mean of the previous day and the
long-term average for the day's date. A day's actual high, low and
average temperatures are then randomly obtained from a normal
distribution centered at this expected (mean) value. (Details can be
found in the source code and in Baskerville and Emin, 1969. Also see
Section V.D, below.)

##### 

##### E.4.2 Today's Oviposition Rate - Subroutine EGG

This routine computes eggpro, the expected *proportion* of her total egg
complement that each surviving female parent will deposit on a
particular day (Coulson et al. 1978; Gagne et al. 1982). This quantity
is determined as a function of the day's temperatures that were
generated by subroutine **TEMP**.

##### 

##### E.5 (Continuation of Section E.4) Command Central for the Simulation Code - 

##### Subroutine UPDATE

Our model differs from others in the literature in the detail to which
it simulates, on a day-by-day, tree-by-tree basis, the fundamental
activities -- mortality, oviposition, juvenile maturation -- that govern
beetle population dynamics and that cause tree damage and death during
the attack process that leads to beetle reproduction. Subroutine
**UPDATE** represents the core of this simulation activity. It calls the
subroutines -- **COLD, DIFFUSE, ATTACK, ENDOFDAY** and **DISPERSE --**
which compute the four daily updates of the beetle populations on each
tree, and the tree damage these populations cause.

##### E.5.1 Imposing Winter Mortality - Subroutine COLD

Prior to a day\'s updates, subroutine **COLD** computes beetle
mortality, if any, caused by overnight low temperatures. Beetle
populations can plummet during very cold temperatures, the pupal stage
being most severely affected. Temperature dependent mortality rates are
modeled using normal distributions (Ungerer et al. 1999). Details are in
Section V.F below.

##### E.5.2 The Spread of SPB Pheromones - Subroutine DIFFUSE

Simulation experiments suggest that occurrence of 'spots' typically
observed in loblolly pine stands requires two phenomena:

-   rapid diffusion of pheromones through the stand (in our simulations,
    'rapid' means 'by the next stand update')

-   for most of the beetles departing from a tree, movement in the same
    direction within the stand

The second of these requirements is achieved in subroutine **DISPERSE**
below. Prior to that, subroutine **DIFFUSE** computes the day's positive
and negative pheromone concentrations surrounding each tree. These
pheromone plumes then spread throughout the stand. Based on observations
obtained with atmospheric tracer gas (Thistle et al. 2004), if pheromone
production at a tree (T) results in concentration there at level C, this
leads in the model to a concentration of 0.9C in the first ring of trees
around T, a level of (0.9)^2^C in the second ring, etc. The total
concentration near any tree, A, is the sum (superposition) of the levels
induced by beetles in each of the trees in the stand. Details appear in
the file **spblob_params_phero_diffusion.doc** contained in the
**Supporting Information** folder (see the introduction to **Section
V.A** below).

##### E.5.3 Beetle Penetration of Tree Bark - Subroutine ATTACK

Because beetles attack only live trees, all those that land by chance on
dead trees during the preceding update must depart. For each live tree,
landed beetles attack or leave independently, with probabilities
determined by tree status. See Section V.H for details.

**E.5.4 Tree Death from Beetle Penetrations**

Following **(**Hodges et al. 1979) we assume the number of beetle
penetrations, per square meter of bole surface area, needed to kill a
tree is a normal random variable, W, having mean μ = 100 and standard
deviation σ = 10. This means, for instance, that

Pr(a tree dies from 110 or fewer penetrations per m^2^ bole surface
area) =

Pr(W ≤ 110) = 0.8413. Daily computations of tree mortality using code
subroutine **UPDATE** proceed as in this

**Example:** Suppose that in preceding days, surviving tree T has
already suffered a total of 100 successful penetrations per unit bole
area. If by day's end the new total is 110, then

Pr(T dies today) = Pr(W ≤110W \>100) = Pr(100 \< W ≤110)/Pr(W \>100)

= (0.8413 - 0.5)/0.5 = 0.6826.

Section V.J.2 contains a more general discussion.

##### 

##### E.5.5 Summing up the Day's Activities - Subroutine ENDOFDAY

For each tree, this subroutine computes

-   The total number of eggs laid during the day \-- based on the
    proportion eggpro, computed above in subroutine **EGG**

-   The number of parents that reemerge \-- based on Coulson et
    al. (1978) and Gagne et al. (1982), the *proportion* of parents that
    reemerge is assumed equal to eggpro

-   The number of remaining parents that survive

-   The new productions of attracting and repelling pheromones
    \--proportional, respectively, to the number of successful tree
    penetrations by female SPB during the day, and to the number of
    males in the tree

-   The day's quantity of resin flow from the tree \-- dependent on tree
    resin level and the number of beetle attacks that day

-   The new maturity levels attained by beetles in each juvenile cohort

-   The cumulative number of newly mature brood adults that emerged
    during the day's four updates

Based on these computations, each tree's list of juvenile cohorts is
updated by inserting the new eggs and removing those juveniles that
emerged upon reaching maturity. (These newly emerged brood adults are
added to the group of departing adults, as described above.) Finally,
the trees that succumb to beetle attack are identified and killed. See
Sections V.J.2 and V.K for details.

##### E.5.6 Moving the Departing Adult Beetles - Subroutine DISPERSE

#####  

##### As suggested in E.5.2 above, development of cohesive infestations depends on movement of most departing beetles toward roughly the same stand location. This is accomplished in DISPERSE by moving all beetles departing from each particular tree (T) to the tree situated at most one row or column removed from T that currently has the highest positive (attracting) pheromone level. Details are in Section V.M.

##### 

##### Coda:

#####  

##### This completes the general description of the SPBLOB source code structure and the tasks performed by individual code subroutines. For user convenience, this source code has been used to generate both an executable interactive code for single stand simulations in which users have wide latitude in choosing values of stand attributes - such as number of rows and columns, and time of harvest - and an executable version that simulates multiple, independent replications of a more limited list of stand scenarios. These two executable codes differ considerably in input data requirements and in output produced, and for this reason are discussed separately. Users interested in the Interactive version should consult Section III, while those wishing multiple replications are referred to Section IV. 

**III. USER GUIDE to Implementation of the Interactive**

**Windows-Based SPBLOB Code**

**A. Introduction and Overview**

**SPBLOB** is a software package that simulates the joint population
dynamics of Southern Pine Beetle and Loblolly Pine, from time of stand
planting until harvest. This User Guide contains a general description
of how the code works, along with notes and suggestions linked to the
on-screen responses needed to use this code interactively.

> An executable version of the **SPBLOB** interactive code is available
> at
>
> **www.srs.fs.usda.gov/4501/SPBLOB/Interactive,**
>
> while an electronic copy of the source code can be downloaded by going
> to

**www.srs.fs.usda.gov/4501/SPBLOB/spblob_1_25_2007.**

In addition, there is a **Multiple-Simulation** version available for
users interested in investigations that require running several
independent replications of particular stand plantings, typically for
comparison of proposed management scenarios or for other research
purposes. See Section IV below for details.

In the rest of this section, the discussion concerns the **Executable
Interactive Version** of the code. Section II above contains background
information concerning the Code Structure. For those users who want to
tailor **SPBLOB** to their own specific uses we include in Section V
below some comments about the **Source Code**.

**B. Downloading and Installing SPBLOB on Your Machine**

This section describes the steps for downloading and installing the
**Windows-based executable Interactive version** of the **SPBLOB** code.
In these notes, the code file itself is called **tree.exe**, while the
file **input.txt** consists of a list of 7000 randomly generated pairs
of seeds useful for running code simulations. Both files are contained
in a folder called **Interactive**. In order for this version of the
code to run properly, these two files *must* remain together in the
**Interactive** folder.

1\. Begin by downloading the folder **Interactive** from the website
cited above

(this folder contains both the **tree.exe** and **input.txt** files) and
save it in a

folder or directory of your choice. For example, this figure shows the
files

are saved in a folder named "interactive" in "D" drive (D directory).

![](media/image2.jpeg){width="6.013888888888889in"
height="3.1805555555555554in"}

You are now set up to run an interactive simulation of **SPBLOB**. On
some machines, you can begin a session by simply clicking on the icon
for \'tree.exe\'. If this succeeds, skip to screen #5 below. Otherwise,
first complete steps 2-4.

2\. To begin a simulation session, click "Run" in the "Start" menu.

![](media/image3.png){width="2.4375in" height="2.1875in"}

3\. In the "Run" window, type "cmd" to run the "Command" window. Click
"OK" to

finish.

![](media/image4.png){width="3.645138888888889in"
height="1.8854166666666667in"}

4\. In the "Command" window, type "D:" to access the "D" drive, then "cd

interactive" to change to the \'interactive\' folder containing the two
code

files.

![](media/image5.jpeg){width="6.013888888888889in"
height="2.986111111111111in"}

5\. As an example of code simulation, type "tree 10" to start the
simulation code.

Upon hitting \"Return\" you see this window on the screen:

![](media/image6.jpeg){width="6.14375in" height="3.2104166666666667in"}

If you answer \"no\" or \"n\", the session terminates.

If you answer \"yes\" or \"y\", you are presented with a series of
on-screen

requests that allow you to specify attributes of the stand you wish to

simulate. These are discussed in the next section.

**Notes**:

#####  1. For convenience in running simulations, the file input.txt contains a list of 

#####  7000 number pairs, each of which consists of two seeds with which a 

#####  simulation can begin. These 7000 pairs were randomly generated and 

#####  provide the potential for a large number of independent simulations. If you 

#####  wish to repeat a particular simulation, you need only begin by typing the 

#####  same tree-space-integer combination you used previously.

2\. The integer "10" typed above indicates you want to begin the
simulation

using the two seeds in the tenth row of **input.txt**. This number can
be

replaced with any integer between 1 and 7000, inclusive, to start
another

simulation which uses a different pair of seeds to initialize the random

> number generators used in the code. As an alternative, you can type
> \"tree\" alone,
>
> with no number following. In this case, the simulation will begin by
> using a randomly chosen line from **input.txt**. (The line chosen will
> be indicated in the output file generated after simulation is
> complete; see **Sections C.2.1.1** and **C.2.1.2** below.)

**C. The Interactive, Input Phase of SPBLOB**

Roughly speaking, an interactive session using **SPBLOB** consists of
two major phases:

-   An **Input Phase**, in which you specify the particular features you
    want your simulated stand to possess

-   An **Output Phase**, in which you receive yearly information about
    the status of the stand as the simulation proceeds, and a detailed
    overall summary at the end. Some of this information appears on the
    screen, but a more complete version is saved in two output files.

**C.1 The Input Phase**

Let us assume you answered \"yes\" to the question in screen \# 5 above.
In order to set up the simulation with the type of stand you wish, you
will be asked to respond to a series of questions that appear on your
monitor. In each case, set your preferred value(s) and press **Return.**

**Note:** You can terminate the session at any time by closing the
window.

Temperature is the most important determinant of SPB vital rates
(movement, mortality, oviposition, survival). For this reason,
**SPBLOB** first asks you to specify the mean annual temperature at your
stand location.

Yearly mean temperature (degrees Celsius) at stand location: e.g. 15.00

Input your value:

Next, there are questions about stand duration and geometry.

Number of Growing Seasons: e.g. 35

Input your value (must be at least 9 years):

Spacing between rows (feet): e.g. 9.0000

Input your value:

Spacing between columns (feet): e.g. 9.0000

Input your value:

As indicated in Section II, in **SPBLOB** each 'tree' is a pair of
parameter

values:

-   **g0** = tree growth rate, the 'site index' associated with that
    tree

-   **RMAX** = the maximum rate at which resin flows from holes in the
    tree bark produced by beetles (decreases when the tree is attacked
    by beetles)

The code always generates resin flow rates using the \'ordinary\'
gamma(α,β) distribution described in Roberds et al. (2003).

On the other hand, growth rates are randomly and independently generated
using a normal distribution for which you can specify the mean μ
(equivalent to a Site Index - see the Note below).

Thus, the next question is:

Stand Site Index (feet): e.g. 80.0000

Input your value:

**Note \--** The \'site index\' of the tree represents the expected
(average) height, in feet, the

tree will attain by age 25

In **SPBLOB**, there are three principal stand types available:

-   \'Ordinary\' stands, planted with improved tree materials derived
    from any source,

> e.g., improved stock from a seed orchard

-   Stands planted using \'clonally generated\' materials, the clonal
    ramets that

make up the stand being randomly mixed before planting

-   Stands planted using clonally generated materials, the clonal ramets

being planted in blocks, all ramets in any particular block belonging

to the same clone

Thus, the next question is:

Will the stand be planted using clonal materials? (Yes or No):

########  If you answer 'no' or 'n', the screen skips to the questions in the subsection headed \'C.1.5 No Clones\', near the end of this section C, and you can now move directly to that point in this guide.

If you answer 'yes' or 'y' (clonal materials should be used), you will
be

asked to respond to a series of questions about these clones. The first
is:

Will the clone ramets be planted in clonal blocks? (Yes or No):

######## If you answer 'no' or 'n', the screen skips to the section below headed \'C.1.1 Number of Clones When Blocks Are Not Used\'. Please move to that point in this guide and continue from there.

########  If you answer 'yes' or 'y', there are questions concerning these blocks:

In any simulation involving clonal blocks, the stand is created by
specifying the size and arrangement of the individual blocks. All blocks
in any particular simulation are the same size; however, this size can
vary from one simulation to another. The questions are:

Single block size

Number of rows: e.g. 16

Input your value:

Number of columns: e.g. 16

Input your value:

Number of blocks

In row direction: e.g. 2

Input your value:

In column direction: e.g. 2

Input your value:

At this point, you see a message on the screen which indicates the
effect on stand size of your choices so far:

\*\*\* Your stand size is: x by y,

where x = (# rows per block)\*(# blocks in the row direction)

and y = (# columns per block)\*(# blocks in the column direction)

You can now skip to the heading 'C.1.2 More Information Related to
Blocks'.

**C.1.1 Number of Clones When Blocks Are Not Used**

At this point, we need to specify the number of clones, and the
parameters that define these clones. Acceptable values of these
quantities depend on stand size and on whether blocks are to be used. If
there are no clonal blocks, you see these questions on the screen:

Stand size:

Number of rows: e.g. 32

Input your value:

Number of columns: e.g. 32

Input your value:

How many clones do you want to use? (must be a factor of \[# rows\]\*\[#
columns\])

Input your value:

**Note**: The above restriction is of no real consequence. For example,
if your stand contained 32 rows and 32 columns, the requirement would
read "must be a factor of 1024", allowing numbers far in excess of the
practical limit of 30-40 clones (see, e.g., Libby 1982; Bishir and
Roberds 1995, 1997, 1999; Roberds and Bishir 1997; Roberds et al. 1990).
It appears only because in **SPBLOB** we want to plant the same number
of ramets of each clone used. This requires the number of clones to be a
divisor of the total number of trees in the stand.

If you are not using clonal blocks, you can skip the next paragraph.

**C.1.2 More Information Related to Blocks**

When your stand is formed using blocks of clones, you see:

How many clones do you want to use? (must be a factor of w)

Input your value:

**Note:** The number w that appears on your screen is equal to

w = total \# of blocks

= (# blocks in the row direction)\*(# blocks in the column direction)

because in **SPBLOB**, if clones are planted in blocks, each clone is
used in the

same number of blocks.

**C.1.3 Clonal Parameter Values**

Now, some questions about clonal parameter values:

Do you want to specify the growth and resin flow parameters for each
clone? (Yes or

No):

########  If you answer 'no' or 'n', you see 

######## \*\*\* The growth and resin flow rates for the clones are:

followed by a list of randomly generated values for your stand.

**Note**: To reflect the fact that clones are typically chosen due to
their elevated

growth rates, the values of g0 are generated by censoring the normal
distribution you specified above so that only the top 15% - those with
growth

rates exceeding μ + 1.065σ - are accepted.

User assigned parameter values, on the other hand, are unrestricted.
Thus,

If you answer \'yes\' or \'y\', indicating that you prefer to set your
own growth and

resin flow rates (g0 and r0, respectively - see Section B.2 of Chapter V
for units),

you see (repeated for each integer x from 1 up to n = number of clones):

Enter your desired parameters for clone \# x:

Input a value for g0:

Input a value for r0:

Finally, whether **SPBLOB** generated the clone parameters, or you did,
clonal ramets must be deployed in the stand.

If your clones are **not** planted in blocks, there is one more
question.

Should the clones be deployed randomly in the stand? (Yes or No):

> If you answer 'yes' or 'y', **SPBLOB** will randomly mix clonal ramets
> throughout the
>
> stand, taking care to plant about the same number of ramets of each
> clone. The
>
> pattern of assignment of clone numbers will then appear on the screen,
> followed
>
> by the beginning of simulation output. (Details are in Section C.2
> below.)

If, however, you answer 'no' or 'n', SPBLOB will ask you to assign a
clone number

to every tree, one at a time. Since this can be quite time-consuming (a
32x32

stand will require over 1000 choices) this choice is to be avoided
unless there are

strong reasons for pursuing it.

In either case, you can skip in this guide to the heading \'C.1.4 Wrapup
When Clones Are Being Used.\'

If, however, your clones **are** to be deployed in blocks, there is one
additional query:

Should the clones be assigned randomly to blocks? (Yes or No):

######## If you answer 'yes' or 'y', the pattern of this placement will appear on the screen, where you will see: 

\*\*\* The arrangement of clonal blocks in your stand is:

######## followed by an array of numbers indicating the stand positions in which **SPBLOB** has randomly placed the clones (each clone will appear in the same number of blocks). 

########  If you answer \'no\' or \'n\', you can assign clones to blocks in whatever pattern and in whatever frequencies you wish. For each block (x,y), you will be asked for the:

Clone number for block x,y:

where x runs from 1 to nrsb = number of rows of blocks

and y runs from 1 to ncsb = number of columns of blocks

**Note:** Now that growth and resin flow rates have been specified for
each clone,

**SPBLOB** will plant the individual trees in the stand according to the
instructions

you have given. Whether ramets are randomly distributed in the stand or
are

planted in clonal blocks, each ramet of a particular clone type is
assigned the

*same* pair of growth and resin flow values as the other ramets within
its clone

(these values differ from one clone to another). However, some variation
is

achieved between ramets of the same clone by adding randomly and

independently distributed deviations to each tree's growth and resin
flow values,

thus modeling heterogeneous environmental factors (site variation) at
different

stand locations.

**C.1.4 Wrapup When Clones Are Being Used**

This ends the input required when clones are used. Whether **SPBLOB**
assigned clones to blocks, or you did, simulation now begins
immediately, and you see these two lines on the screen:

Stand setup is complete, simulation in progress

Clone \# = 0 indicates the whole stand

succeeded by a report on stand status as of the end of the juvenile
stage (end of year 8). Following this, there will be pauses between the
subsequent yearly reports as computation proceeds.

A detailed discussion of simulation output appears in Section C.2 below.
Here, we note only that output lines labeled \'clone \# = 0\' contain
reports for the whole stand, while lines designated \'clone \# = n\',
where n is positive, report output for clone \# n alone.

**C.1.5 No Clones**

When a stand is planted using improved materials, but not clones, we

need only specify the size of the stand:

Overall stand size:

Number of rows: e.g. 32

Input your value:

Number of columns: e.g. 32

Input your value:

**C.2. Simulation Output**

**SPBLOB** will now plant the stand you have specified, assigning
randomly and independently chosen growth and resin flow values, as well
as adding randomly and independently distributed deviations to each
tree's growth and resin flow values, thus modeling heterogeneous
environmental factors (site variation) throughout the stand. On the
screen, you will see:

Stand setup is complete, simulation in progress

Clone \# = 0 indicates the whole stand

followed by a report on stand status as of the end of the juvenile stage
(end of year 8).

**Note:** Clone \# = 0 is used to indicate the entire stand, whether or
not clones are

actually used in stand setup.

The first line on the screen precedes the data from year 8 and looks
like:

**Stand Age Clone \# Merch. Vol. Percent Surv. Trees Immi**

This line will be followed by yearly reports on the screen about the
status of the stand as the simulation proceeds. In these reports,

-   **Stand age** is expressed in **years** since stand planting

-   **Clone \#** = 0 corresponds to the entire stand; if clones are
    used, remaining lines

> run from 1 to numcln (the total number of clones used in this
> simulation)

-   **Merch. vol**. is the total merchantable timber volume (ft^3^/ac)
    currently available in the stand (**Clone \#** = 0) or from trees in
    each clone (**Clone \#** = 1 to numcln)

-   **Percent surv. trees** is the percentage of that clone\'s planted
    trees that still survive

-   **Immi** = \'false\' appears in the first line (**Clone \#** = 0) if
    there was no beetle immigration that year; if beetles did arrive,
    then **Immi** = \'true\' in that line. The value of **Immi** does
    not appear in the lines corresponding to individual clones.

**Notes:**

-   All the above information is also written into an output file - see
    Section C.2.1 below.

-   On a PC, a 32x32 stand will require 1-10 minutes of machine time.

Average run time is roughly proportional to the total number of trees

in the stand. Time differences between simulations of same-size stands

are due primarily to random variations in the number of stand years
during

which beetles migrate into the stand.

-   In principle, there is no limit on stand size. However, practical
    limitations on the amount of output storage space and on the time
    required for simulation, mitigate against very large stands.

```{=html}
<!-- -->
```
-   **A caveat** - We consider the model and its accompanying simulation
    code primarily as investigative tools, rather than accurate
    predictors of stand timber volume or degree of beetle activity in
    *particular* stand scenarios. The reason is that primary beetle
    activities - egg production, development in juvenile stages,
    movement, and mortality - are largely controlled by temperature
    (Payne 1980; Coulson 1980; Wagner et al. 1983/1984; Wagner et al.
    1984; Stephen and Lih 1985; Ungerer et al.1999) which, in turn, is
    subject to unpredictable chance fluctuations. In addition, field
    observations indicate that in some stands beetles can appear in
    considerable numbers (e.g., via migration from outside) by age 10 or
    before, while others remain free of beetle infestations throughout
    their existence. Thus, model comparisons of management strategies
    (e.g., clonal plantings vs stands generated from seedlings or, in
    clonally generated stands, block plantings vs random mix of ramets)
    should be based, not on single stand simulations, but rather on
    averages obtained from multiple, independent replications (discussed
    in Section IV below) of any proposed strategy. Comparison of
    individual replications can, however, provide information about the
    variation in stand output to be expected among plantings. And single
    stand simulations by themselves can yield useful insights into
    beetle activities and their consequences as a stand ages.

**C.2.1 What Output Information Is Saved to Files**

For each stand simulation, the code saves output information in two
files that are automatically created and placed in the folder
\'**Interactive**.\' The first of these, identified by \"t\", followed
by the simulation number (line \# from **input.txt**), stores
information about each individual tree in the simulated stand. The
second output file, labeled \"s\", followed by the simulation #,
contains yearly summaries about the stand. These outputs include all the
information that appears on your screen (see above) plus several
additional entries. Here are the details:

**C.2.1.1 The Tree-information File**

Each row in a \"t-file\" contains 9 entries that correspond to a
particular tree in the simulation. These entries appear in the sequence:

**\[seed pair, i, j, RMAX(i,j), g0(i,j), clntyp(i,j), yrdied(i,j),
daydied(i,j),**

**lightning(i,j), vmerch(i,j)\]**

where

**seed pair** is the **line number** used from the file **input.txt**

**i** is the tree\'s **row number** in the simulated stand

**j** is the **column number** of the tree

**RMAX(i,j)** = **maximum resin flow rate** (including site variation)
for tree (i,j)

**g0(i,j)** = **growth rate** (including site variation) for tree (i,j)

**clntyp(i,j)** = **number of the clone** to which tree (i,j) belongs

(= 0 if clones are not used in this simulation)

**yrdied(i,j)** = **stand age** (years) at which tree (i,j) died (= 0 if
tree is still

alive at harvest)

**daydied(i,j)** = **day of year** (1 = Jan 1; 365 = Dec 31) on which
tree (i,j) died

(= 0 if tree died from competition or is still alive at harvest)

**lightning(i,j)** = **cause of death** of tree (i,j) (see Notes below)

**vmerch(i,j)\] = merchantable volume** (ft^3^) of tree (i.j) (= 0 if
the tree is dead)

**Notes Related to the Tree-information File:**

\-- In each t-file, there is a line for each tree, each with 10 entries
(roughly 1000 lines for

a 32x32 stand).

\-- Values of **RMAX, g0** and **clntyp** are assigned at the beginning
of the simulation and constitute part of the stand initialization. On
the other hand, entries in **yrdied,**

**daydied** and **lightning** are determined as the simulation
computations proceed,

and are not known for all trees until the simulation ends.

\-- Values of **yrdied, daydied and lightning** contain more information
than

is at first apparent. For instance, **yrdied(i,j)** = 0 only if tree
(i,j) is still

alive at time of harvest, **yrdied(i,j)** = 8 indicates a tree that died
during the

8-year tree juvenile period during which no beetles were present in the

stand, and **yrdied(i,j)** = k \> 8 means that tree (i,j) died during
stand year k.

\-- Moreover, **daydied(i,j)** \> 0 if tree(i,j) was killed by beetle
penetrations **and/or**

by lightning. For example, **(yrdied(i,j)** = 21, **daydied(i,j)** =
154,

> **lightning(i,j)** = 16) will be counted as meaning tree(i,j) was
> killed by beetles **and** by lightning - on day 154 (June 3) of stand
> year 21, 16 days after being struck. On the other hand,
> **(yrdied(i,j)** = 21, **daydied(i,j)** = 154, **lightning(i,j)** =
> 30) means the tree was killed by lightning but not by beetles, while
> **(yrdied(i,j)** = 21, **daydied(i,j)** = 0) means that tree (i,j)
> died from tree competition at the beginning
>
> of year 21, prior to that year\'s daily updates, and was not killed by
> beetles.

\-- To clarify further the information supplied by the **lightning**
matrix,

**lightning(i,j)** = 0 if tree (i,j) was never struck by lightning
during the entire

course of the simulation; **lightning(i,j)** = n, where 0 \< n \< 30, if
tree (i,j)

WAS struck by lightning, then died n days later as a consequence of

beetle attacks; and **lightning(i,j)** = 30 means the tree was struck by

lightning, not killed by beetles during the following 30 days, and died

(by default, in the code) due to the original lightning strike.

\-- In short, from values of **yrdied, daydied** and **lightning** we
can determine

*when* (year and day) and *why* (competition, lightning, beetles) each
tree

died, if at all. Here is the list of possibilities:

-   **yrdied(i,j)** indicates the year of death, if any: **yrdied(i,j)**
    = 0 if tree (i,j)

survives until harvest

-   if tree(i,j) dies, then **yrdied(i,j)** \> 0; in this case,

    -   **daydied(i,j)** = 0 means tree(i,j) was killed by competition

    -   **daydied(i,j)** \> 0, **lightning(i,j)** = 0 means killed by
        beetles, not lightning

    -   **daydied(i,j)** \> 0, **lightning(i,j)** = n (0 \< n \< 30)
        means killed by beetles **and** lightning

    -   **daydied(i,j)** \> 0, **lightning(i,j)** = 30 means killed by
        lightning, not beetles

**C.2.1.2 The Stand-information File**

In an \"s-file\", each row contains 18 entries that correspond to a
particular clone in a particular year, within a particular simulation.
These entries are arranged in the sequence:

**\[seed pair, age, immi, clone #, %treesurv, %treekill_B,
%treekill_L,**

**%treekill_C, avdiam_surv, avdiam_KB, avdiam_KC, avhgt_surv,**

**avhgt_KB, avhgt_KC, avres_surv, avres_KB, merchvol, basal_area\]**

where **seed pair** is as above, and

**age** = stand age in years

**immi** = 1 if there is beetle immigration this year, and = 0 if no
immigration

**clone \#** = number of the \'clone\' about which info is reported (0 =
entire

> stand)

**%treesurv** = for each clone, the percentage of planted trees that
still survive

**%treekill_B** = for each clone, the percentage of planted trees killed
by beetles

*this* year

**%treekill_L** = for each clone, the percentage of planted trees killed
by

lightning *this* year

**%treekill_C =** for each clone, the percentage of planted trees killed
by

competition *this* year

**avdiam_surv** = for each clone, the average diameter (in) for
surviving trees

**avdiam_KB =** for each clone, the average diameter (in) for trees
killed by

> beetles *this year* (\*\*)

**avdiam_KC =** for each clone, the average diameter (in) for trees
killed by

competition *this year* (\*\*)

**avhgt_surv** = for each clone, the average height (ft) for surviving
trees

**avhgt_KB =** for each clone, the average height (ft) for trees killed
by beetles

*this year* (\*\*)

**avhgt_KC =** for each clone, the average height (ft) for trees killed
by

> competition *this year* (\*\*)

**avres_surv =** for each clone, the average resin flow (RMAX: gm/24h
from a

> 0.5\" diameter wound**)** for surviving trees

**avres_KB =** for each clone, the average resin flow (RMAX: gm/24h**)**
in trees

> killed by beetles up to now

**merchvol** = for each clone, the total merchantable timber volume
(ft^3^/acre) in

> its surviving trees

**basal_area** = for each clone, the basal area (ft^2^/acre) in its
surviving trees

**Notes Related to Stand-information Files:**

\-- In each s-file, in each year from age 8 to harvest, there is a line
for each clone.

With an average of 15 lines (1 + \# of clones) per year and 28 years of
output

(if harvest age = 35), an average s-file contains approximately 15x28 =
420 lines,

each with 18 entries.

\-- In the code, although competition deaths, tree diameters and tree
heights are

computed at the *beginning* of each year, all information recorded in
s-files

is current as of the *end* of each year, in order to also take account
of the

year\'s tree mortality due to beetles.

\-- The four percentages listed above can total more than 100. The
reason is

that trees can be killed by both beetles and by lightning (in the code,
this

occurs when a lightning struck tree is subsequently killed by beetles

before the default 30 day period ends - see the 4th note related to the

t-files, above).

\-- Juvenile tree deaths are listed in addition to those due to beetles,
competition

and lightning in order for the percentages to include all planted trees.

\-- Items marked by (\*\*) above differ from other closely related items
in that they report

on the current year, rather than being cumulative from time of stand
planting.

> This is because heights and diameters of killed trees can be
> meaningfully

compared to those of surviving trees only within the year of death.

\-- The merchantable volume in the stand is the sum of volumes of all
*surviving* trees,

> computed on a per acre basis. For an individual clone, the same
> computation is made, though here the divisor used is the area occupied
> only by the originally planted ramets of *that* clone, not the area of
> the entire stand.

\-- The basal area in the stand is the sum of all the cross-sectional
areas of the

> surviving trees, computed on a per acre basis. For an individual
> clone, the same computation is made, though here the divisor used is
> the area occupied only by the originally planted ramets of *that*
> clone, not the area of the entire stand.

\-- **Warning**: When clonal ramets are randomly mixed before planting,
relative

magnitudes of the per-acre volumes or basal areas reported may, due to
chance

proximity of clones with different growth rates, not reflect true
comparisons of

volumes these clones would produce if planted in individual blocks.

##### IV. GUIDE TO USING THE EXECUTABLE MULTIPLE-SIMULATION 

#####  VERSION OF SPBLOB

##### 

##### A. How to Run Batch Simulation Jobs

Four files are required in order to run the version of **SPBLOB**
designed for batch simulation jobs. These files, contained in a folder
called **multiple**, are:

**configs.txt input.txt**

**run.txt tree.exe**

> An executable version of this Multiple-Simulation code is available at
>
> **www.srs.fs.usda.gov/4501/SPBLOB/Multiple**
>
> while an electronic copy of the source code can be downloaded by going
> to

**www.srs.fs.usda.gov/4501/SPBLOB/splob_1_25_2007.**

(This source code is the same as that used to generate the executable
Interactive Code discussed in Section III. If you wish to modify the
source code, e.g., to study effects of changing certain parameter
values, open the relevant subroutine files and make the desired changes
by editing them as you would any text file. After making your changes,
set the variable \"**batch**\" in the subroutine **MAIN** to \"1\" to
rebuild the **Multiple** simulation version. Then compile the resulting
new **Source Code**. See Section V for full details)

The first task is to place the \'**multiple\'** folder in a directory.
Any directory can be used; for convenience, we shall refer to it as
\'Directory D.\' In Section A.4 below we indicate how to run batch
simulations using this setup. First, though, a brief description of the
contents of each file.

**A.1 The File Configs.txt**

Each of the 23 rows of the file **configs.txt,** shown here as **Table
1**, contains 13 entries. The first two values serve as **row labels**,
while the last eleven are input parameters that constitute the
particular **stand setup**, or **configuration**, in that row. The
**labels** are, respectively, the row number of the configuration and a
code related to the overall pattern in which trees are planted. In this
code, the 1 or 2 digits (0, 1, 2, 6, 12, 18, 30 or 50) *preceding* the
last digit represent the number of clones, if any, that are used in
planting the stand (0 = no clones). If the number of clones is 0 or 1,
the last

digit identifies the mean yearly temperature (3^rd^ column of the
table). Otherwise, the *last digit* indicates the pattern in which
clonal materials are deployed:

1= all trees are randomly mixed before planting

2 = for each clone, all trees are planted in a single block

3 or 4 = trees of each clone are planted in more than one block

> (see parameter bspc in the list below)

##### Table 1. List of preprogrammed configurations for SPBLOB simulations

  ------- --------- ------- ---- ---- -------- ------ ------ ------ ------ ------ ---- ----
  setup   pattern   tmean   nr   nc   numcln   rspb   cspb   nrsb   ncsb   bspc   px   py

  1       01        15      32   32   0        0      0      0      0      0      3    3

  2       02        17      32   32   0        0      0      0      0      0      3    3

  3       03        19      32   32   0        0      0      0      0      0      3    3

  4       11        15      32   32   1        32     32     1      1      1      3    3

  5       12        17      32   32   1        32     32     1      1      1      3    3

  6       13        19      32   32   1        32     32     1      1      1      3    3

  7       21        17      32   32   2        1      1      32     32     512    3    3

  8       22        17      32   32   2        32     16     1      2      1      3    3

  9       23        17      32   32   2        16     16     2      2      2      3    3

  10      61        17      36   30   6        1      1      36     30     180    3    3

  11      62        17      36   30   6        12     15     3      2      1      3    3

  12      63        17      36   30   6        9      10     4      3      2      3    3

  13      64        17      36   30   6        6      5      6      6      6      3    3

  14      121       17      36   30   12       1      1      36     30     90     3    3

  15      122       17      36   30   12       9      10     4      3      1      3    3

  16      123       17      36   30   12       6      5      6      6      3      3    3

  17      181       17      36   30   18       1      1      36     30     60     3    3

  18      182       17      36   30   18       6      10     6      3      1      3    3

  19      183       17      36   30   18       4      5      9      6      3      3    3

  20      301       17      36   30   30       1      1      36     30     36     3    3

  21      302       17      36   30   30       6      6      6      5      1      3    3

  22      501       17      40   30   50       1      1      40     30     24     3    3

  23      502       17      40   30   50       4      6      10     5      1      3    3
  ------- --------- ------- ---- ---- -------- ------ ------ ------ ------ ------ ---- ----

The eleven **input parameters** are:

-   tmean = average annual temperature (^o^C) at the stand location

-   nr = total number of rows in the stand

-   nc = total number of columns in the stand

-   numcln = number of clones used in planting the stand (0 = no clones)

-   rspb = number of rows of trees in each clonal block

-   cspb = number of columns of trees in each clonal block

-   nrsb = number of rows of blocks, if any, in the stand

-   ncsb = number of columns of blocks, if any, in the stand

-   bspc = for each clone, the number of blocks, if any, in which it is
    planted

-   px = distance (m) between columns in the stand

-   py = distance (m) between rows in the stand

**Notes:**

-   In simulating these configurations, clonal and/or individual tree
    values of

> g0 and RMAX are generated randomly (see the discussions in
>
> Sections II.C.2 above and V.B below)

-   numcln = 0 indicates that no clones are used

-   numcln \> 0 implies use of clones in planting the stand

-   If clonal blocks are used, all blocks in a particular simulation
    have the

same size and shape, and each clone is planted in the same

> number of blocks

-   The entries in **Table 1** can be varied to allow simulation of
    stands with

> different overall sizes and internal geometries. For instance, if you
>
> want to run simulations of stand setup \# 2, but for a larger stand
>
> with 60 rows and 50 columns, simply replace nr = 32 and nc = 32
>
> by, respectively, nr = 60 and nc = 50. Then save the file
>
> **configs.txt** and run the simulations in the usual way, as described
>
> in **Section A.4** below.

-   **Warning** - There are restrictions on the changes allowed in
    **configs.txt.**

> For instance, if you want a stand with 7 rows of clonal blocks, and
>
> 6 columns of blocks, each block containing 5 rows and 6 columns
>
> of trees, then your stand size must be 7x5 = 35 rows by 6x6 = 36
>
> columns. Further, the total number of blocks in the stand
>
> is 7x6 = 42, which requires the number of clones to be a divisor of
>
> 42 (one of the numbers 2, 3, 6, 7, 14, 21 or 42) in order for each
>
> clone to appear in the same number of blocks. If any of these
>
> restrictions are violated, the code gives an error message and
>
> simulation will not proceed.

**A.2 The File Input.txt**

File **input.txt** contains a list of 7000 randomly and independently
generated pairs of integers, the first 15 of which are exhibited in
**Table 2**. As described below, batch simulations are constructed by
choosing groups of pairs (i.e., rows) from this table. Each particular
pair leads to a different simulation, the two seeds serving to
initialize random number streams that generate 1) the growth and resin
flow parameter values (g0 and RMAX) associated with each tree and 2) the
random sequence of daily temperatures, for the simulation.

##### Table 2. The first 15 rows in the file input.txt

  --------- ------------------------------ ------------------------------
  pair      first seed                     second seed

  1         769507989                      391973004

  2         786410668                      766832987

  3         557965049                      908347825

  4         143157914                      657742113

  5         921209698                      883037717

  6         24862566                       688162812

  7         584741353                      715680511

  8         785792706                      505503855

  9         980069746                      105311691

  10        520294406                      844090295

  11        189630999                      754658824

  12        519544887                      172030176

  13        99964747                       909842805

  14        985824992                      19186678

  15        272832311                      687032720
  --------- ------------------------------ ------------------------------

**A.3 The File Run.txt**

When generating multiple, independent replications of a particular stand
setup, it is convenient to enter a single command rather than having to
specify each simulation separately. The file **run.txt,** shown here as
**Table 3,** facilitates this approach.

**Table 3. The file run.txt. Beginning with the third row of the table**

> **(setup \# 1) row entries under the headings \"First Line\"**
>
> **and \"Last Line\" designate the line numbers in the file**
>
> **input.txt that are reserved for simulation of the stand**
>
> **setup indicated in the row.**

  -------- ---------- ----------------- ---------------- -----------------
  setup    pattern    Lines Reserved in                  Harvest Age
                      input.txt                          

                      First Line        Last Line        

  1        01         1001              1200             35

  2        02         1201              1400             35

  3        03         1401              1600             35

  4        11         1601              1800             35

  5        12         1801              2000             35

  6        13         2001              2200             35

  7        21         2201              2400             35

  8        22         2401              2600             35

  9        23         2601              2800             35

  10       61         2801              3000             35

  11       62         3001              3200             35

  12       63         3201              3400             35

  13       64         3401              3600             35

  14       121        3601              3800             35

  15       122        3801              4000             35

  16       123        4001              4200             35

  17       181        4201              4400             35

  18       182        4401              4600             35

  19       183        4601              4800             35

  20       301        4801              5000             35

  21       302        5001              5200             35

  22       501        5201              5400             35

  23       502        5401              5600             35
  -------- ---------- ----------------- ---------------- -----------------

**Note -** The lines from **input.txt** actually used re a particular
setup can be

varied by editing entries in Table 3. For instance, if you wish to run

a set of 50 independent replicates of Setup \# 15, beginning with line
(seed

pair) \# 3842, and with harvest age = 40 years, simply change the row
entries

to read: 15 122 3842 3891 40 and re-save the file.

**A.4 The File Tree.exe**

This is the executable code file that does the work of running repeated,
independent stand simulations. (**Warning -** Though it has the same
name, this file is *not* the same as the file *tree.exe* contained in
the folder \'**interactive**\' discussed in Section III above. To avoid
confusion, and to obtain proper results, these files must remain in
their respective folders.) This executable code also writes simulation
results to output files, the contents of which are discussed in
**Section B**, below. Here, we describe how to run a batch of
simulations.

As indicated in **Table 3**, seed pairs from file **input.txt** are
dedicated to each of the 23 stand scenarios listed in **Table 1**. With
this arrangement, a single command suffices to run all the replicates
associated with a particular scenario of interest. (See the **Note**
immediately preceding this section.)

For example, suppose you want to run independent replications of stand
scenario number 15. The first step is to move to the directory
containing the four code files necessary for batch simulation. If, as
above, these files are saved in "D" drive (D directory)

-   Click "Run" in the "Start" menu

![](media/image3.png){width="2.4375in" height="2.1875in"}

-   In the "Run" window, type "cmd" to open the "Command" window. then

> click "OK" to finish.

![](media/image4.png){width="3.645138888888889in"
height="1.8854166666666667in"}

-   In the Command Window, type \"D:\"; then press the enter key to
    access

> the "D" drive

-   On the screen you will see **D:\\\>** . Now type **cd multiple** to
    enter the

> folder called \'**multiple**\' that contains the four files. Finally,
> type
>
> \"**tree 15**\" (the word \"tree\" followed by the number of the
> scenario
>
> you want to simulate) and press \"enter\".

This last action starts the simulations. On the screen, a heading and
one or more rows of output data will follow immediately. The heading
looks like this:

**Stand Age Clone \# Merch. Vol. Percent Surv. Trees Immi**

The output rows that follow appear in groups, each group corresponding
to one year of stand simulation. The first group you will see represents
a summary of stand status at the end of the 8-year juvenile tree stage.

A glance at **Table 1** shows that stand setup \# 15 involves 12 clones,
planted in separate blocks containing 9 rows and 10 columns each. In
this case, each yearly group of output consists of 13 rows, the first
having 5 entries, the others 4 entries each. In these rows,

-   **Stand Age** is expressed in **years** since stand planting

-   **Clone \#** = 0 corresponds to the entire stand; remaining values
    indicate

> clones 1 through 12

-   **Merch. Vol**. is the total merchantable timber volume (ft^3^/ac)
    currently available

> in the stand (clone \# = 0) or from trees in individual clones (1-12)

-   **Percent Surv. Trees** is the percentage of a clone\'s originally
    planted

> trees that still survive

-   **Immi** = \"false\" appears in the first line (clone \# = 0) if
    there was no beetle

> immigration that year; if beetles did arrive, then **Immi** = \"true\"
> in that line. The value of **Immi** is not repeated in the lines
> corresponding to individual clones.

**Notes:**

-   On a PC, a 32x32 stand will require 1-10 minutes of machine time.

Average run time is roughly proportional to the total number of trees

in the stand. Time differences between simulations of same-size stands

are due primarily to random variations in the number of stand years
during

which beetles migrate into the stand.

-   In principle, there is no limit on stand size. However, practical
    limitations on the amount of output storage space and on the time
    required for simulation, mitigate against very large stands.

```{=html}
<!-- -->
```
-   **A caveat** - We consider the model and its accompanying simulation
    code primarily as investigative tools, rather than accurate
    predictors of stand timber volume or degree of beetle activity in
    *particular* stand scenarios. The reason is that primary beetle
    activities - egg production, development in juvenile stages,
    movement, and mortality - are largely controlled by temperature
    (Payne, 1980; Coulson, 1980; Wagner et al. 1983/1984; Wagner et al.
    1984; Stephen and Lih 1985; Ungerer et al.1999) which, in turn, is
    subject to unpredictable chance fluctuations. In addition, field
    observations indicate that in some stands beetles can appear in
    considerable numbers (e.g., via migration from outside) by age 10 or
    before, while others remain free of beetle infestations throughout
    their existence. These random fluctuations also appear in model
    simulations. Thus, model comparisons of management strategies (e.g.,
    clonal plantings vs those generated from seedlings or, in clonally
    generated stands, block plantings vs random mix of ramets) should
    not be based on single stand simulations, but rather on multiple,
    independent replications of any proposed strategy. Single stand
    simulations can, however, yield useful insights into beetle
    activities and their consequences as a stand ages.

**B. The Simulations Themselves - What Output Information Is Saved**

As described above, **Table 1** enables us to use a single setup number
to represent the list of 11 stand parameters corresponding to an
individual simulation. Together with **Table 2**, it also indicates the
overall plan for the simulations themselves. Specifically, each
simulation is uniquely identified by its setup number and the line used
from **input.txt**.

Following this format for identifying simulation inputs, for each stand
simulation the code saves output information in two files that are
automatically created and placed in the folder \'multiple.\' The first
of these, identified by t, followed by the simulation number (line \#
from **input.txt**), stores information about each individual tree in
the simulated stand. The second output file, labeled s, followed by the
simulation number, contains yearly summaries about the stand. These
outputs include all the information that appears on your screen (Section
A.4) plus several additional entries. Here are the details:

**B.1 The Tree-information Files**

For each stand simulation, a \"t-file\" is produced at harvest time, In
this file, each row contains 11 entries that correspond to a particular
tree in the simulation. These entries appear in the sequence:

**\[setup #, seed pair, i, j, RMAX(i,j), g0(i,j), clntyp(i,j),
yrdied(i,j),**

**daydied(i,j), lightning(i,j), vmerch(i,j)\]**

where

**setup \#** is the **setup number** for the simulation

(one of the entries in the first column of **Table 1**)

**seed pair** is the **line number** used from the file **input.txt**

**i** is the tree\'s **row number** in the simulated stand

**j** is the **column number** of the tree

**RMAX(i,j)** = **maximum resin flow rate** (including site effect) for
tree (i,j)

**g0(i,j)** = **growth rate** (including site effect) assigned to tree
(i,j)

**clntyp(i,j)** = **number of the clone** to which tree (i,j) belongs

(= 0 if clones are not used in this simulation)

**yrdied(i,j)** = **stand age** (years) at which tree (i,j) died (= 0 if
tree is still

> alive at harvest)

**daydied(i,j)** = **day of year** (1 = Jan 1; 365 = Dec 31) on which
tree (i,j) died

(= 0 if tree died from competition or is still alive at harvest)

**lightning(i,j)** = **cause of death** of tree (i,j) (see Notes below)

**vmerch(i,j)\] = merchantable volume** (ft^3^) of tree (i.j) (= 0 if
the tree is dead)

**Notes Related to the Tree-information Files:**

\-- In each t-file, there is a line for each tree, between 1000 and 1200
lines, each

with 10 entries, if a stand setup is chosen from Table 1.

\-- Values of **RMAX, g0** and **clntyp** are assigned at the beginning
of the

simulation and constitute part of the stand initialization. On the other

hand, entries in **yrdied, daydied** and **lightning** are determined as
the

daily simulation computations proceed, and are not known for all trees

> until the simulation terminates.

\-- Values of **yrdied, daydied and lightning** contain more information
than

is at first apparent. For instance, **yrdied(i,j)** = 0 only if tree
(i,j) is still

alive at time of harvest, **yrdied(i,j)** = 8 indicates a tree that died
during the

8-year tree juvenile period during which no beetles were present in the

stand, and **yrdied(i,j)** = k \> 8 means that tree (i,j) died during
stand year k.

\-- Moreover, **daydied(i,j)** \> 0 if tree(i,j) was killed by beetle
penetrations **and/or**

by lightning. For example, **(yrdied(i,j)** = 21, **daydied(i,j)** =
154,

> **lightning(i,j)** = 16) means tree(i,j) was killed by beetles **and**
> by lightning - on day 154 (June 3) of stand year 21, 16 days after
> being struck. On the other hand, **(yrdied(i,j)** = 21,
> **daydied(i,j)** = 154, **lightning(i,j)** = 30) means the tree was
> killed by lightning but not by beetles, while **(yrdied(i,j)** = 21,
> **daydied(i,j)** = 0) means that tree (i,j) died from tree competition
> at the beginning of year 21, prior to that year\'s daily updates, and
> was not killed by beetles.

\-- To clarify further the information supplied by the \"lightning\"
matrix,

**lightning(i,j)** = 0 if tree (i,j) was never struck by lightning
during the entire

course of the simulation; **lightning(i,j)** = n, where 0 \< n \< 30, if
tree (i,j)

WAS struck by lightning, then died n days later as a consequence of

beetle attacks; and **lightning(i,j)** = 30 means the tree was struck by

lightning, not killed by beetles during the following 30 days, and died
(by

default, in the code) due to the original lightning strike.

\-- In short, from values of **yrdied, daydied** and **lightning** we
can determine

*when* (year and day) and *why* (competition, lightning, beetles) each
tree

died, if at all. The possibilities are:

-   **yrdied(i,j**) indicates the year of death, if any: **yrdied(i,j)**
    = 0 if tree ij

survives until harvest

-   if tree(i,j) dies, then **yrdied(i,j)** \> 0; in this case,

    -   **daydied(i,j)** = 0 means tree(i,j) was killed by competition

    -   **daydied(i,j)** \> 0, **lightning(i,j)** = 0 means killed by
        beetles, not lightning

    -   **daydied(i,j)** \> 0, **lightning(i,j)** = n (0 \< n \< 30)
        means killed by beetles

> **and** lightning

-   **daydied(i,j)** \> 0, **lightning(i,j)** = 30 means killed by
    lightning, not beetles

**B.2 The Stand-information files**

At the beginning of each simulation, an \"s-file\" is created, and
information is added to the file at the end of each simulated year. Each
row in the file contains 19 entries that correspond to a particular
clone, in a particular year, within a particular simulation. These
entries are arranged in the sequence:

**\[setup #, seed pair, age, immi, clone #, %treesurv, %treekill_B,**

**%treekill_L, %treekill_C, avdiam_surv, avdiam_KB, avdiam_KC,**

**avhgt_surv, avhgt_KB, avhgt_KC, avres_surv, avres_KB, merchvol,**

**basal_area\]**

where **setup \#** and **seed pair** are as above, and

**age** = stand age in years

**immi** = 1 if there was beetle immigration this year, and = 0 if no

> immigration

**clone \#** = number of the \'clone\' about which info is reported (0 =
entire

> stand)

**%treesurv** = for each clone, the percentage of planted trees that
still survive

**%treekill_B** = for each clone, the percentage of planted trees killed
by beetles

*this* year

**%treekill_L** = for each clone, the percentage of planted trees killed
by

lightning *this* year

**%treekill_C =** for each clone, the percentage of planted trees killed
by

competition *this* year

**avdiam_surv** = for each clone, the average diameter (in) for
surviving trees

**avdiam_KB =** for each clone, the average diameter (in) for trees
killed by

beetles *this year* (\*\* see the first Note below)

**avdiam_KC =** for each clone, the average diameter (in) for trees
killed by

competition *this year* (\*\*)

**avhgt_surv** = for each clone, the average height (ft) for surviving
trees

**avhgt_KB =** for each clone, the average height (ft) for trees killed
by beetles

*this year* (\*\*)

**avhgt_KC =** for each clone, the average height (ft) for trees killed
by

competition *this year* (\*\*)

**avres_surv =** for each clone, the average resin flow (RMAX: gm/24h
from a

> 1.275 cm diameter wound) for surviving trees

**avres_KB =** for each clone, the average resin flow (RMAX: gm/24h) in
trees

> killed by beetles up to the present

**merchvol** = for each clone, the total merchantable timber volume
(ft^3^/acre) in its

surviving trees

**basal_area** = for each clone, the total basal area (ft^2^/acre) in
its surviving trees

**Notes Related to Stand-information Files:**

\-- Items marked by (\*\*) above differ from other closely related items
in that they report

on the current year, rather than being cumulative from time of stand
planting.

This is because heights and diameters of killed trees can be
meaningfully

compared to those of surviving trees only within the year of death.

\-- In each s-file, in each year from age 8 to harvest, there is a line
for each clone. With

an average of 15 lines (1 + \# of clones) per year, and harvest time of
35 years,

an average s-file would contain about 15x28 = 420 lines, each with 19
entries.

\-- In the code, although competition deaths, tree diameters and tree
heights are

computed at the *beginning* of each year, all information recorded in
s-files

is current as of the *end* of each year, in order to also take account
of the

year\'s tree mortality due to beetles.

\-- The four percentages listed above can total more than 100. The
reason is

that trees can be killed by *both* beetles and lightning (in the code,
this

occurs when a lightning struck tree is subsequently killed by beetles

before the default 30 day period ends - see the 4th note related to the

t-files, above).

\-- Juvenile tree deaths are listed in addition to those due to beetles,
competition

and lightning in order for the percentages to include all planted trees.

\-- The merchantable volume in the stand is the sum of volumes of all
*surviving* trees,

> computed on a per acre basis. For an individual clone, the same
> computation is made, though here the divisor used is the area occupied
> only by the originally planted ramets of *that* clone, not the area of
> the entire stand.

\-- The basal area in the stand is the sum of all the cross-sectional
areas of the

> surviving trees, computed on a per acre basis. For an individual
> clone, the same computation is made, though here the divisor used is
> the area occupied only by the originally planted ramets of *that*
> clone, not the area of the entire stand.

\-- **Warning**: When clonal ramets are randomly mixed before planting,
relative

magnitudes of the per-acre volumes or basal areas reported may, due to
chance

proximity of clones with different growth rates, not reflect true
comparisons of

volumes these clones would produce if planted in individual blocks.

##### 

#####  

##### V. THE SPBLOB Model SOURCE CODE, WITH EMPHASIS ON 

#####  Parameters, THEIR Values and THEIR Effects ON 

#####  MODEL OUTPUT

**A. Introduction**

This section V of the **SPBLOB** Documentation is offered for those
users who want to tailor **SPBLOB** to their own uses. It explains the
principal computations in the **SPBLOB** source code, as well as
describing the major parameters, their values and their sources in the
literature. More detailed descriptions of a few particularly complex
computations can be found in the folder

**SPBLOB_Supporting_Info_re_Params_Values**,

available at

**www.srs.fs.usda.gov/4501/SPBLOB.**

This folder will hereinafter be referred to as comprising the
'**Supporting Information**' folder.

Before continuing here, it may be helpful to review Section II, which
contains an overview of the **SPBLOB** code and provides a descriptive
summary of each of the major subroutines.

**A.1 The SPBLOB Source Code**

The **Source Code** is primarily in Fortran with, as indicated below, a
few random number generators in C. If you want to build your own version
of this code, begin by downloading a copy from:

**www.srs.fs.usda.gov/4501/SPBLOB/spblob_1_25_2007.**

This code is in the form of a folder, called **spblob_1_25_2007,** in
which individual files constitute the subroutines, functions, etc.
necessary for compilation and running of the code. In addition to these
files, this folder contains the **Interactive** and **Multiple**
subfolders that constitute the executable versions of the **SPBLOB**
code described in Sections III and IV. Properly speaking, these
subfolders are not part of the source code, but have been created from
it.

**A.1.1 Modifying and Recompiling the Source Code**

To modify the source code, open the relevant subroutine files and make
the desired changes, editing them as you would any text file. After
making your changes, if you want to re-build the executable
**Interactive** code, set the variable \"**batch**\" in the subroutine
**MAIN** to \"0\". To build the **Multiple** simulation version, set
\"**batch**\" equal to \"1\".

Various Fortran compilers are available in the public domain, and we
assume you have your own favorite. (For instance, users of Unix/Linux
compilers will probably create a Makefile for the project.) We ourselves
used the **Open Watcom** compiler, which can be downloaded from:

[**http://www.openwatcom.org/**](http://www.openwatcom.org/)

The current (February, 2007) version is 1.6, released in December, 2006.
Two installation files are required, one for C/C++ and one for Fortran,
since we use both languages in the **SPBLOB** code.

**If** you apply the **Open Watcom** compiler software to open the
project file (**tree.wpj**) it creates, you find three types of source
files (other compilers may build different components):

-   Fortran source code (for which subroutine names end with **.for**)

-   C source code (ending with **.c**)

-   Compiled lib files (ending with **.lib**)

The Fortran source codes are the major simulation subroutines. Several
of these require generation of random values from uniform, gamma,
normal, and/or exponential probability distributions. For this purpose,
we include in the code folder the random number generator library that
can be downloaded from:

[**http://biostatistics.mdanderson.org/SoftwareDownload/SingleSoftware.aspx?Software_Id=27**](http://biostatistics.mdanderson.org/SoftwareDownload/SingleSoftware.aspx?Software_Id=27)

For us, some features of the Fortran version of this library did not
work as expected and we used the C version instead. To call the random
number generator library in C source code requires use of the '\$pragma"
syntax to let the **Open Watcom** compiler know that some routines
called by the Fortran codes are written in C. If you use a compiler
other than **Open Watcom**, the '\$pragma' syntax probably will need to
be changed to fit your compiler\'s requirements. On the other hand, if
you use the Fortran version of the random number library indicated
above, the '\$pragma' syntax is not required.

The subroutines **'get_command_argument'** and
'**command_argument_count'** called, e.g., in subroutine **MAIN** are
subroutines from the lib file **f2kcli.lib**. The source file for this
lib file, along with different binary versions for different operating
systems and compilers, can be downloaded from:

[**http://www.winteracter.com/f2kcli/**](http://www.winteracter.com/f2kcli/)

If you use this library with a different compiler than **Open Watcom**,
consult the **f2kcli** website for information concerning the proper
documents to download. In **SPBLOB**, this library is used to interpret
and act on command lines such as 'tree 1' ('1' is a command line
parameter). For some compilers this library may not be necessary since
they may already include similar functions.

##### A.2 SPBLOB Parameters: Values, Sources, and Effects on Model Output

Symbols used in the **SPBLOB** code are either variables, parameters or
constants:

**Variable** A quantity that can change over the course of a simulation

**Parameter** A quantity that does not change value within a particular

simulation of a particular stand, but may change from

> one simulation to another

**Constant** A quantity that has the same numerical value in all model
simulations

Most variables and parameters are numerical, but a few are logical
quantities.

**Parameters** serve to clarify and give flexibility to the model code,
and to calibrate simulation output. Most of the parameter values used in
**SPBLOB** were obtained from field studies reported in the literature.
Values not available from published sources were then chosen to match
model output as closely as possible to field observations. The
**Glossary** in Section VI contains definitions of all the important
symbols, including parameters. Here, we focus on sources from which the
principal parameters values were obtained and, for those parameters for
which values were not available directly, the methods used for their
computation.

##### 

##### B. Parameters Associated With Stand Planting

As noted in Section II.C.2 above, in **SPBLOB** a 'tree' is a pair of
parameters that govern tree interactions with their neighbors and with
the environment. These fundamental **SPBLOB** parameters are

-   **g0**, the tree\'s growth rate, a combination of inherent genetic
    factors and the site variation associated with the particular
    location in which it is planted; during simulation, this value,
    together with the g0 values assigned to the tree\'s neighbors,
    controls the tree\'s survival and growth in the face of intrastand
    competition.

-   **RMAX**, the tree\'s inherited maximum resin flow rate, plus site
    variation, which governs its success in repelling beetle attacks.

##### 

##### B.1 The Parameter g0

All g0 values are assigned in subroutine **DEPLOY**, then used in
subroutine **JUVEN (**Burkhart et al. 1987; also see Section II) as the
stand is advanced through its 8-year juvenile growth stage. These values
also play a central role in subroutines **COMP2** and **GROW2** (derived
from Burkhart et al. 1987) that govern the competitive yearly survival
and growth of individual trees.

If clones are not used in establishing the stand, tree seedlings are
generated randomly and independently from improved stock, their g0
values coming from a normal distribution having mean 80 (average tree
height, in feet, at age 25) and standard deviation 5. These values were
chosen to represent tree growth on a site of moderately high
productivity (Schultz 1997).

If clones *are* to be used, the g0 values can be assigned by the user,
if desired. When clone values are not chosen by the user, the code
censors the Normal(80,25) distribution by accepting only those g0 values
that fall in the upper 15% of the distribution, these censored values
being independently assigned to clones. For each particular clone, all
ramets are assigned the same g0 value, whether these ramets are deployed
together in blocks or are randomly mixed with those of other clones.
Environmental variation is added to each value as the ramets are
planted.

##### B.2. The Parameter RMAX

#####  Based on a large-scale study of resin flow distributions (Roberds et al. 2003), we assume that RMAX follows a gamma distribution with parameters α = 2.5 and β = 0.6. (Equivalently, the mean of RMAX is αβ = 1.5 while the standard deviation is sqrt(αβ^2^) ≅ 0.95.) The units are gm/24h from a hole 1.27 cm in diameter.) This distribution is used whether or not the stand is generated from clonal materials. 

Predicated on results reported in Roberds et al. (2003), we assume in
assigning RMAX values to seedlings or clonal ramets that r0 and g0
values are positively correlated, with correlation coefficient ρ = 0.25.
Each g0 value is assigned first, after which the code calls subroutine
**CORRELATE** to determine the mean and variance of the gamma
distribution from which the corresponding RMAX will be chosen.
Subroutine **CLONE** then randomly chooses an RMAX value from this
distribution. Full details of this computational process appear in the
file **spblob_params_g0&r0.doc,** contained in the '**Supporting
Information**' folder (see Section A).

##### 

##### B.3 Assignment of Phenotypic Values to Trees 

Pairs of g0 and RMAX values represent tree genotypes. Tree planting then
results in altered phenotypic values that combine genetic influences
with environmental effects associated with the site at which the tree is
located. In code subroutine **DEPLOY**, we assign individual tree,
broad-sense heritability (Balocchi et al. 1993; Falconer and Mackay
1996; Roberds et al. 2003) equal to hsqg0 = 0.4 (a value consistent with
a modest level of within-site environmental variability in loblolly pine
stands). Here hsqg0 = g0var/(g0var + siteg0) is the proportion of
variation in growth rate attributable to genetic sources, of the
combined genetic and site variances. To obtain a tree\'s phenotypic g0,
we first solve this equation for the variance, siteg0, associated with
environmental effects:

siteg0 = g0var\*(1-hsqg0)/hsqg0,

then randomly sample a N(0,siteg0) distribution to generate a random
environmental variation, E. Adding E to the genotypic g0 value (Section
B.1) produces the new, phenotypic value representing the total growth
rate of the tree. This value, usually dubbed the \'site index\' of the
tree, represents the expected (average) height, in feet, the tree will
attain by age 25.

Following the same procedure, the code computes a phenotypic value for
RMAX by adding a random environmental effect to the genotypic value
generated in Section B.2, assuming hsqr0 = 0.6 (Roberds et al. 2003). As
the stand ages, the phenotypic value of RMAX is recalibrated each year,
via subroutine **BEGYR.**

**B.4 Edge Effects**

To diminish edge effects associated with small stands, subroutine
**DEPLOY** augments the stand with 4 additional rows above and below and
4 columns on each side. For these border trees, RMAX and g0 values are
assigned in ways that match them with the trees situated on the edges of
the stand proper. In particular, if parameters for the trees in the
stand are generated randomly and independently from fixed distributions,
the same procedure is used for the added trees. If trees in the stand
are arranged in clonal blocks, each block situated on an edge of the
stand is extended perpendicular to that edge through 4 more rows or
columns of trees having the same pair of g0 and RMAX values assigned to
trees in the block. Corners are extended so the rectangular shape of the
original stand is preserved. The augmented stand is then simulated.
However, all reported code outputs relating to trees (percent survival,

numbers killed by beetles, total merchantable volume at harvest, etc.)
or to beetle populations (average numbers of penetrations per attacked
tree, average brood success, etc.) pertain only to the *original* stand.

##### 

##### 

##### C. Parameters Associated with Beetle Population Dynamics

Once the stand is planted, most of the code computation is related to
changes that occur within beetle populations, and the effects that these
changes have on trees in the stand.

#####  

##### 

##### C.1 Computation of SPB Vital Rates - Subroutine PARBTL

#####  While SPB vital rates are influenced by a variety of external factors, the most important of these is temperature (Payne 1980; Coulson 1980; Wagner et al. 1983/1984; Wagner et al. 1984; Stephen and Lih 1985; Ungerer et al. 1999). In SPBLOB simulations, prior to commencement of yearly and daily updates of beetle populations, subroutine PARBTL computes these rates as functions of temperature, using formulas described in Wagner et al. (1983/84) and Ungerer et al.(1999). These temperature-dependent functions fall into two categories. The first computes the per-day parent female oviposition rate, as a *proportion* of the average total number of eggs a female parent deposits while in a tree, *if* temperature were to remain constant throughout the day. Formulas in the second category compute the *proportion* of total development within each juvenile life stage (egg, larva, pupa, preadult) that would accrue over the course of a day, again assuming temperature remains constant throughout the day. For each temperature from 1^o^C to 45^o^C, values of these five rates are stored in matrix DEVRAT, to be used throughout the remainder of the simulation (in particular, see Sections G and K, below.) 

**C.2 Beetle Immigration and Emigration**

To some extent, most loblolly plantations have low-level, endemic beetle
populations, even in absence of infestations that kill trees. These
beetles can result from chance events such as wind-blown arrival from
infested stands, but appear primarily to be associated with weakened and
lightning-struck trees within the stand (see, e.g., Blanche et al. 1983
and references therein). Occasionally, favorable weather or high
immigration will initiate stand infestation consisting of enough beetles
to kill healthy trees, the progeny spreading to nearby trees and
creating the \'spots\' of dead trees familiar in infested stands.

**C.2.1 Immigration**

In the southeastern United States, beetle populations build, on average,
about once in 10-12 years. Based on this observation, **SPBLOB** models
the immigration process as a Markov chain (Bishir et al. 2004) having a
probability of approximately 1/10 of occurring in any given year. In
years when immigration occurs, the per day rate of beetle arrival to the
stand is modeled as

avimmi = parimmi\*ntrees\*exp(-(y\*\*2)/2)

where

parimmi = maximum arrival rate per tree per day,

ntrees = number of live trees in the stand,

y = (nday - 110)/50

and

nday = day number during the year (January 1 is day 1, etc.).

The default value of parimmi is 0.04, chosen to approximate 'spot'
activity observed in the field, and leading to an average total
immigration of about 4,200 adult females during migration years, into a
stand containing 1000 trees. (As is the case for all parameters, this
value - assigned in subroutine **SIMUL** - can be modified by the user.)

**C.2.2** **Emigration**

Beetle departure from the stand is possible whenever SPB are present in
the stand, not just during immigration years, and is modeled in a
similar way to immigration. Here, the probability that an *individual*
beetle leaves the stand during a particular update period (these occur 4
times per day) is

premi = paremi\* exp(-(y\*\*2)/2)/(1.0 + (yrnum/15.0)\*\*4)

where

paremi = maximum departure rate per beetle, per update period,

yrnum = stand age in years,

y = (nday - 110)/50

and

nday = day number during the year (January 1 is day 1, etc.).

The default value of paremi is 0.001, leading to a maximum departure
probability, per update period, of 0.00088 (about 1 beetle in 1,130) in
year 9, declining to 0.000033 (about 1 beetle in 30,000) per update
period in year 35. At these rates, the probability a beetle in a 9 year
old stand would leave *sometime* during the year is about 0.5, a figure
that drops to about 0.02 at stand age 35. The exit probability decreases
as the stand matures because beetles appear to prefer taller, denser and
older stands of loblolly pine (Coster and Searcy 1981; Belanger et al.
1993).

Beetle migration itself is computed in subroutine **DISPERSE**,
departing beetles being a portion of the dispersing group, while
immigrants are added to those beetles that will constitute the
\'landed\' category in the next update period. See Section M.

**D. Simulation of Temperatures**

In the southeastern and southern United States, the region of interest
in this study, daily long-term average temperatures over the course of
the year roughly follow a sine curve, the lowest point occurring about
January 20 and the highest about July 20. Such a curve is completely
determined by TMEAN, the overall yearly average, and TDEV, the
difference between the average for July 20 and the yearly average,
TMEAN. Typically, both these parameters are set by the user. For any
particular location, these two quantities can be found by going to
www.noaa.gov and clicking on NOAA Climate Data at a Glance. Once these
parameters are specified, the long-term average for any particular day
can be computed during the course of each simulation.

Actual daily averages can vary considerably from the corresponding
long-term averages. Also, successive daily readings tend to be
positively correlated (e.g., if one day is warmer than average, the next
tends to be, also). In the **SPBLOB** code, these aspects of daily
temperatures are modeled in two steps:

1\) (Correlation) av(n) = \[av(n-1) + long(n)\]/2,

where

av(n) = \'expected\' average temperature on day n

av(n -1) = average observed temperature on day n-1

and

long(n) = long-term average temperature for day n.

2\) (Randomization) The simulated \'actual\' average, high, and low

temperatures on day n are generated using normal distributions

centered at av(n).

Once the average and high temperatures are known for the day,
temperatures over the course of that day can then be approximated by
another sine curve (Baskerville and Emin, 1969). In **SPBLOB**, we use
subroutine **TEMP** to compute these temperatures every three hours,
then read corresponding temperature-dependent SPB rates from matrix
DEVRAT (Section C.1, above) for each three hour period, assuming a
constant temperature level over the period.

##### 

##### E. Frequency of Lightning Strikes 

#####  Lightning plays an important role in the loblolly pine/southern pine beetle system. For instance, Hodges and Pickard (1971) estimated that 31% of the spots that occur in loblolly stands originate in lightning-struck trees. In SPBLOB we follow Coulson et al. (1983) in modeling the average number of thunderstorms on each day during the year as a multiple of a normal density having mean 220 and standard deviation 60. Since day \# 220 is August 8, this means that storm activity is highest in early August, although storms are common over most of the summer and fall seasons. The multiple depends on the geographic location, being highest in Florida and diminishing as we move north (Coulson et al. 1983, Fig. 2). A multiple of 60, e.g., results in an average of about 50 storms per year. Assuming 2 strikes per square mile, per storm (Coulson et al. 1983, p. 186) we estimate that a stand planted using 3m spacing and having 32 rows and columns (about 2.28 acres, or 0.9 ha) would experience about 12 strikes between ages 9 and 35 years, the default simulation period in SPBLOB. More extensive computational details appear in the file spblob_params_lightning.doc, contained in the Supporting Information folder.

##### 

**F. Beetle Mortality Due to Extremely Cold Temperatures**

Following Ungerer et al. (1999) beetle mortality engendered by cold
temperatures are modeled in subroutine **COLD** using normal
distributions. Suppose X is a normal random variable with mean μ and
standard deviation σ.  In model simulations, if a day's simulated
minimum temperature, tmin, is less than --2^o^C, then

-   Adult landed and parent beetles (see Section II.E.4) independently
    survive

> with probability Pr(X \< tmin), where μ = -12.9^o^C and σ = 2.9^o^C

-   Eggs and larvae survive independently with probability Pr(X \<
    tmin),

where μ = -11.5^o^C and σ = 2.91^o^C

-   Pupae survive independently with probability Pr(X \< tmin), where

μ = -9.8^o^C and σ = 2.6^o^C

Thus, for example, a single overnight temperature of -15 ^o^C kills
about 77% of adults, 88% of larvae and almost 98% of the pupae in each
tree, mortality levels which can profoundly decrease potential tree
damage in the following year.

**G. Juvenile Development and Mortality over the Course of a Day**

In addition to the special beetle mortality discussed in Section F,
there is mortality associated with all beetle activities. The remainder
of this section discusses juvenile development as well as mortality,
while adult mortalities are the subjects of Sections K.3 and M below.

For each juvenile life stage, subroutine **UPDATE** starts with the
temperatures generated every three hours by subroutine **TEMP**, reads
from array DEVRAT (Sections C.1 and D, above) - for that life stage -
the per-day development proportion corresponding to each of these 8
temperatures, multiplies each proportion by 1/8 to convert to a 3-hour
basis, and sums values for the eight proportions to obtain the
***proportion*** of total development in that life stage during the day.
Subroutine **ENDOFDAY** (Sections J and K below) then checks the current
life stage of each juvenile cohort and adds today's additional
proportion to the proportion already accumulated. If the new proportion
is less than 1.0, the cohort remains in the same stage. Otherwise, the
code moves the cohort to the next stage, and assigns its state of
maturity in the new stage as the accumulated proportion in the preceding
stage, less 1.0.

The code then assesses the proportion of beetle natural mortality in
each stage, over and above that due to cold temperatures. Following
Coulson et al. (1977), Fargo et al. (1979) and Coulson (1980, Figs. 5-8,
5-9), it applies a constant survival rate of 0.98 *per 1% of **total**
(egg to maturity) development*, a figure chosen to result in a number of
brood adults totaling about 13% of the original number of new eggs. In
addition, based on Thatcher and Pickard (1967), additional larval
mortality is assigned on unusually hot days. Details concerning both
these computations appear in the file
**spblob_params_SPB_natural_mort.doc** contained in the **Supporting
Information** folder.

**H. The Process of Beetle Attack on Trees - Subroutine ATTACK**

To reproduce, female SPB must successfully penetrate the bark of a tree,
create egg galleries in the cambium, and oviposit. Such attacked trees
must be alive at the time of attack and their ability to repel an attack
must have been overcome. The code models this reproductive process in
successive steps:

-   If a live tree has been struck by lightning, the probability that an
    arriving beetle

attacks is 1.0 (that is, all landed beetles attack)

-   Otherwise, the probability a landed beetle attacks is equal to the
    larger of

p1 = (0.37 - 0.345\*exp(-RATK6\*phrpos))\*exp(-RATK3\*phrneg)

and

R = (npen/bole)/100,

where

> exp denotes the exponential function,

phrpos = concentration of positive pheromones near the tree

> phrneg = concentration of negative pheromones near the tree
>
> npen = current total number of penetrations of the tree, and
>
> bole = tree bole area, in square meters.

Because successful beetle reproduction depends on tree death, this form
for R is based on an assumption that beetles are more likely to attack
as tree death approaches. The values of parameters RATK6 AND RATK3 are
set in subroutine **ATTACK**, being chosen so average levels of tree
mortality seen in code simulations seem reasonable. As more information
becomes available, it can be incorporated as desired.

The probability that an attacking beetle successfully penetrates to the
cambium depends on the tree\'s resin flow rate, and is modeled as:

prob = 0.94/(1.0 + (resin(i,j)/ratk50)^m^),

the values ratk50 = 2.3 and m = 2 being chosen to reflect experimental
observations, primarily Hodges et al. (1979)**.** Details are in the
file **spblob_params_SPB_attack_success.doc** contained in the
**Supporting Information** folder. The numerator value 0.94 derives from
the observation in Bunt et al. (1980) that 6% of beetles landing on a
tree during mass attack succumbed to predators.

##### J. Tree Response to Beetle Attack - Subroutine ENDOFDAY

Trees defend themselves from SPB primarily by exuding resin at attack
sites on the bole. Trees with sufficient resin 'pitch out' beetles
attempting to enter. Those with naturally low resin flow rates, and
others that lose resin due to lightning strikes or drought, are more
susceptible to attack and often serve as initiation points for spot
formation (Gara et al. 1965; Hodges and Pickard 1971).

##### J.1 Resin Loss from Flow Through Attack Holes

In the code, trees that are attacked by beetles but have not been struck
by lightning lose resin at a rate proportional to both the number of
penetration holes and the tree\'s resin flow rate, according to the
formula

new_resin_level = previous_level\*(1.0 - 0.02\*natk/bole)

where

natk = total number of penetration holes through which resin is flowing

and

bole = total bole area of the tree.

The multiplier 0.02 corresponds to an average resin flow rate per hole.
Resistant trees respond rapidly to this loss by producing new resin
(Ruel et al. 1998). However, any tree can be overcome by a sufficiently
large number of beetle attacks.

Lightning struck trees are modeled differently, in that there is no
resin replenishment. Following Hodges and Pickard (1971), who reported
flow rates in struck trees being only 5-15% of those in similar control
trees, we model tree resin level as

new_resin_level = 0.905\*previous_level.

The multiplier 0.905 leads to a 95% loss of tree resin within 30 days,
leaving these trees highly susceptible to attack.

**J.2 Tree Mortality Due to Beetle Penetrations**

Based on Hodges et al. (1979) we assume that, on average, the number of
successful beetle penetrations, per m^2^ of bole area, required to kill
a tree is about 100\*r0, where r0 is the tree\'s resin flow rate. (For
explanation of this formula, see the Supporting File entitled
**spblob_params_SPB_attack_success.doc)**. In the code, the *actual*
number of penetrations needed in any particular tree is modeled as a
normal random variable having mean mu = 100\*r0 and standard deviation
sig = 10\*r0. Thus, for example, the probability that a tree having
resin flow value r0 = 1.0 dies from 110 or fewer penetrations equals
Pr(W ≤ 110) = 0.8413.

Daily code computations appear in subroutine **UPDATE**, and proceed as
in this

**Example**:

Suppose that in preceding days, a surviving tree T - with inherent resin
flow value

1.0 - has already suffered a total of x successful penetrations per unit
bole area.

If by today's end the new total is y, then

> Pr(T dies today) = Pr(W ≤ yW \> x) = Pr(x \< W ≤ y)/Pr(W \> x).

For instance, if x = 100 and y = 110, then

Pr(T dies today) = Pr(100 \< W ≤ 110)/Pr(W \> 100)

= (0.8413 - 0.5)/0.5 = 0.6826.

For a more general discussion, go to the **Supporting Information**
folder and consult the file **spblob_params_tree_mort_from_SPB.doc**.

##### 

##### 

##### K. Activities of Parent Beetles - Subroutine ENDOFDAY (continued)

#####  

##### In SPBLOB, all parent beetle activities are directly or indirectly based on parameter values drawn from the literature (Wagner et al 1983/1984; Ungerer et al. 1999). 

##### 

##### K.1. *Proportion* of Her Total Eggs That a Female Parent Deposits in a Day

##### 

#####  Starting with the temperatures generated every three hours by subroutine TEMP, subroutine EGG reads from array DEVRAT (Section C.1) the *daily* oviposition proportion corresponding to each temperature, multiplies this proportion by 1/8 to convert to a 3-hour basis, and sums values for the eight proportions to obtain eggpro, the *proportion* of the total egg production a parent female will deposit that day.

**K.2 Total *Number* (newegg) of Eggs Laid in Each Tree**

Coulson et al. (1976, Fig. 3) model the total egg gallery length that a
single female parent creates during the course of her stay in a
particular tree as a function of the density of parents presently in the
tree. In symbols,

G = 66.41exp(-0.245A),

where

G = average total gallery length (cm) created per female parent,

A = average number of female parents per m^2^ of bole area,

and

exp denotes the exponential function.

Using the Foltz et al. (1976) estimate that about 1.59 eggs are
deposited per cm of gallery length, and replacing A by B = average
number of female parents per *square ft* of bole area, the formula for G
translates into

TE = 105.6\*exp(-0.0264B),

where

TE = average total number of eggs a female will deposit in the tree.

Since the code tracks only the *female* portion of the beetle
population, we assume a 1:1 sex ratio (e.g., Coulson et al. 1976; Bunt
et al. 1980) and approximate the total number of female eggs that
ultimately will be deposited in the tree by ALL the parents now in the
tree as

nnnegg = int(0.5\*PAR\*TE),

where

PAR = total number of parents in the tree

and the function int returns the largest integer less than or equal to
its argument.

Finally, on a particular day during which the egg proportion **eggpro**
(Section K.1) is deposited, the total number of new eggs deposited in
the tree is

newegg = bin(nnnegg, eggpro)

In words, newegg is a randomly generated value from a binomial
distribution (a \'binomial deviate\') with nnnegg \'trials\' and
per-trial \'success\' probability equal to eggpro.

**K.3. Daily Parent Survival and Reemergence**

In apparent absence of in-vivo observations of the SPB parent survival
rate, the **SPBLOB** code approximates this rate indirectly, based on
three facts:

-   When a parent completes her oviposition in a tree, the cumulative
    proportion of her total egg complement deposited is 1.0

-   Estimates of reemergence rates of parent SPB vary widely (Coulson
    1980, p. 77), the average of these estimates being about 0.7

-   Algebraically, exponents are added when products of exponential
    quantities are computed: e.g., a^x^a^y^ = a^x+y^

Following the observation in Gagne et al. (1982, p. 1220) that \"the
pattern of oviposition rate \... is surprisingly similar to the pattern
of reemergence,\" we compute the proportion of parents that reemerge on
a particular day as *equal* to eggpro, the proportion of her total eggs
a female deposits that day (Section K.1). (This observation, and our
application of it, are consistent with the supposition that parents
reemerge *because* their oviposition is complete. As a consequence, the
total number of reemerging parents is computed as

RE = bin(PAR, eggpro),

where, as in Section K.2,

PAR is the total number of parents in the tree,

and

bin denotes a binomial deviate.

To ensure a reemergence rate consistent with the literature, we begin
with the fact that the total probability a parent survives her stay in a
tree is the product of the daily survival probabilities encountered
along the way. Based on the three \'facts\' listed above, we therefore
model the probability that a parent survives on a particular day as

ps = 0.7^eggpro^.

The *definition* of eggpro (Section K.1) implies that a surviving parent
should reemerge when her eggpro values total 1.0. Thus, the formula for
ps, together with 'fact' \# 3 above, produces a total survival
probability (equivalently, a reemergence probability) of 0.7^1.0^ = 0.7,
consistent with the second \'fact\' above.

Finally, we update the parent population itself using the formulas

NEWPAR = PAR - RE + PEN

and

NEWPAR = bin(NEWPAR, ps),

where

PAR, RE, ps and bin are as above,

NEWPAR is the new number of parents,

and

PEN is the total number of successful penetrations by landed beetles
today.

##### L. Pheromone Production, Concentrations and Diffusion 

Female beetles that successfully penetrate tree bark emit pheromones
that attract both sexes; this pheromone production ceases when bark
penetration is complete (Payne, 1980). In the model, each female that
enters is joined by a male, and the two work to create galleries along
which the female deposits eggs. Males that join successful females in a
penetration hole emit pheromones that repel other beetles, an emission
that continues throughout gallery construction and oviposition. The
result is an attack process that initially increases in intensity, as
more beetles are attracted to female pheromones, then diminishes as
increasing numbers of males arrive (Payne, 1980; Coulson, 1980).

**L.1 Production - Subroutines UPDATE and ENDOFDAY**

In the code, positive and negative pheromone concentrations are updated
each day. The formulas, which apply to each live tree, are:

newpos = RPHER1\*PEN \[in **UPDATE**, 4 times per day\]

newneg = RPHER2\*nmale \[in **ENDOFDAY**\]

where, for each particular tree,

PEN = the number of attacking females that have successfully penetrated
the

tree\'s bark \'today\'

nmale = the number of males (cumulative, up to the present) that have

joined successful females in attack holes

newpos = the positive pheromone concentration resulting from female

production

newneg = the negative pheromone concentration resulting from male

production

and

RPHER1 and RPHER2 are the per individual productions from females

and males, respectively.

Because in the model beetle movement depends only on *relative*
pheromone concentrations on nearby trees (Section M, below), values of
parameters RPHER1 and RPHER2 are arbitrarily set equal to 1.0. (These
parameters thus help clarify the structure of the model, but do not
represent actual pheromone levels observed in field studies.)

**L.2 Diffusion and Resulting Concentrations - Subroutine DIFFUSE**

The pheromone production streams newpos and newneg diffuse rapidly (in
the model, within the same update period) through the stand, leading to
adjusted concentrations that affect beetle attack and dispersal
processes. The *concentration* of positive pheromone surrounding tree T
is computed as

pherpos(T) = newpos(T) + Σ newpos(B)\*(((0.9)\*\*m)/(8\*m)),

where

pherpos(T) = total concentration around tree T

newpos(B) = positive pheromone produced at tree B (Section L.1)

and

m = larger of the number of rows and the number of columns separating
the

two trees.

The sum extends over all trees in the stand, with the exception of tree
T itself. The \'dissipation\' rate 0.9 approximates measurements
reported in Thistle et al. (2004) and discussed above in Section
II.E.5.2, while 8m is the number of trees located in the m-th ring of
trees surrounding tree T. \[A similar formula is used to compute the
concentration of negative pheromone. See the **Source Code** for
details.\]

**M. Beetle Dispersal - Subroutine DISPERSE**

In the model, beetle dispersal from a tree is governed by the *relative*
pheromone concentrations associated with trees that are adjacent to the
source tree (Section L.2). More precisely, all beetles that depart from
a tree (T) move to the tree adjacent to T that has the largest positive
pheromone concentration.

Based on Coulson (1980) and Pope et al. (1980), the temperature
dependent mortality rate, per day, experienced by dispersing beetles is
modeled in subroutine **UPDATE** as

DFLY = 0.03748 + 0.00562\*exp(0.11644\*avtemp),

where

avtemp = average temperature for the day.

See the file **spblob_params_dispersal_mort.doc** in the **Supporting
Information** folder for additional details.

**N. Tree Decay and Fall - Subroutine SIMUL**

In the code, trees that die, whether from competition, lightning or
beetle attack, are categorized as \'fallen\' after being dead 5 years.
The impact is to shift over to standing trees the beetle immigrants that
otherwise would land on dead trees.

**P. Merchantable Volume - Subroutine OUTPUT**

Beginning at stand age 9 years, and continuing until tree harvest,
subroutine **OUTPUT** computes the total cubic foot volume of
merchantable timber available in the stand. While the original version
of this **PTAEDA** subroutine (Daniels and Burkhart (1975); Burkhart et
al. 1987) offers several choices for the meaning of \'merchantable\',
**SPBLOB** reports only 'inside bark' volume of the portion of the tree
bole between stump height and the upper point at which the diameter =
4", taken from all live trees whose diameters exceed 4.5 inches (in
**PTAEDA** notation, dbh \> 5.0). However, these default choices can be
altered by the user, if desired.

# **VI. GLOSSARY FOR THE SPBLOB SOURCE CODE** 

#  In this Glossary of variables and parameters used in the **SPBLOB** code, **PTAEDA2** is the Burkhart et al. (1975, 1987) tree competition and growth code, portions 

# of which are incorporated by permission into **SPBLOB**. In this Glossary:

**Variables** are quantities that change over time within a particular
simulation

**Parameters** are quantities that do not change value within a
particular simulation

of a particular stand

# 

# **Constants** are quantities whose values remain the same in all model simulations

# 

# Most of these quantities are numerical, but a few have logical or character type. 

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

**acres** The number of acres in the simulated stand.

**alpha** First parameter of a gamma distribution for resin flow that is

correlated with a previously generated g0 value (see **beta**

below).

**alphares** First parameter of the gamma distribution for tree resin
flow rates

reported in Roberds et al. (2003) (see **betares**).

**amerch** Average merchantable timber volume per live tree.

**ampl** The deviation, on a particular day, between the simulated

'actual' average temperature for that day, and the simulated

high (or low) temperature for the day.

**avdeath** Average number of penetrations, per m^2^ of bole area,
needed to kill

a tree. The number needed is assumed to be a N(100,100)

random variable.

**avdiam** A vector, in which the k-th entry is the average diameter
(in) of all live

trees of clone type k 'this' year.

**avdiamkb** A vector, in which the k-th entry is the average diameter
(in) of all

trees of clone type k that were killed by beetles 'this' year.

**avdiamkc** A vector, in which the k-th entry is the average diameter
(in) of all

trees of clone type k that were killed by competition 'this' year.

**avhgt** A vector, in which the k-th entry is the average height (ft)
of all live

trees of clone type k 'this' year.

**avhgtkb** A vector, in which the k-th entry is the average height (ft)
of all

trees of clone type k that were killed by beetles 'this' year.

**avhgtkc** A vector, in which the k-th entry is the average height (ft)
of all

trees of clone type k that were killed by competition 'this'

year.

**avimmi** The average number of beetles, per (alive or dead) standing
tree, that

migrate into the stand on a particular day.

**avres** A vector, in which the k-th entry is the average resin flow
rate

(gm/24h) of all live trees of clone type k 'this' year.

**avreskb** A vector, in which the k-th entry is the average resin flow
rate of all

trees of clone type k that were killed by beetles 'this' year.

**avstrk** Average number of lightning strikes in the stand on a
particular day.

**avtemp** The simulated 'actual' average temperature on a particular
day.

**basal** A vector, in which the k-th entry is the total basal area
(ft^2^/ac) of all

live trees of clone type k 'this' year.

**batch** A variable which has the value zero (0) when the Source Code
is

being used to generate the 'Interactive' user code, and the value one

\(1\) when generating the 'Multiple-Simulations' code.

**beetles** A logical variable which is true if there are SPB in the
stand 'this' year,

and otherwise false.

**beta** Second parameter of a gamma distribution for resin flow that is

correlated with a previously generated g0 value (see **alpha**).

**betares** Second parameter of the gamma distribution for tree resin
flow rates

reported in Roberds et al. (2003) (see **alphares**).

**bole** A matrix, in which the i,j entry is the total bole surface area
(ft^2^) of the

tree in the i,j position in the stand.

**bspc** In simulated stands involving clonal blocks, **bspc** is the
number

of columns of blocks in the stand.

**bspr** In simulated stands involving clonal blocks, **bspr** is the
number

of rows of blocks in the stand.

**cenG0MEAN** The mean of the censored normal distribution from which
clonal g0

values are obtained.

**cenG0VAR** The variance of the censored normal distribution from which
clonal g0

values are obtained.

**cip** A matrix, the i,j entry of which is the 'pine competition index'

associated with the tree in the i-th row and j-th column of the

stand (see Burkhart et al. 1987).

**clonetyp** A matrix, in which the i,j entry is the clone number of the
tree in the

i-th row and the j-th column of the stand (the value is always zero if

clones are not used to plant the stand).

**clotyp** A matrix, in which the i,j entry is the number of the clone
in the clonal

block located at the intersection of the i-th row of blocks and the j-th

column of blocks in the stand (the value is always zero if clones are

not used to plant the stand; if there are clones, but no blocks, then

each tree constitutes a 'block').

**colda** Survival rate for landed and parent beetles as a function of
low winter

temperatures.

**coldel** Survival rate for eggs and larvae as a function of low winter

temperatures.

**coldp** Survival rate for pupae as a function of low winter
temperatures.

**crnlgt** A matrix, in which the i,j entry is the crown length
associated with the

tree in the i-th row and j-th column of the stand.

**cspb** In simulated stands involving clonal blocks, the number of
columns

of trees in each block (see **rspb**).

**daydied** A matrix, in which the i,j entry is the calendar day (1 =
Jan 1, etc) on

which tree i,j died, if it was killed by beetles or by lightning during
the

course of a simulation. (If it did not die, or if it died because of
intra-

stand competition, the entry is 0.) This matrix, together with matrix

**yrdied**, below, enables us to track the course of beetle-caused
deaths

during each stand year, from age 9 until harvest.

**depart** A matrix, in which the i,j entry is the total number of
beetles leaving

tree i,j at the end of one of the 4 per day update computations.

**devpar** A matrix containing parameter values from Ungerer et al.
(1999) that

govern temperature related rates of beetle oviposition, and of

maturation rates in the egg, larval, pupal and teneral adult stages

(see **devrat**).

**devrat** A matrix containing temperature dependent rates of beetle

oviposition, and of maturation rates in the egg, larval, pupal

and teneral adult stages (see **devpar**).

**DFLY** Probability of beetle survival during dispersal -- depends on

temperature.

**diam** A matrix, in which the i,j entry is the bole diameter (in) of
the tree

in the i,j position in the stand.

**diff** In the censored normal distribution used for clonal g0 (tree
growth)

rates, the difference between the smallest acceptable g0 and

the mean value of the original distribution. See **ming0**.

**distan** In **PTAEDA2** (Burkhart et al. 1987)**,** the distance
between two

competing trees.

**dvelop** A vector, the k-th entry of which corresponds to the k-th
juvenile life

stage (egg, larva, pupa, teneral adult) and represents the

proportion of total development within that stage accrued by

juveniles in each cohort that is currently in that stage.

**dvjuv2** A 3-D matrix which, together with matrices **mjuv1** and
**njuv3**,

contains information relating to all juvenile beetles in the stand.

The i,j,k entry of dvjuv2 is the proportion of total development within

its current developmental class attained thus far by the juvenile cohort

oviposited k days earlier in tree i,j.

**eggpro** The temperature dependent proportion of her total complement
of

eggs each parent adult produces on a particular day. Used to

calculate both the number of eggs each female parent oviposits on a

particular day, and, following Coulson et al. (1978) and Gagne et al.

(1982), the number of female parent beetles that reemerge that day.

**factor** A variable that depends on stand age -- For each surviving
tree in the

stand at the beginning of 'this' year, the code computes its current

maximum resin flow rate as

RMAX = factor\*(the maximum rate at time of stand planting).

**flow** For each tree, a temporary variable equal to today's total
amount of

resin lost through beetle attack holes.

**g0** A matrix, in which the i,j entry is the overall growth rate
associated

with the tree in the ith row and jth column of the simulated stand.

('Overall' rate = genetic rate, plus an added environmental amount,

which can be positive or negative, associated with stand location i,j.)

**g0gain** The difference between the mean of the censored normal
distribution

used for clonal g0 (tree growth) rates, and the mean of the original

normal distribution.

**G0MEAN** The mean of the normal distribution that is sampled to obtain
tree

growth rates when a stand is planted with improved tree

materials. When clonal ramets are used in stand planting,

this same distribution is sampled and the values censored to

generate the clonal g0 values (see **gc,** below).

**G0VAR** The variance of the underlying normal distribution that is
sampled to

obtain tree growth rates (see **G0MEAN** above).

**gc** A matrix, in which the i,j entry is the genetic component of the

growth rate for the clone planted in the clonal block located

at the intersection of the i-th row of blocks and the j-th

column of blocks in the stand (see **rc**).

**gengam** Random deviate from a gamma distribution having specified
alpha

and beta parameter values.

**gennor** Random deviate from a normal distribution having specified
mean

and variance.

**genunf** Random deviate from a uniform distribution on a specified
interval.

**gp** A temporary variable used in subroutine **DEPLOY** to fill a
clonal block

with trees that all have the same *genetic* value of growth rate g0.

**hd** A matrix, in which the entries are computed initially in
subroutine

**JUVEN**. At the beginning of each subsequent year, the

preceding year's values are used in **GROW2** to calculate the

new height increments for all surviving trees, after which **hd**

is updated for use in the following year.

**height** A matrix, in which the i,j entry is the current height of
tree i,j.

**hsqg0** The heritability (broad-sense) associated with tree growth
rates.

(Used in computing the entries in matrix **g0**, above.)

**hsqr0** The heritability (broad-sense) associated with tree resin flow
rates.

**ignbin** Random deviate from a binomial distribution having specified
number

of 'trials' and 'success' probability.

**ignpoi** Random deviate from a Poisson distribution having a specified
mean.

**immi** A 'character' used in Interactive simulations to indicate on
the screen

whether there is immigration in each particular year.

(immi = 'true' if there is immigration 'that' year, and

immi = 'false' if no immigration occurs that year).

**killb** A vector, the k-th entry of which counts, for any particular
year, the

number of trees of clone number k killed by beetles that year.

**killc** A vector, the k-th entry of which counts, for any particular
year, the

number of trees of clone number k killed by competition that year.

**killl** A vector, the k-th entry of which counts, for any particular
year, the

number of trees of clone number k killed by lightning that year.

**landed** A matrix, in which the i,j entry is the current number of
beetles on the

bole surface of tree i,j. Recomputed 4 times per simulated day.

**light** A matrix, in which the i,j entry is the current number of days
since tree

i,j was struck by lightning (the value = zero if the tree has never been

struck).

**lotemp1** A variable equal to the lowest temperature encountered so
far in a

particular year. Reset to zero at the beginning of each year, then

updated daily within the year.

**lotemp2** A variable equal to the second lowest temperature
encountered so far

in a particular year. Reset to zero at the beginning of each year, then

updated daily.

**lotemp3** A variable equal to the third lowest temperature encountered
so far in

a particular year. Reset to zero at the beginning of each year, then

updated daily.

**maxday** A numerical constant equal to 469, and representing the last

day in a simulated model year (see **minday**, below).

**merch** A vector, the k-th entry of which is the total merchantable
timber

volume (ft^3^/ac) in surviving trees of clone type k (see **vmerch)**.

**minday** A numerical constant representing the first day in a
simulated

model year. Set equal to 105 in order to begin the \'season\'

after the last cold winter temperatures (value will change

depending on the geographical location of the stand being

simulated - see **maxday**, above).

**ming0** In the censored normal distribution used for clonal g0 (tree
growth)

rates, the smallest acceptable g0 value.

**mjuv1** A 3-D matrix which, together with matrices **dvjuv2** and
**njuv3**,

contains information relating to all juvenile beetles in the stand. The

i,j,k entry of mjuv1 is the current developmental class (egg, larva,
pupa

or teneral adult) attained by the juvenile cohort oviposited k days

earlier in tree i,j.

**mtemp** Average temperature during a 3-hour update period.

**mu** The mean value of the normal random variable representing the

number of successful beetle penetrations per m^2^ of bole

area required to kill a tree (see **sig).**

**nalive** The current total number of live trees in the stand.

**natk** A matrix, in which the i,j entry is the cumulative number of
beetle

attacks on tree i,j during a day.

**nc** The number of columns of trees in the simulated stand. This

parameter is read in **MAIN** from an input file.

**nc0** Temporary value of **nc**.

**nclone** A vector, the k-th entry of which is the total number of
trees of clone

type k planted in the stand.

**nday** During a simulation, nday is the current day of the year for
which

computations are being made.

**ndead** The total number of dead trees in the stand.

**nemrge** A temporary variable equal to the number of juveniles that
mature and

emerge from a particular tree 'today.' Added to the departing group

prior to daily dispersal.

**nenter** A matrix, the i,j entry of which is the cumulative number of
successful

beetle penetrations 'today'.

**newegg** For each tree i,j in the code computations, the name used
temporarily

for the total number of eggs oviposited in that tree on a particular
day.

The value is later stored as the i,j,1 entry in matrix **njuv3.**

**newneg** A matrix, the i,j entry of which is the current negative
(male, repelling)

pheromone concentration currently being generated in tree i,j.

**newpos** A matrix, the i,j entry of which is the current positive
(female,

attracting) pheromone concentration currently being generated on

and/or in tree i,j.

**nfallen** In a simulation, the cumulative number of 'fallen' trees
(trees that

died 5 or more years ago) in the stand.

**nimmi** A variable which equals zero if there is no beetle migration
into the

stand 'this' year, and equals 1 if there is migration.

**njuv3** A 3-D matrix which, together with matrices **mjuv1** and
**dvjuv2**,

contains information relating to all juvenile beetles in the

stand. The i,j,k entry of njuv3 is the number of survivors in

the juvenile cohort oviposited k days earlier in tree i,j.

**nmale** The current number of male SPB in a tree. Used only to compute
the

current concentration of male, repelling pheromone being produced.

**nnpen** A matrix, in which the i,j entry is the cumulative number of
beetle

penetrations in tree i,j so far this year, as of the beginning of

'today.' Used to compute the probability tree i,j dies from

beetle penetrations today. See **npen**.

**npen** A matrix, in which the i,j entry is the cumulative number of
beetle

penetrations in tree i,j so far this year, by today\'s end. Used to

compute the probability tree i,j dies from beetle penetrations today.

See **nnpen**.

**nr** The number of rows of trees in the stand. This parameter is read
in

**MAIN** from from an input file.

**nr0** Temporary value of **nr**.

**nreemr** A temporary variable whose value is the number parent beetles
that

reemerge from a particular tree 'today.'

**nsim** The number of the simulation currently being processed, when

running a group of simulations chosen from the list **input.txt.**

**nstay** A temporary variable whose value is the number parent beetles
that

do not reemerge from a particular tree 'today' (i.e., remain in the tree

at day's end).

**nstrks** The number of lightning strikes in the stand today (a Poisson
deviate).

**ntrees** The total number of trees planted in a simulated stand (equal
to

nr\*nc).

**numcln** In simulations involving clones, the number of clones used in
planting

the stand.

**numdevs** A parameter: In the censored normal distribution used for
clonal g0

(tree growth) rates, the number of standard deviations

between the mean value of the original distribution and the

smallest acceptable value of g0. See **diff** and **ming0**.

**numtre** A variable, the value of which is the total number of live
trees in a

simulated stand at the beginning of 'this year.'

**nyears** The harvest age of a simulated stand.

**oalive** A numerical constant having value +2. Used in the code to
indicate a

tree that is currently still alive.

**odead** A numerical constant having value zero. Used in the code to
indicate

a tree that died, due to beetle attacks, lightning or intrastand

competition.

**paremi** A parameter whose value is used to calibrate the probability
that a

particular beetle emigrates from the stand during a particular

update period.

**parent** A matrix, the i,j entry of which is the current number of
nesting

(ovipositing) adult females in tree i,j.

**parimmi** A parameter, used to calibrate the number of beetles that
migrate into

a stand during an update period.

**phrmax** A temporary variable, whose value is the highest positive
pheromone

concentration among the 9 trees located no more than one row or

column away from a particular tree of interest.

**phrneg** A matrix, the i,j entry of which is the present concentration
of

repelling (male) pheromone (from all sources) surrounding tree i,j.

**phrpos** A matrix, the i,j entry of which is the present concentration
of

attracting (female) pheromone (from all sources) surrounding tree i,j.

**PI** A constant, the value of which is 3.141592653585.

**pra** A temporary variable, used in computing the probability a
particular

tree dies from beetle penetrations 'today'.

**prb** A temporary variable, used in computing the probability a
particular

tree dies from beetle penetrations 'today'.

**premi** For each individual 'dispersing' beetle, the current
probability that that

beetle migrates out of the stand during an update period.

**propor** The proportion of beetle attacks on a tree that lead to
successful

penetration.

**px** The distance (m) between adjacent columns in a simulated stand.

**px0** The assigned value of variable **px**: read from an input file
at the

beginning of a simulation.

**py** The distance (meters) between adjacent rows in a simulated stand.

**py0** The assigned value of variable **py**: read from an input file
at the

beginning of a simulation.

**qblock** A parameter, equal to 1 if a simulated stand involves clonal
blocks,

and 2 if there are no blocks.

**qclone** A parameter, equal to 1 if a simulated stand involves clones,
and zero

if there are no clones.

**qimmi** A logical variable, equal to 1 if there is beetle migration
into the

stand 'this' year, and zero if there is no immigration this year.

**qspb** A logical parameter, equal to 1 if a simulated stand involves
beetles,

and zero if there are no beetles.

**r0** A matrix, the i,j entry of which is the resin flow of the tree in
the ith row

and jth column of the stand.

**r0g0mean** The \'correlated mean\' of resin flow rate r0, given a
previously

determined value of tree growth rate g0.

**r0g0var** The \'correlated variance\' of resin flow rate r0, given a
previously

determined value of tree growth rate g0.

**r0mean** The mean of the underlying distribution from which simulated
resin

flows are randomly generated (see Roberds et al. 2003).

**r0var** The variance of the underlying distribution from which
simulated resin

flows are randomly generated (see **r0mean**).

**RATK3** A parameter that relates to the relative decrease in rate of
transition

of landed beetles (see **landed**) to the attacking stage, per

unit increase in **phrneg** (repelling pheromone).

**ratk50** A parameter, used to calibrate the proportion of beetle
attacks that

result in successful penetration of a tree.

**RATK6** A parameter that relates to the relative increase in rate of
transition of

landed beetles (see l**anded**) to the attacking stage, per unit
increase in

**phrpos** (attracting pheromone level).

**rbeg** The first line chosen from file **input.txt** when simulating a
\'batch\' of

stands (see Section IV of this Documentation - also see **rend**).

**rc** A matrix, the i,j entry of which is the genetic component of the
resin

flow rate associated with the clone planted in the clonal block located

at the intersection of the i-th row of blocks and the j-th column of

blocks in the stand (see **gc**).

**rend** The last line chosen from file **input.txt** when simulating a
\'batch\' of

stands (see Section IV of this Documentation - also see **rbeg**).

**resin** A matrix, the i,j entry of which is the basic resin level for
tree i,j,

expressed as flow per day, per currently open beetle penetration hole.

In model simulations, holes remain open for only one day.

**rho** The coefficient of correlation between same-tree values of the
tree

growth rate g0 and the tree resin flow rate r0.

**RMAX** A matrix, the i,j entry of which is the maximum resin flow
rate, at the

present stand age, of the tree in the ith row and jth column of the

stand. Based on Roberds et al. (2003).

**rp** A temporary variable used to fill a clonal block with trees which
all

have the same *genetic* value of resin flow rate r0.

## 

## RPHER1 A constant, the value of which is the increase in concentration of 

##  positive pheromone (see newpos above) around a tree, per 

##  successfully attacking female.

## 

## RPHER2 A constant, the value of which is the increase in concentration of 

##  negative pheromone (see newneg above) being produced in a tree, 

##  per male consort of a female that has successfully penetrated the tree.

**rspb** In simulated stands involving clonal blocks, the number of rows
of

trees in each block (see **cspb**).

**savdiam** Average diameter (in) of all live trees 'this' year.

**savdiamkb** Average diameter (in) of all trees that were killed by
beetles 'this' year.

**savdiamkc** Average diameter (in) of all trees that were killed by
competition 'this'

year.

**savhgt** Average height (ft) of all live trees 'this'year.

**savhgtkb** Average height (ft) of all trees that were killed by
beetles 'this' year.

**savhgtkc** Average height (ft) of all trees that were killed by
competition 'this'

year.

**savres** Average resin flow rate (gm/24h) of all live trees 'this'
year.

**savreskb** Average resin flow rate (gm/24h) of all trees that were
killed by

beetles 'this' year.

**sbasal** Total basal area (ft^2^/ac) of all live trees in the stand
'this' year.

**seed** The number of the line from the file **input.txt** that is used
to initiate a

particular simulation.

**seed1** In model simulations, two seeds are used to generate random
values

associated with temperature, trees and beetles, so these quantities

can be randomized within and between simulations.

**seed2** See **seed1**, above.

**setup** The number of the line from the file **config.txt** that is
used in running

a particular simulation.

**sig** The standard deviation of the normal random variable
representing

the number of successful beetle penetrations per m^2^ of bole

area required to kill a tree **(**see **mu).**

**sim** A temporary variable, the value of which is the current
simulation (line

number in the file input.txt) in a batch sequence of simulations.

**siteg0** The environmental variation added to a tree\'s genetic
g0-value during

the planting process.

**siter0** The environmental variation added to a tree\'s genetic
r0-value during

the planting process.

**sjuv** The probability of juvenile survival on a particular day, as a
function of

the amount of maturation on that day.

**skillb** A variable which, for any particular year, counts the total
number of

trees in the stand killed by beetles that year.

**skillc** A variable which, for any particular year, counts the total
number of

trees in the stand killed by competition that year.

**skilll** A variable which, for any particular year, counts the total
number of

trees in the stand killed by lightning that year.

**smerch** The total volume (ft^3^/ac) of merchantable timber in
surviving trees in

the stand (see **vmerch**).

**ssurv** A variable which, for any particular year, counts the total
number of

surviving trees in the stand at the end of that year.

**staytokill** A logical variable which, if equal to \'true,\' indicates
that beetles don\'t

leave when tree death is imminent. In this circumstance, the

proportion of landed beetles that attack on a particular day equals the

ratio of variable **npen** (above) to the \'expected\' total number of

beetle penetrations that will ultimately be required to kill the tree.

**surv** A vector, the kth element of which is the number of surviving
ramets

of clone k.

**surv1** The probability a juvenile beetle is alive at the beginning of
a day, as

a function of its level of maturation.

**surv2** The probability a juvenile beetle is alive at the end of a
day, as a

function of its level of maturation.

**survd** The probability a beetle that departs from one tree survives
to land

on another, not necessarily different, tree in the stand.

**survp** The probability a parent beetle survives throughout a
particular day,

as a function of the proportion of her total egg production

deposited that day.

**t0** A matrix, the i,j entry of which is the resin toxicity of the
tree in the

ith row and jth column of the stand. This parameter is not used in

current (1/07) source code computations.

**T0MEAN** The mean value of the resin toxicity distribution. Not used
in current

(1/07) source code computations.

**T0VAR** The variance of the resin toxicity distribution. Not used in
current

(1/07) source code computations.

**TAMPL** In model simulations, the amplitude of the maximum deviation
of daily

average temp from the yearly mean. Computed as July average

minus **TMEAN** (see below).

**tavg** Long term average temperature on a particular day.

**TC** A matrix, the i,j entry of which is the genetic value of the
resin toxicity

of the clone planted in the clonal block located in the ith row and jth

column of blocks in the stand. Not used in current (12/06) source

code computations.

**tdev** Simulated maximum temperature deviation away from the simulated

average for a particular day.

**tdmean** Mean of **tdev** (see above) - this varies over the year.

**TDVAR** The variance of **tdev** around its mean, **tdmean** (see
above).

**THRESH** Temperature threshold below which no beetle activity occurs,
except

possible mortality, and no degree-days are accumulated.

**tmax** The simulated highest temperature on a particular day.

**TMEAN** The long term yearly mean average temperature (oC). A major
input

parameter.

**tmerch** The total volume (ft^3^/ac) of merchantable timber in the
stand.

**tmin** The simulated lowest temperature on a particular day.

**totegg** The total number of SPB eggs deposited in a particular tree
on a

particular day.

**tp** The number of trees planted per acre.

**tree** A matrix, in which the i,j entry indicates the current status
of the tree

in the i,j position of the stand:

(tree(i,j) = 2 if the tree is alive,

= 1 if the tree died earlier this year from beetle

penetrations,

= 0 if the tree died in a preceding year,

= -1 if the tree is 'fallen' (i.e., died over 5 years

ago and is no longer standing).

**TVAR** The variance of **avtemp** (see above) around its mean,
**tavg** (above).

**urand** Random deviate from a uniform distribution on the interval
\[0,1\].

**vmerch** A matrix, the i,j entry of which is the current number of
cubic feet of

merchantable timber in tree i,j (adapted from **PTAEDA2** -

Burkhart et al. 1987). Only live trees have positive **vmerch**.

**voldbh** Minimum dbh (in) considered in merchantable volume
computation

(see **PTAEDA2** - Burkhart et al. 1987).

**voltop** Minimum upper bole diameter (in) included in merchantable
volume

computation (see **PTAEDA2** - Burkhart et al. 1987).

**voltyp** Indicates the type of merchantable volume reported (see
**PTAEDA2** -

Burkhart et al. 1987). In **spblob**, volume is \'inside bark\'.

**xbeg** Maturity level (cumulative) of a juvenile cohort at the
beginning of a

simulated day.

**xcoord** A matrix, in which the i,j entry indicates the horizontal
distance (in

feet) between tree i,j and the upper left corner of the stand.

**xend** Maturity level (cumulative) of a juvenile cohort at the end of
a

simulated day.

**xxxjuv** Average number of juvenile beetles, per attacked tree.

**xxxland** Average number of landed beetles, per attacked tree.

**xxxpar** Average number of parent beetles, per attacked tree.

**ycoord** A matrix, in which the i,j entry indicates the vertical
distance (in feet)

between tree i,j and the upper left corner of the stand.

**yrdied** A matrix, in which the i,j entry indicates the stand year in
which tree i,j

died. This entry is 0 if the tree lives until harvest. (Trees that die

during the juvenile phase are assigned value 8.) This matrix, together

with matrix **daydied** above, enables us to track the course of beetle-

caused deaths throughout the life of the stand, from planting until

harvest.

**yrnum** A variable whose value is the current stand age in a
simulation.

**VII. REFERENCES**

Baskerville, G.L. and P. Emin (1969). Rapid estimation of heat
accumulation

from maximum and minimum temperatures. *Ecology* **50**(3): 514-517.

Balocchi, C.E., F.E. Bridgwater, B.J. Zobel and S. Jahromi (1993). Age
trends in genetic

parameters for tree height in a nonselected population of loblolly pine.
*Forest*

*Science* **39**: 231-251.

Belanger, R.P., R.L. Hedden and P.L. Lorio, Jr. (1993). Management
strategies to reduce

losses from the southern pine beetle. *Southern Journal of Applied
Forestry*

**17**(3): 150-154.

A.A. Berryman, A.A. (1972). Resistance of conifers to invasion by bark
beetle-fungus

association. *BioScience*, **22**: 598-602.

Bishir, J. and J. Roberds (1995). Analysis of failure time in clonally
propagated plant

populations. *Mathematical Biosciences* **125**: 109-125.

Bishir, J. and J. Roberds (1997). Limit theorems and a general framework
for risk analysis

in clonal forestry. *Mathematical Biosciences* **142**: 1-11.

Bishir, J. and J. Roberds (1999). On numbers of clones needed for
managing risks in

clonal forestry, *Forest Genetics* **6**(3): 149-155.

Bishir, J., J. H. Roberds and B. L. Strom (2004). On-bark behavior of
*Dendroctonus* *frontalis*: A Markov chain analysis. *Journal of Insect
Behavior* **17**:281-301.

# 

Bishir, J., J.H. Roberds, B.L. Strom and X. Wan (2007). SPBLOB: A
simulation model of

the interaction between loblolly pine and southern pine beetle, over the
lifetime of a

plantation. Preprint, available at **www.**srs.fs.usda.gov/4501.

Blanche, C. A., J. D. Hodges, T. E. Nebeker and D. M. Moehring (1983).
Southern

pine beetle: the host dimension. Mississippi Agricultural and Forestry

Experiment Station, Bulletin 917.

Bunt, W.D., J.E. Coster and P.C. Johnson (1980). Behavior of The
Southern

Pine Beetle on the Bark of Host Trees during Mass Attack. *Ann.
Entomol.*

*Soc. Am.* **73**: 647-652.

Burkhart, H.E., K.D. Farrar, R.L. Amateis, and R.F. Daniels (1987).
Simulation

of individual tree growth and development in loblolly pine plantations
on

cutover, site-prepared areas. Publication No. FWS-1-87, School of

Forestry and Wildlife Resources, Virginia Polytechnic Institute and
State

University, Blacksburg, VA.

Coster, J. E. and J. L. Searcy (eds) (1981). Site, stand, and host
characteristics

of southern pine beetle infestations. USDA Combined Forest Pest Research

and Development Program, Technical Bulletin No. 1612.

Coulson, R.N. (1980). Population dynamics. pp. 71-105, In: R. C.
Thatcher,

J. E. Searcy, J. E. Coster and G. D. Hertel (Eds.). *The Southern Pine*

*Beetle*. U.S. Dept.of Agriculture Forest Service, Science and Technical

Administration Technical Bulletin 1631.

Coulson, R.N., A.M. Mayyasi, J.L. Foltz, F.P. Hain and W.C. Martin
(1976).

Resource utilization by the southern pine beetle, *Dendroctonus
frontalis*

(Coleoptera: Scolytidae). *Canadian Entomologist* **108**:353-362.

Coulson, R.N., P.E. Pulley, J.L. Foltz, W.C. Martin and C.L Kelley
(1977).

Survival models for within-tree populations of *Dendroctonus frontalis*

(Coleoptera: Scolytidae). *Can. Ent*. **109**: 1071-1077.

Coulson, R.N., W.S. Fargo, P.E. Pulley, J.L. Foltz, D.N. Pope, J.V.
Richerson

> and T.L. Payne (1978). Evaluation of the re-emergence process of
> parent
>
> adult *Dendroctonus frontalis* (Coleoptera: Scolytidae). *Can. Ent*.
>
> **110**: 475-486.

Coulson, R.N., P.B. Hennier, R.O. Flamm, E.J. Rykiel, L.C. Hu, and T.L
Payne

(1983). The role of lightning in the epidemiology of the Southern Pine

Beetle. *Z. angewande Ent.* **96**: 182-193.

Coulson, R. N., R. M. Feldman, P. J. H. Sharpe, P. E. Pulley, T. L.
Wagner and

T. L. Payne (1989). An overview of the TAMBEETLE model of

*Dendroctonus frontalis* population dynamics. *Holarctic Ecology*
**12**:445-450.

Coyne, J. F. and L. H. Lott (1976). Toxicity of substances in pine
oleoresin to

southern pine beetles. *J. Georgia Entomological Society*
**11**(4):301-305.

Curry, G.L., R.M. Feldman and K.C. Smith (1978). A stochastic model of a
temperature-

dependent population. *Theoretical Population Biology* **13**:197-212.

Daniels, R.F. and H.E. Burkhart (1975). Simulation of individual tree
growth

and stand development in managed loblolly pine plantations. Publication

FWS-5-75, Division of Forestry and Wildlife Resources, Virginia

Polytechnic Institute and State University, Blacksburg, VA.

Fargo, W.S., R.N. Coulson, J.A. Gagne and J.L. Foltz (1979). Correlation
of

southern pine beetle attack density, oviposition, and generation
survival

with host tree characteristics and preceeding beetle life stages within
the

host. *Environmental Entomology* **8**:624-628.

Falconer, D.S. and T.F.C. Mackay (1996). Introduction to Quantitative
Genetics, 4^th^ Ed.,

Longman Group Limited.

Foltz, J.L., A.M. Mayyasi, F.P. Hain, R.N. Coulson, and W.C. Martin
(1976).

Egg-gallery length relationship and within-tree analyses for the
southern

pine beetle, *Dendroctonus frontalis* (Coleoptera: Scolytidae).

*Can. Ent.* **108**:341-352.

Gagne, J.A., T.L. Wagner, P.J.H. Sharpe, R.N. Coulson and W.S. Fargo
(1982).

Reemergence of *Dendroctonus frontalis* (Coleoptera: Scolytidae) at

constant temperatures. *Environ. Entomol*. **11:** 1216-1222.

Gara, R.L., J.P. Vite' and H.H. Cramer (1965). Manipulation of
*Dendroctonus*

*frontalis* by use of a population aggregating pheromone.
*Contributions*

*from Boyce Thompson Institute for Plant Research* **23**:55-66.

Hodges, J.D. and L.S. Pickard (1971). Lightning in the ecology of the

southern pine beetle, *Dendroctonus frontalis* (Coleoptera: Scolytidae).

*Can. Ent*. **103**: 44-51.

Hodges, J.D., W.W. Elam, W.F. Watson and T.E. Nebeker (1979). Oleoresin

characteristics and susceptibility of four southern pines to southern
pine

beetle (Cloeoptera: Scolytidae) attacks. *Canadian Entomologist*

**111**: 889-896.

Libby, W.J. (1982). What is a safe number of clones? In: Resistance to
Disease

and Pests in Forest Trees. (H. Heybroek, B. Stephan and K. von

Weissenberg, eds.). Proceedings of the 3^rd^ International Workshop on
the

Genetics of Host-Parasite Interactions in Forestry. Pp. 342-360.

Wageningen, The Netherlands, Sept. 1980.

Lih, M.P., F.M. Stephen, and K.J. Jeng (1995). *SPBMODEL A User\'s
Guide.*

Offprint, Department of Entomology, University of Arkansas,

Fayetteville, AR, USA.

Lorio, P.L, Jr. (1986). Growth-differentiation balance: a basis for
understanding

southern pine beetle-tree interactions. *Forest Ecology and Management*,

**14**: 259-273.

P.L. Lorio, Jr., R.A. Sommers, C.A. Blanche, J.D. Hodges, and T.E.
Nebeker (1990).

Modeling pine resistance to bark beetles based on growth and
differentiation

balance principles. pp. 402-409, In: R.K. Dixon, R.S. Meldahl, G.A.
Ruark, and

W.G. Warren (Eds.). *Process Modeling of Forest Growth Responses to*

*Environmental Stress*.

Muhs, H.-J. (1993). Policies, regulations and laws affecting clonal
forestry.

In: M. Ahuja and W.J. Libby (Eds.). *Clonal Forestry Vol. 2*, pp
215-227.

Springer-Verlag.

Payne, T. L. (1980). Life history and habits. pp. 7-28, In: R. C.
Thatcher,

J. E. Searcy, J. E. Coster and G. D. Hertel (Eds.). *The Southern Pine*

*Beetle*. U.S. Dept.of Agriculture Forest Service, Science and Technical

Administration Technical Bulletin 1631.

##### Pope, D.N., R.N Coulson, W.S. Fargo, J.A. Gagne and C.W. Kelly (1980). The 

##### allocation process and between-tree survival probabilities in *dendroctonus* 

##### *frontalis* infestations*. Researches in Population Ecology* 22:197-210. 

Press, W.H., S.A. Teukolsky, W.T. Vetterling and B.P. Flannery (1992).
Numerical

Recipes in FORTRAN, 2^nd^ Ed. Cambridge University Press. New York, NY.

Roberds, J.H., G. Namkoong and T. Skroppa (1990). Genetic analysis of
risk in

clonal populations of forest trees. *Theor. Appl. Genet*.
**79**:841-848.

Roberds, J.H. and J. Bishir (1997). Risk analyses in clonal forestry.
*Can. J. For.*

*Res*. **27**:425-432.

Roberds, J.H., B.L. Strom, F.P. Hain, D.P. Gwaze, S.E. McKeand and L.H.
Lott

(2003). Estimates of genetic parameters for oleoresin and growth traits
in

juvenile loblolly pine. *Can. J. For. Res.* **33**: 2469-2476.

Ross, S. (2002). A First Course in Probability, 6^th^ Ed. Prentice-Hall,
Inc. Upper

Saddle River, NJ.

Ruel, J. J., M. P. Ayers, and P. L. Lorio (1998). Loblolly pine responds
to

mechanical wounding with increased resin flow. *Can. J. For. Res.*

**28**: 596-602.

Schultz, R.P. (1997). Loblolly Pine: The Ecology and Culture of Loblolly
Pine

(Pinus taeda L.). Agricultural Handbook 713, U.S. Department of

Agriculture, Forest Service, Washington, D.C.

Sharpe, P.J.H., G.L. Curry, D.W. DeMichele and C.L. Cole (1977).
Distribution

model of organism development times. *Journal of Theoretical Biology*

**66**:21-38.

Stephen, F.M. and M.P. Lih. (1985). A *Dendroctonus frontalis*
infestation

growth model: organization, refinement, and utilization, pp. 186-194,

In: S.B. Branham and R.C. Thatcher (eds.). *Integrated Pest Management*

*Research Symposium: The Proceedings.* U.S. Dept. of Agriculture Forest

Service, Asheville, NC., USA, General Technical Report SO-56.

Thatcher, R.C. (1967). Winter brood development of southern pine beetle
in

southeast Texas. *J Economic Entomology* **60**(2): 599-600.

Thistle, H.W., H. Peterson, G. Allwine, B. Lamb, T. Strand, E.H. Holsten
and

P.J. Shea (2004). Surrogate pheromone plumes in three forest trunk

spaces: Composite statistics and case studies. *Forest Science*
**50**(5):

610-625.

Ungerer, M.J., M.P. Ayres and M.J. Lombardero (1999). Climate and the

northern distribution limits of *Dendroctonus frontalis* Zimmermann

(Coleoptera: Scolytidae). *Journal of Biogeography* **26**:1133-1145.

Wagner, T.L., J.A. Gagne, P.J.H. Sharpe and R.N. Coulson (1983/1984). A

biophysical model of Southern Pine Beetle, *Dendroctonus frontalis*

Zimmermann (Coleoptera: Scolytidae) development. *Ecological*

*Modelling* **21**:125-147.

Wagner, T.L., J.A. Gagne, P.J.H. Sharpe and R.N. Coulson (1984). Effects

of constant temperature on longevity of adult Southern Pine Beetles

(Coleoptera: Scolytidae). *Environmental Entomology* **13**:1125-1130.

**VIII. APPENDIX - LEGALITIES**

**A. Introduction**

The compilers, random number generators and command files used in
connection with the **SPBLOB** code have been downloaded from sites in
the public domain. Some citations and credits appear in **Section V** of
this Documentation. This Appendix contains detailed information
concerning Disclaimers and Legal Requirements associated with these
uses, and also with the use of the **SPBLOB** model and code. Briefly,
we place the model and code in the public domain. However, some of the
code related to random number generation originates in Association for
Computing Machinery (**ACM**) publications. **ACM** retains copyright to
this code, but allows use in any computer for any purpose, with the
restriction that distribution for commercial purposes requires
permission from the **ACM**, and perhaps payment.

**B. Web Sources**

Technical material downloaded for use with the spblob code originates
from three web sites:

-   Subroutines **'get_command_argument'** and
    '**command_argument_count'**

> from the lib file **f2kcli.lib**
>
> [**http://www.winteracter.com/f2kcli/**](http://www.winteracter.com/f2kcli/)

-   **OpenWatcom** Fortran and C code compilers from **SYBASE, INC**

[**http://www.openwatcom.org/**](http://www.openwatcom.org/)

-   Code files for random number generators (**RANDLIB**) from the
    **Software**

> **Download Site of the Department of Biostatistics & Applied
> Mathematics**
>
> **of the University of Texas M.D. Anderson Cancer Center**
>
> [**http://biostatistics.mdanderson.org/SoftwareDownload/SingleSoftware.a
> spx?Software_Id=27**](http://biostatistics.mdanderson.org/SoftwareDownload/SingleSoftware.aspx?Software_Id=27)

**C. Legalities**

**C.1 f2kcli**

While no legal statement is required by **Interactive Software Services
Limited**,

we wish to acknowledge with thanks our use of their **f2kcli** command
line interface.

**C.2 SYBASE, INC.**

**OpenWatcom** Fortran and C code compilers have been used in compliance
with

the **Sybase Open Watcom Public License version 1.0** to compile our
source

code and produce the executable **Interactive** and **Multiple
Simulation** files

included in the **SPBLOB** code folder **spblob_1_25_2007**.

**C.3 RANDLIB**

**C.3.1 From the RANDLIB90 README FILE**

RANDLIB90 is a set of Fortran 90 routines that generate random numbers

from the usual statistical distributions. The complete description of

functionality is found in the documentation files.

======================== FILES =============================

The following files should be present when you unpack this code.

INSTALLATION: Instructions on compiling and installing the code. The

instructions are necessarily vague because almost any computer with

any operating system could contain these programs.

\-\-\-\--

LEGALITIES: Describes the legalities of use of the code. Briefly: We

place our efforts in the public domain, but much of the code

originates from ACM publications. ACM retains copyright to this code

but allows the use in any computer for any purpose. However,

distribution for commercial gain requires permission from the ACM and

perhaps payment.

**C.3.2 The RANDLIB90 LEGALITIES FILE**

Legalities

==========

We place our efforts in writing this package in the public domain.

However, code from ACM publications is subject to the ACM policy

(below).

References

==========

Base Generator

\-\-\-\-\-\-\-\-\-\-\-\-\--

The base generator and all code in the ecuyer_cote_mod come from the

reference:

P. L\'Ecuyer and S. Cote. (1991) \`\`Implementing a Random Number

Package with Splitting Facilities.\'\' ACM Trans. on Math. Softw. 17:1,
pp

98-111.We transliterated the Pascal of the reference to Fortran 95.

The Beta Random Number Generator

================================

R. C. H. Cheng (1978) \`\`Generating Beta Variates with Nonintegral

Shape Parameters.\'\' Communications of the ACM, 21:317-322 (Algorithms

B and BC)

The Binomial Random Number Generator

====================================

Kachitvichyanukul, V. and Schmeiser, B. W. (1988) \`\`Binomial Random

Variate Generation.\'\' Communications of the ACM, 31: 216. (Algorithm

BTPE.)

The Standard Exponential Random Number Generator

================================================

Ahrens, J.H. and Dieter, U. (1972) \`\`Computer Methods for Sampling

from the Exponential and Normal Distributions.\'\' Communications of the

ACM, 15: 873-882. (Algorithm SA.)

The Standard Gamma Random Number Generator

==========================================

Ahrens, J.H. and Dieter, U. (1982) \`\`Generating Gamma Variates by a

Modified Rejection Technique.\'\' Communications of the ACM, 25: 47-54.

(Algorithm SA.)

The Standard Normal Random Number Generator

===========================================

Ahrens, J.H. and Dieter, U. (1973) \`\`Extensions of Forsythe\'s Method

for Random Sampling from the Normal Distribution.\'\' Math. Comput.,

27:927 - 937. Algorithm FL (method=5)

ACM Policy on Use of Code

=========================

Here is the software Policy of the ACM.

Submittal of an algorithm for publication in one of the

ACM Transactions implies that unrestricted use of the

algorithm within a computer is permissible. General

permission to copy and distribute the algorithm without

fee is granted provided that the copies are not made or

distributed for direct commercial advantage. The ACM

copyright notice and the title of the publication and

its date appear, and notice is given that copying is by

permission of the Association for Computing

Machinery. To copy otherwise, or to republish, requires

a fee and/or specific permission. Krogh, F. (1997)

\`\`Algorithms Policy.\'\' ACM Tran. Math. Softw. 13,

183-186.

We do not know the policy of the Royal Statistical Society; they have

discontinued publishing algorithms. However, they made a number of

these programs available on Statlib on condition that there be no

charge for their distribution.

Here is our standard disclaimer.

NO WARRANTY

> WE PROVIDE ABSOLUTELY NO WARRANTY OF ANY KIND EITHER EXPRESSED OR
> IMPLIED, INCLUDING BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
> MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK
> AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD
> THIS PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY
> SERVICING, REPAIR OR CORRECTION. IN NO EVENT SHALL THE UNIVERSITY OF
> TEXAS OR ANY OF ITS COMPONENT INSTITUTIONS INCLUDING M. D. ANDERSON
> HOSPITAL BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY LOST PROFITS,
> LOST MONIES, OR OTHER SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES
> ARISING OUT OF THE USE OR INABILITY TO USE (INCLUDING BUT NOT LIMITED
> TO LOSS OF DATA OR DATA OR ITS ANALYSIS BEING RENDERED INACCURATE OR
> LOSSES SUSTAINED BY THIRD PARTIES) THE PROGRAM.

(Above NO WARRANTY modified from the GNU NO WARRANTY statement.)