Attribute VB_Name = "Module1"
'----------------------------------------------------------------
'Routines in this modules are linked with entry point
'at the RunModel routine. System projects stand tables for
'forest strata and calculates harvestable volumes
'This module is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Option Explicit
Const ProgID = "MYRLIN#3 1.05b"

'****** global variable declarations ********
Dim Yws As Worksheet            '{yield} worksheet used for outputs
Dim St() As Double              'stand table (classes x groups x strata)
Dim DCwidth() As Double         'diameter class widths
Dim DCmidBA() As Double         'diameter class mid-point BA
Dim DCdiam() As Double          'diameter class mid-diameter
Dim Nst As Integer              'number of forest strata
Dim Ndc As Integer              'number of diameter classes
Dim Ngp As Integer              'number of species groups
Dim Tp As Integer               'time period
Dim DmgF As Double              'logging damage factor
Dim FellCyc As Integer          'felling cycle
Dim Tlim As Integer             'simulation time limit
Dim AreaC As Double             'area of periodic coupe
Dim LuL As Integer              'index of last unit logged
Dim aTol As Double              'area tolerance
Dim RecF As Double              'recruitment factor

Type StratumDef
   Name As String               'stratum name
   SeqNo As Integer             'sequence number
   Area As Double               'area
   Nloss As Double              'tree losses from mortality (N/km2)
   Vst As Double                'current total volume
   Vcom As Double               'current commercial volume (before harvest)
   Vh As Double                 'current harvested volume
End Type

Dim Strata() As StratumDef           'stratum codes
Dim StList As Variant

Type GroupDef
   Name As String
   D95 As Double
   Dinc As Double
   Amr As Double
   Dlim As Double
   Hpct As Double
   FormHt As Double
End Type
Dim Grps() As GroupDef             'groups codes
Dim GpList As Variant
Dim opr As Integer                 'output row

Public Sub RunModel()
'----------------------------------------------------------------
'Sequences the main code blocks in the simulation model
'This module is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim t As Integer
Application.StatusBar = ProgID
ReadGroups
ReadStandTable
PrepOutputs
For t = 0 To Tlim Step Tp
   GrowStand t
   DoHarvest t
   WriteSummary t
Next t
End Sub

Private Sub PrepOutputs()
'prepares output sheet
Dim nr As Integer
Set Yws = Sheets("Yield")
Yws.Activate
nr = Yws.UsedRange.Rows.Count
If nr < 2 Then nr = 2
Yws.Range(Yws.Cells(2, 1), Yws.Cells(nr, 6)).ClearContents
opr = 1
With Sheets("Table1")
    nr = .UsedRange.Rows.Count
    If nr < 2 Then nr = 2
    .Range(.Cells(2, 1), .Cells(nr, 6)).ClearContents
End With
End Sub

Private Sub ReadGroups()
'----------------------------------------------------------------
'Reads species groups and model data, and also stratum and
'area data into internal data structures
'This module is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim wks As Worksheet
Dim k As Integer, j As Integer
Dim hr As Integer
'initialise sensitive globals
LuL = 0
GpList = Empty
'read species group information
Set wks = Sheets("Models")
'general factors
FellCyc = wks.Cells(2, 2)
DmgF = wks.Cells(3, 2)
Tp = wks.Cells(4, 2)
Tlim = wks.Cells(5, 2)
aTol = wks.Cells(6, 2)
RecF = wks.Cells(7, 2)
hr = 10
Do While wks.Cells(k + hr + 1, 1).text > ""
   k = k + 1
   ReDim Preserve Grps(1 To k)
   Grps(k).Name = wks.Cells(k + hr, 1)
   Grps(k).D95 = wks.Cells(k + hr, 2)
   Grps(k).Dinc = wks.Cells(k + hr, 3)
   Grps(k).Amr = wks.Cells(k + hr, 4)
   Grps(k).Dlim = wks.Cells(k + hr, 5)
   Grps(k).Hpct = wks.Cells(k + hr, 6)
   Grps(k).FormHt = wks.Cells(k + hr, 7)
   j = FindA(Grps(k).Name, GpList)
   If j <> k Then Stop  'algorithm check
