## 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))