{\rtf1\ansi\ansicpg1252\cocoartf949\cocoasubrtf540
{\fonttbl\f0\fswiss\fcharset0 Helvetica;}
{\colortbl;\red255\green255\blue255;}
\margl1440\margr1440\vieww17100\viewh10780\viewkind0
\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\ql\qnatural\pardirnatural

\f0\fs24 \cf0 ## to get data for a specific year ##\
data2.2002 = read.table("w6 2002.txt", sep = ",", header = T, na.strings="-")\
\
##define e values##\
e.values = read.csv("e values.csv", header=FALSE)\
\
## make the species codes and vigor codes text ##\
data2.2002[,3]=as.character(data2.2002[,3])\
\
## read in parameters table and separate the character column from the numeric columns ##\
parameters.2002 = read.csv("parameters2.csv")\
parameters1.2002 = as.character(parameters.2002[,1])\
parameters2.2002 = as.matrix(parameters.2002[,2:length(parameters.2002[1,])])\
\
## matrix for all parameters ##\
parameters.all.2002 = matrix(0,length(data2.2002[,1]), length(parameters2.2002[1,]))\
\
## fill the parameters matrix ##\
for(j in 1:length(data2.2002[,1]))\{\
parameters.all.2002[j,] = subset(parameters2.2002, parameters1.2002==data2.2002[j,3])\
\}\
\
## read in vigor parameters table and separate the character column from the numeric columns ##\
VC.2002 = read.csv("vigor parameters.csv")\
VC1.2002 = as.character(VC.2002[,1])\
VC2.2002 = as.matrix(VC.2002[,2:length(VC.2002[1,])])\
\
## matrix for all vigor codes ##\
VC.all.2002 = matrix(0,length(data2.2002[,1]), length(VC2.2002[1,]))\
\
## fill the parameters matrix ##\
for(j in 1:length(data2.2002[,1]))\{\
VC.all.2002[j,] = subset(VC2.2002, VC1.2002==data2.2002[j,8])\
\}\
\
## define DBH matrix ##\
DBH.2002 = matrix(0,length(data2.2002[,1]),1000)\
\
## make 1000 estimates of DBH [cm] for each tree ##\
for(i in 1:length(data2.2002[,1]))\{\
DBH.2002[i,] = (rnorm(1000, 0, 1)) * 0.05 + data2.2002[i,7]\
\}\
\
## define SB##\
SB.2002  = sqrt(parameters.all.2002[,22]+((log10(data2.2002[,7])-parameters.all.2002[,20])^2/parameters.all.2002[,21]))\
\
\
## define correction factors##\
correction.factor.stem.wood.2002 = parameters.all.2002[,23]\
correction.factor.stem.bark.2002 = parameters.all.2002[,24]\
correction.factor.branches.2002 = parameters.all.2002[,25]\
correction.factor.leaves.2002 = parameters.all.2002[,26]\
correction.factor.twigs.2002 = parameters.all.2002[,26]\
correction.factor.roots.2002 = parameters.all.2002[,27]\
correction.factor.dead.branches.2002 = parameters.all.2002[,55]\
\
##define leaf ratios##\
leaf.ratio.leaves.2002 = parameters.all.2002[,19]\
leaf.ratio.twigs.2002 = (1-(parameters.all.2002[,19]))\
\
\
##define area matrix##\
plot.area.2002 = matrix(0,length(data2.2002[,1]),1)\
\
\
##fill area matrix. areas are in m^2##\
for(i in 1:length(data2.2002[,1]))\{\
if(data2.2002[i,7]>=10)\{plot.area.2002[i]=130000\} else \
if(data2.2002[i,7]<10)\{plot.area.2002[i]=15702.48\}\
\}\
\
\
## matrices for biomass variables ##\
height.2002 = matrix(0,length(data2.2002[,1]),1000)\
stem.wood.2002 = matrix(0,length(data2.2002[,1]),1000)\
stem.bark.2002 = matrix(0,length(data2.2002[,1]),1000)\
branches.2002= matrix(0,length(data2.2002[,1]),1000)\
leaves.2002 = matrix(0,length(data2.2002[,1]),1000)\
twigs.2002 = matrix(0,length(data2.2002[,1]),1000)\
roots.2002 = matrix(0,length(data2.2002[,1]),1000)\
dead.branches.2002 =  matrix(0,length(data2.2002[,1]),1000)\
\
## calculate biomass compartments 1000 times with the DBH, a, b, and E##\
for(i in 1:1000)\{\
height.2002[,i] =(10^(parameters.all.2002[,1] + (parameters.all.2002[,2] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,3])*e.values[i,1]*SB.2002)) \
stem.wood.2002[,i] = (10^(parameters.all.2002[,4] + (parameters.all.2002[,5] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,6])*e.values[i,2]*SB.2002))*correction.factor.stem.wood.2002\
stem.bark.2002[,i] = (10^(parameters.all.2002[,7] + (parameters.all.2002[,8] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,9])*e.values[i,3]*SB.2002))*correction.factor.stem.bark.2002\
branches.2002[,i] = (10^(parameters.all.2002[,10] + (parameters.all.2002[,11] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,12])*e.values[i,4]*SB.2002))*correction.factor.branches.2002\
leaves.2002[,i] = (10^(parameters.all.2002[,13] + (parameters.all.2002[,14] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,15])*e.values[i,5]*SB.2002))*correction.factor.leaves.2002*leaf.ratio.leaves.2002\
twigs.2002[,i] = (10^(parameters.all.2002[,13] + (parameters.all.2002[,14] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,15])*e.values[i,6]*SB.2002))*correction.factor.twigs.2002*leaf.ratio.twigs.2002\
roots.2002[,i] = (10^(parameters.all.2002[,16] + (parameters.all.2002[,17] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,18])*e.values[i,7]*SB.2002))*correction.factor.roots.2002\
dead.branches.2002[,i] = (10^(parameters.all.2002[,52] + (parameters.all.2002[,53] * log10(DBH.2002[,i])) + log10(parameters.all.2002[,54])*e.values[i,20]*SB.2002))*correction.factor.dead.branches.2002\
\}\
\
## matrices for vigor codes and area (g/m2) ##\
stem.wood.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
stem.bark.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
branches.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
leaves.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
twigs.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
roots.2002.gm2 = matrix(0,length(data2.2002[,1]),1000)\
dead.branches.2002.gm2 =  matrix(0,length(data2.2002[,1]),1000)\
\
##calculate biomass compartments 1000 times with vigor code##\
for(i in 1:1000)\{\
stem.wood.2002.gm2[,i] = stem.wood.2002[,i]*VC.all.2002[,8]*VC.all.2002[,1]/plot.area.2002\
stem.bark.2002.gm2[,i] = stem.bark.2002[,i]*VC.all.2002[,8]*VC.all.2002[,2]/plot.area.2002\
branches.2002.gm2[,i] = branches.2002[,i]*VC.all.2002[,8]*VC.all.2002[,3]/plot.area.2002\
leaves.2002.gm2[,i] = leaves.2002[,i]*VC.all.2002[,8]*VC.all.2002[,4]/plot.area.2002\
twigs.2002.gm2[,i] = twigs.2002[,i]*VC.all.2002[,8]*VC.all.2002[,5]/plot.area.2002\
roots.2002.gm2[,i] = roots.2002[,i]*VC.all.2002[,8]*VC.all.2002[,6]/plot.area.2002\
dead.branches.2002.gm2[,i] = dead.branches.2002[,i]*VC.all.2002[,8]*VC.all.2002[,7]/plot.area.2002\
\}\
\
\
##sum columns for total biomass (gm2) for all trees, 1000 Monte Carlo estimates##\
biomass.sums.2002=colSums(stem.wood.2002.gm2 + stem.bark.2002.gm2 + branches.2002.gm2 + leaves.2002.gm2 + twigs.2002.gm2+ roots.2002.gm2 + dead.branches.2002.gm2)\
\
##convert g/m2 to kg/Ha##\
biomass.sums.2002.MgHa=biomass.sums.2002/100\
\
##find 95%CI ( these will be the 26th highest and 26th lowest values)##\
quantile(biomass.sums.2002.MgHa, 0.975)\
quantile(biomass.sums.2002.MgHa, 0.025)\
\
\
## matrices for nutrient mass (N) ##\
stem.wood.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
stem.bark.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
branches.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
leaves.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
twigs.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
roots.nutrient.2002.N = matrix(0,length(data2.2002[,1]),1000)\
dead.branches.nutrient.2002.N =  matrix(0,length(data2.2002[,1]),1000)\
\
## calculate nutrient compartments 1000 times with the biomass, nutrient content, and nutrient content E (in grams)##\
for(i in 1:1000)\{\
stem.wood.nutrient.2002.N[,i] = stem.wood.2002.gm2[,i] *((parameters.all.2002[,28] + (parameters.all.2002[,29] *e.values[i,8]))/100)\
stem.bark.nutrient.2002.N[,i] = stem.bark.2002.gm2[,i] *((parameters.all.2002[,30] + (parameters.all.2002[,31] *e.values[i,9]))/100)\
branches.nutrient.2002.N[,i] = branches.2002.gm2[,i] *((parameters.all.2002[,32] + (parameters.all.2002[,33] *e.values[i,10]))/100)\
leaves.nutrient.2002.N[,i] = leaves.2002.gm2[,i] *((parameters.all.2002[,34] + (parameters.all.2002[,35] *e.values[i,11]))/100)\
twigs.nutrient.2002.N[,i] = twigs.2002.gm2[,i] *((parameters.all.2002[,36] + (parameters.all.2002[,37] *e.values[i,12]))/100)\
roots.nutrient.2002.N[,i] = roots.2002.gm2[,i] *((parameters.all.2002[,38] + (parameters.all.2002[,39] *e.values[i,13]))/100)\
dead.branches.nutrient.2002.N[,i] = dead.branches.2002.gm2[,i] *((parameters.all.2002[,32] + (parameters.all.2002[,33] *e.values[i,21]))/100)\
\}\
\
##sum columns for total nutrients (gm2) for all trees, 1000 Monte Carlo estimates##\
nutrient.sums.2002.N=colSums(stem.wood.nutrient.2002.N + stem.bark.nutrient.2002.N + branches.nutrient.2002.N + leaves.nutrient.2002.N + twigs.nutrient.2002.N + roots.nutrient.2002.N + dead.branches.nutrient.2002.N)\
\
##convert g/m2 to kg/Ha##\
nutrient.sums.2002.kgHa.N=nutrient.sums.2002.N*10\
\
\
##sort from highest to lowest to find 95%CI ( these will be the 26th highest and 26th lowest values)##\
quantile(nutrient.sums.2002.kgHa.N, 0.975)\
quantile(nutrient.sums.2002.kgHa.N, 0.025)\
\
## matrices for nutrient mass (P) ##\
stem.wood.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
stem.bark.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
branches.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
leaves.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
twigs.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
roots.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
dead.branches.nutrient.2002.P = matrix(0,length(data2.2002[,1]),1000)\
\
## calculate nutrient compartments 1000 times with the biomass, nutrient content, and nutrient content E (in grams)##\
for(i in 1:1000)\{\
stem.wood.nutrient.2002.P[,i] = stem.wood.2002.gm2[,i] *((parameters.all.2002[,40] + (parameters.all.2002[,41] *e.values[i,14]))/100)\
stem.bark.nutrient.2002.P[,i] = stem.bark.2002.gm2[,i] *((parameters.all.2002[,42] + (parameters.all.2002[,43] *e.values[i,15]))/100)\
branches.nutrient.2002.P[,i] = branches.2002.gm2[,i] *((parameters.all.2002[,44] + (parameters.all.2002[,45] *e.values[i,16]))/100)\
leaves.nutrient.2002.P[,i] = leaves.2002.gm2[,i] *((parameters.all.2002[,46] + (parameters.all.2002[,47] *e.values[i,17]))/100)\
twigs.nutrient.2002.P[,i] = twigs.2002.gm2[,i] *((parameters.all.2002[,48] + (parameters.all.2002[,49] *e.values[i,18]))/100)\
roots.nutrient.2002.P[,i] = roots.2002.gm2[,i] *((parameters.all.2002[,50] + (parameters.all.2002[,51] *e.values[i,19]))/100)\
dead.branches.nutrient.2002.P[,i] = dead.branches.2002.gm2[,i] *((parameters.all.2002[,44] + (parameters.all.2002[,45] *e.values[i,22]))/100)\
\}\
\
##sum columns for total nutrients (gm2) for all trees, 1000 Monte Carlo estimates##\
nutrient.sums.2002.P=colSums(stem.wood.nutrient.2002.P + stem.bark.nutrient.2002.P + branches.nutrient.2002.P + leaves.nutrient.2002.P + twigs.nutrient.2002.P + roots.nutrient.2002.P + dead.branches.nutrient.2002.P)\
\
##convert g/m2 to kg/Ha##\
nutrient.sums.2002.kgHa.P=nutrient.sums.2002.P*10\
\
\
##sort from highest to lowest to find 95%CI ( these will be the 26th highest and 26th lowest values)##\
quantile(nutrient.sums.2002.kgHa.P, 0.975)\
quantile(nutrient.sums.2002.kgHa.P, 0.025)\
\
par(mfrow=c(1,3))\
hist(biomass.sums.2002.MgHa, main = "Biomass iterations (2002)")\
hist(nutrient.sums.2002.kgHa.N, main = "N iterations (2002)")\
hist(nutrient.sums.2002.kgHa.P, main = "P iterations (2002)")\
par(mfrow=c(3,3))\
}