Loop
Ngp = k
'read forest stratum data
Set wks = Sheets("Areas")
hr = 2
k = 0
AreaC = 0
Do While wks.Cells(k + hr + 1, 1).text > ""
   k = k + 1
   ReDim Preserve Strata(1 To k)
   With Strata(k)
      .Name = wks.Cells(k + hr, 1)
      .SeqNo = wks.Cells(k + hr, 2)
      .Area = wks.Cells(k + hr, 3)
      .Nloss = 0
   End With
   AreaC = AreaC + Strata(k).Area
   FindA Strata(k).Name, StList
Loop
AreaC = AreaC / FellCyc * Tp           'coupe area
Nst = k
End Sub

Private Sub ReadStandTable()
'----------------------------------------------------------------
'Reads in a general stand table into memory
'This module is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim Dclass As Variant       'list of diameter class lower bounds
Dim Tbl As Worksheet        'stand table sheet
Dim d As Integer            'diameter class index
Dim g As Integer            'species group index
Dim u As Integer            'stratum index
Dim r As Integer            'table row index
Dim k As Integer            'secondary diameter class index
Dim Lk As Integer           'last k value with valid class
Dim r0 As Integer           'initial output row for a table
Set Tbl = Sheets("StandTables")
Tbl.Activate
'set diameter class limits.  Halt if classes not well
'organized (ascending sequence, no overlaps)
Dclass = GetClasses([c2])
Ndc = 0
For k = 1 To UBound(Dclass)
   If IsArray(Dclass(k)) Then
      Ndc = Ndc + 1
      Lk = k
      ReDim Preserve DCwidth(1 To Ndc)
      ReDim Preserve DCmidBA(1 To Ndc)
      ReDim Preserve DCdiam(1 To Ndc)
      DCwidth(Ndc) = Dclass(k)(1) - Dclass(k)(0) + 1
      If DCwidth(Ndc) <= 0 Then ErrorMsg "Bad diameter class: " + Cells(2, 3 + k).text
      DCdiam(Ndc) = (Dclass(k)(0) + DCwidth(Ndc) / 2)
      DCmidBA(Ndc) = DCdiam(Ndc) ^ 2 * 0.00007854
      If Ndc > 1 Then
         If DCmidBA(Ndc) <= DCmidBA(Ndc - 1) Then
            ErrorMsg "Diameter class out of sequence: " + Cells(2, 3 + k).text
         End If
      End If
   End If
Next k
'adjust the BA mid point of last class if too wide
If DCwidth(Ndc) > 10 Then
   DCmidBA(Ndc) = (Dclass(Lk)(0) + 5) ^ 2 * 0.00007854
End If
'set up stand table arrays
ReDim St(1 To Ndc, 1 To Ngp, 1 To Nst)
'loop through forest strata
r0 = 3: r = 3
Do While Cells(r0, 1).text > ""
   u = FindB(Cells(r0, 1).text, StList)
   If u <= 0 Then ErrorMsg "Stratum " + Cells(r0, 1).text + " is not defined in stratum list" + _
         "  (error noted at cell " + Cells(r0, 1).Address + ", sheet " + ActiveSheet.Name + ")"
   'loop through groups
   Do
      g = FindB(Cells(r, 2).text, GpList)
      If g <= 0 Then ErrorMsg "Group " + Cells(r, 2).text + " is not defined in group list" + _
         "  (error noted at cell " + Cells(r, 2).Address + ", sheet " + ActiveSheet.Name + ")"
      For d = 1 To Ndc
         On Error Resume Next
         St(d, g, u) = Cells(r, d + 2)
         If Err.Number > 0 Then
            Err.Clear
            ErrorMsg "Invalid stand table data. N/km2 values expected." + _
               "  (error noted at cell " + Cells(r, d + 2).Address + ", sheet " + ActiveSheet.Name + ")"
         End If
      Next d
      r = r + 1
   Loop Until Cells(r, 2).text = "" Or Cells(r, 1).text > ""
   r0 = r
Loop
End Sub

Private Sub WriteSummary(t As Integer)
'output stratum data
Dim u As Integer, r As Integer
Dim Vst As Double, Vc As Double, Vh As Double, aTot As Double
For u = 1 To Nst
   Yws.Cells(opr + u, 1) = t
   Yws.Cells(opr + u, 2) = Strata(u).Name
   Yws.Cells(opr + u, 3) = Strata(u).Area
   Yws.Cells(opr + u, 4) = Strata(u).Vst * Strata(u).Area
   Yws.Cells(opr + u, 5) = Strata(u).Vcom * Strata(u).Area
   If Strata(u).Vh > 0 Then
      Yws.Cells(opr + u, 6) = Strata(u).Vh * Strata(u).Area
   End If
   aTot = aTot + Strata(u).Area
   Vst = Vst + Strata(u).Vst * Strata(u).Area
   Vc = Vc + Strata(u).Vcom * Strata(u).Area
   Vh = Vh + Strata(u).Vh * Strata(u).Area
   Strata(u).Vh = 0
Next u
With Sheets("Table1")
    r = t / Tp + 2
    .Cells(r, 1) = t
    .Cells(r, 2) = Vst / aTot
    .Cells(r, 3) = Vc
    .Cells(r, 4) = Vh
    .Cells(r, 5) = Vh / aTot / Tp
End With
opr = opr + Nst
End Sub

Private Sub GrowStand(t As Integer)
'calculate growth, mortality and recruitment on stand
Dim d As Integer            'diameter class index
Dim g As Integer            'species group index
Dim u As Integer            'stratum index
Dim r As Integer            'relative output row
Dim fm As Double            'moving fraction in a class
Dim tm As Double            'moving trees in a class
Dim fd As Double            'dead trees from a class
Dim Vst As Double           'volume in a stratum
Dim Vcom As Double          'commercial volume
Dim V As Double             'volume of a table cell
Dim Nlost As Double         'number of trees lost in a stratum
Dim GrpWt() As Double       'group weights for recruitment
Dim nTot As Double          'debug output - total N/km, unit 1
Dim vTot() As Double        'do. - total m3/ha, unit 1 by spp.
For u = 1 To Nst
   Nlost = Strata(u).Nloss  'harvest and damage losses last cycle
   Strata(u).Nloss = 0      'reset
   ReDim GrpWt(0 To Ngp)
   ReDim vTot(0 To Ngp)
   Vst = 0: Vcom = 0
   nTot = 0
   For g = 1 To Ngp
      'register species weights for recruitment
      For d = 1 To Ndc
         GrpWt(g) = GrpWt(g) + St(d, g, u)
         GrpWt(0) = GrpWt(0) + St(d, g, u)
      Next d
      'adjust stand table for mortality
      fd = 1 - (1 - Grps(g).Amr) ^ Tp
      For d = 1 To Ndc
         Nlost = Nlost + St(d, g, u) * fd
         St(d, g, u) = St(d, g, u) * (1 - fd)
      Next d
      'adjust stand table for moving trees
      For d = Ndc - 1 To 1 Step -1
         fm = Grps(g).Dinc * Tp / DCwidth(d)
         If fm > 1 Then
            ErrorMsg "MYRLIN is a simple stand projection tool that " + _
               "cannot work if trees move more than one diameter " + _
               "class in a growth period.  Either increase class " + _
               "width or reduce time step.  Check increment is " + _
               "correct for species '" + Grps(g).Name + "'"
         End If
         tm = fm * St(d, g, u)               'moving trees
         St(d + 1, g, u) = St(d + 1, g, u) + tm
         St(d, g, u) = St(d, g, u) - tm
      Next d
   Next g
   If GrpWt(0) <= 0 Then ErrorMsg "Stratum " + Strata(u).Name + " is in the {Areas} list," + _
      " but is not in the data.  Please remove from the list before running the model."
   'add recruitment to first class and add volumes
   For g = 1 To Ngp
      St(1, g, u) = St(1, g, u) + Nlost * GrpWt(g) / GrpWt(0) * RecF
      For d = 1 To Ndc
         'volumes for output
         V = St(d, g, u) * DCmidBA(d) * Grps(g).FormHt * 0.01
         Vst = Vst + V
         If (DCdiam(d) - DCwidth(d) / 2) > Grps(g).Dlim And Grps(g).Hpct > 0 Then
            Vcom = Vcom + V
         End If
         'totals for debugging
         nTot = nTot + St(d, g, u)
         vTot(g) = vTot(g) + V
         vTot(0) = vTot(0) + V
      Next d
   Next g
   Strata(u).Vst = Vst
   Strata(u).Vcom = Vcom
   'output debugging info only for unit 1
   #If 0 Then
   If u = 1 Then
      With Sheets("debug")
         r = t / Tp + 2
         .Cells(r, 1) = t
         .Cells(r, 2) = nTot
         For g = 0 To Ngp
            .Cells(r, g + 3) = vTot(g)
         Next g
      End With
   End If
   #End If
Next u
End Sub

Private Sub DoHarvest(t As Integer)
'calculate growth, mortality and recruitment on stand
Dim d As Integer            'diameter class index
Dim g As Integer            'species group index
Dim u As Integer            'stratum index
Dim r As Integer            'relative output row
Dim fm As Double            'moving trees in a class
Dim fd As Double            'dead trees from a class
Dim Vh As Double            'volume harvested
Dim V As Double             'volume of a table cell
Dim AreaH As Double         'area harvested in this operation
Dim Nh As Double            'number harvested in a table cell
Dim df As Double            'damage factor
'check the area of the next unit to log
If AreaC <= 0 Then Exit Sub   'no harvesting required
'continue harvesting until area cut within tolerance of coupe area
AreaH = 0
Do While Abs(AreaH - AreaC) / AreaC > aTol And AreaH < AreaC
   'get the next unit in sequence
   If LuL >= Nst Then LuL = 0
   u = LuL + 1
   'check if the area to be logged would exceed tolerance
   If Strata(u).Area + AreaH > AreaC * (1 + aTol) Then
      'sub-divide next coupe to leave unlogged portion
      SplitCoupe t, u, AreaC - AreaH
   End If
   'do harvesting on this unit
   Vh = 0
   Strata(u).Nloss = 0
   For g = 1 To Ngp
      If Grps(g).Hpct > 0 Then
         For d = Ndc To 1 Step -1
            If (DCdiam(d) - DCwidth(d) / 2) > Grps(g).Dlim Then
               V = DCmidBA(d) * Grps(g).FormHt * 0.01
               Nh = St(d, g, u) * Grps(g).Hpct
               Strata(u).Nloss = Strata(u).Nloss + Nh
               Vh = Vh + V * Nh
               St(d, g, u) = St(d, g, u) - Nh
            End If
         Next d
      End If
   Next g
   'adjust residual stand for logging damage
   If Vh > 0 And Strata(u).Vst > 0 Then
      df = Vh / Strata(u).Vst * DmgF
      For g = 1 To Ngp
         For d = 1 To Ndc
            Nh = St(d, g, u) * df
            Strata(u).Nloss = Strata(u).Nloss + Nh
            St(d, g, u) = St(d, g, u) - Nh
         Next d
      Next g
   End If
   AreaH = AreaH + Strata(u).Area
   Strata(u).Vh = Vh
   LuL = LuL + 1
Loop
End Sub

Private Sub SplitCoupe(t As Integer, uS As Integer, Area As Double)
'makes two copies of harvesting coupe
Dim d As Integer            'diameter class index
Dim g As Integer            'species group index
Dim u As Integer            'stratum index
'insert a new coupe ID into the list
Nst = Nst + 1
ReDim Preserve Strata(1 To Nst), St(1 To Ndc, 1 To Ngp, 1 To Nst)
For u = Nst To uS + 1 Step -1
   Strata(u) = Strata(u - 1)
   For g = 1 To Ngp
      For d = 1 To Ndc
         St(d, g, u) = St(d, g, u - 1)
      Next d
   Next g
Next u
Strata(uS).Area = Area
Strata(uS).Name = Strata(uS).Name + "{" + Trim(CStr(t)) + "}"
If Area >= Strata(uS + 1).Area Then Stop
Strata(uS + 1).Area = Strata(uS + 1).Area - Area
End Sub

Private Sub ErrorMsg(text As String)
'show error message and terminate program
MsgBox text, vbExclamation, ProgID
End
End Sub

Private Function FindA(e As String, a As Variant) As Integer
'----------------------------------------------------------------
'Searches for index of item e in array a and returns its
'index.  Adds it to the list if not already in.
'This routine is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim i As Integer
If Not IsArray(a) Then
    'list not yet started - add e as first item
    a = Array(e)
    FindA = 1
Else
    For i = 0 To UBound(a)
        If a(i) = e Then
            FindA = i + 1
            Exit Function
        End If
    Next i
    'not yet in list - add it
    i = UBound(a) + 1
    ReDim Preserve a(i)
    a(i) = e
    FindA = i + 1
End If
End Function

Private Function FindB(e As String, a As Variant) As Integer
'----------------------------------------------------------------
'Searches for index of item e in array a and returns its
'index.  If not found, returns 0.
'This routine is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim i As Integer
If Not IsArray(a) Then
    'no list defined
    FindB = 0
Else
    For i = 0 To UBound(a)
        If a(i) = e Then
            FindB = i + 1
            Exit Function
        End If
    Next i
    'not in list - return zero
    FindB = 0
End If
End Function


Private Function GetClasses(cell As Range) As Variant
'----------------------------------------------------------------
'Sets diameter classes in a range of cells according to protocol:
'  a-b  converted to array(a,b)
'  a+   converted to scalar a
'Reads across from <cell> of active sheet until blank encountered.
'Halts with error message if invalid diameter class encountered.
'
'This routine is part of the MYRLIN toolkit.
'For information and conditions of use see http://www.myrlin.org.
'----------------------------------------------------------------
Dim ClassSpec As Range      'range of cells to right of <cell>
Dim ThisCell As Variant     'current cell being examined
Dim Classes As Variant      'array of class specs
Dim a As Integer, b As Integer  'class bounds decoded in cell
Dim k As Integer            'position of - or + character in cell text
Dim j As Integer            'class index
Dim nc As Integer           'class count
Set ClassSpec = Range(cell, cell.End(xlToRight))
nc = ClassSpec.Columns.Count
Classes = Array()
ReDim Classes(1 To nc)
For Each ThisCell In ClassSpec
    k = InStr(ThisCell.text, "-")
    If k > 0 Then
        'normal diameter class range
        a = 20000: b = 0
        On Error Resume Next
        a = CInt(Left(ThisCell.text, k - 1))
        b = CInt(Mid(ThisCell.text, k + 1)) + 1
        On Error GoTo 0
        If b <= a Then GoTo BadClassSpec
        j = j + 1
        Classes(j) = Array(a, b)
    Else
        k = InStr(ThisCell.text, "+")
        If k > 0 Then
            'cumulative diameter class lower bound
            a = 0
            On Error Resume Next
            a = CInt(Left(ThisCell.text, k - 1))
            On Error GoTo 0
            If a <= 0 Then GoTo BadClassSpec
            j = j + 1
            Classes(j) = a
        Else
            GoTo BadClassSpec
        End If
    End If
Next ThisCell
'return list of class bounds
GetClasses = Classes
Exit Function
BadClassSpec:
    MsgBox "Bad diameter class specification at " + _
    ThisCell.Address + ": '" + ThisCell.text + "'", _
    vbExclamation + vbOKOnly, "MYRLIN stand table tool"
    ThisCell.Select
    End
End Function