V6 roads

Download Model
Download the .scn, .sel and .lse files by clicking on the following link:

Model Code Exploration
In the following sections we will examine all of the model files for this model. Note that instead of downloading the zip file above, you could just copy the text in the boxes below into a text editor and save it with the appropriate name (Section title). Opening the resulting .scn file in the SELES simulator would run this model.

BaseScenario.scn
SELES Scenario $gisData$ = ..\gisData\cell // Static layers Ecoregion = $gisData$\ecoregion Elevation = $gisData$\elevation Aspect = $gisData$\aspect ForestStatus = $gisData$\ForestStatus StreamClass = $gisData$\water SiteType = $gisData$\SiteEcoregion // dynamic forest layers initialAge = $gisData$\StandAgeEst initialSpp = $gisData$\DomSpecies // dynamic road layers initialDistanceToRoads = $gisData$\Dist2ActiveRoad initialNearestRoadLoc = $gisData$\NearestRoadSegmentLoc Model Dimensions: Ecoregion Model.sel SimPriority Low Priority // Set low priority for simulation engine // Set up display Minimize Static Minimize Initial State Tile // Move to output folder cwd output // create output directories mkdir cell mkdir cellhd // Base Time step (in years) BaseTimestep = 10 // Make it spatial: 10-100ha blocks BlockSizeMin = 10 BlockSizeMax = 100 // No explicit adjacency GreenupYears = 0 // Use area based harvesting (default) AreaBasedHarvest = TRUE // Make it aspatial: No road access restrictions MinNoEffectDist2Rd = 200 MaxDist2Road = 2000

Scenario1.scn
SELES Scenario $gisData$ = ..\gisData\cell // Static layers Ecoregion = $gisData$\ecoregion Elevation = $gisData$\elevation Aspect = $gisData$\aspect ForestStatus = $gisData$\ForestStatus StreamClass = $gisData$\water SiteType = $gisData$\SiteEcoregion // dynamic forest layers initialAge = $gisData$\StandAgeEst initialSpp = $gisData$\DomSpecies // dynamic road layers initialDistanceToRoads = $gisData$\Dist2ActiveRoad initialNearestRoadLoc = $gisData$\NearestRoadSegmentLoc Model Dimensions: Ecoregion Model.sel SimPriority Low Priority // Set low priority for simulation engine // Set up display Minimize Static Minimize Initial State Tile // Move to output folder cwd output // create output directories mkdir cell mkdir cellhd // Base Time step (in years) BaseTimestep = 10 // Make it spatial: 10-100ha blocks BlockSizeMin = 500 BlockSizeMax = 1000 // No explicit adjacency GreenupYears = 0 // Use area based harvesting (default) AreaBasedHarvest = TRUE // Make it aspatial: No road access restrictions MinNoEffectDist2Rd = 200 MaxDist2Road = 2000 SimStart 300 1

Expt1.scn
Seles Scenario $gisData$ = ..\gisData\cell Ecoregion = $gisData$\ecoregion Elevation = $gisData$\elevation Aspect = $gisData$\aspect ForestStatus = $gisData$\ForestStatus StreamClass = $gisData$\water initialAge = $gisData$\StandAgeEst initialSpp = $gisData$\DomSpecies SiteType = $gisData$\SiteEcoregion initialDistanceToRoads = $gisData$\Dist2ActiveRoad initialNearestRoadLoc = $gisData$\NearestRoadSegmentLoc Model Dimensions: Ecoregion Model.sel SimPriority Low Priority // Set low priority for simulation engine // Set up display Minimize Static Minimize Initial State Tile // Move to output folder cwd outputExpt // Base Time step (in years) BaseTimestep = 10 // No explicit adjacency GreenupYears = 0 // Use area based harvesting (default) AreaBasedHarvest = TRUE // Make it aspatial: No road access restrictions MinNoEffectDist2Rd = 200 MaxDist2Road = 2000 for ($min$ = 0:100,50) BlockSizeMin = $min$ for ($max$ = 100:200,100) BlockSizeMax = $max$ cwd max$max$_min$min$ // create output directories system "mkdir cell" system "mkdir cellhd" SimStart 50 1 cwd .. end end

Model.sel
Seles Model Time Units: Year Century 100 260 Landscape Events: // Disturbance event Fire.lse DEBUG // Succession events Succession.lse DEBUG // Management events GrowingStock.lse DEBUG Logging.lse DEBUG // Reporting events ReportResults.lse DEBUG // Legends: defined constants Legends: ////////////////////////////////////////  ForestStatusLegend = "..\gisData\cats\ForestStatus" // For succession SppLegend = "..\gisData\cats\DomSpecies" // For volume yields SiteTypeLegend = ..\gisData\cats\SiteEcoregion // For potential treatment types TreatmentTypeLegend = "cats\PotentialTreatmentType" // For limiting constraints ConstraintTypes = {0:lTotal, 1:lTHLB, 2:lAvailable, 3:lMHA, 4:lDist, 5:lAdjGreenup, 6:lAvailable } Global Constants: ///////////////////////////////////////////////////////////////  // General constants and parameters CellWidth = CELL WIDTH(Ecoregion) HaPerCell = (CellWidth^2) / 10000 MaxStandAge = MAX(initialAge) MaxVolPerHa = 350 MaxTimeSinceDisturbance = 200 NumSppCodes = ROWS(SppLegend) - 1 // Management status: Cell availability and potential treatment types NumTreatmentTypes = ROWS(TreatmentTypeLegend) - 1 // For planting and MHA SiteInfo = ".\inputFiles\SiteInfo.txt" rMHA = 0 rCMAI = 1 NumHarvestableCells = NUM>1(ForestStatus) // Need to link spatial constants with an input raster view Spatial Constants: StudyArea = Ecoregion ///////////////////////////  // Landscape physiography ///////////////////////////  Elevation Aspect StreamClass ///////////////////////////  // Management information (static) ///////////////////////////  ForestStatus SiteType // Need to link spatial variables with an output raster view, and // specify initial conditions (raster, 0, or constant) and output raster bounds Spatial Variables: ///////////////////////////////////////////////////////////////  // Dynamic stand information StandAge[MaxStandAge] <= initialAge Spp[NumSppCodes] <= initialSpp OldForest[1] ///////////////////////////  // Management information (dynamic) ///////////////////////////  VolPerHa[100*MaxVolPerHa] SalvageableVolPerHa[100*MaxVolPerHa] // For stratifying the landscape into availability and potential treatment types PotentialTreatmentType[NumTreatmentTypes] // Road information DistanceFromRoad = DistanceToRoads <= initialDistanceToRoads NearestRoadLoc[NumCells] <= initialNearestRoadLoc ///////////////////////////////////////////////////////////////  // Time since disturbance information Burnt[MaxTimeSinceDisturbance] Logged[MaxTimeSinceDisturbance] VolHarvTotal[10*100*MaxVolPerHa*HaPerCell + 1] // 2 decimal places Global Constants: // Volume tables VolTable = ".\inputFiles\Volume.txt" VolCols = COLS(VolTable) MaxDecade = VolCols - 1 // Minimum yields MinYieldGreenPerHa = 100 // For Constraints NumConstraintTypes = ROWS(ConstraintTypes) // For age class distribution AgeClassSize = 10 MaxAgeClass = 16 rContributing = 0 rNonContributing = 1 // For block size dist BlockClassSize = 10 MaxBlockSize = 100 NumBlockClasses = (MaxBlockSize / BlockClassSize) + 1 // For stream crossings NumStreamClasses = MAX(StreamClass) Global Variables: ///////////////////////////////////////////////////////////////  // Parameters BaseTimestep = 1 // Base Time step (in years) ///////////////////////////  // Fire related parameters ///////////////////////////  ReturnInterval = 50 MeanFiresPerYear = 1 ///////////////////////////  // Logging related parameters ///////////////////////////  HarvestRotation = 100  // Years, for area-based harvest AAC = 100000 // m3 for volume-based harvest // Adjacency GreenupYears = 0 //Default: 0 no explicit adj.   BlockSizeMin = 5 BlockSizeMax = 40 // Whether this is volume based or area based harvesting AreaBasedHarvest = FALSE // For controlling salvage: to turn off salvage, set this to larger than MaxVolPerCell MinYieldSalvagePerHa = 50 //5000000 SalvagePreference = 1 // Range between start of deline in prob.  MinNoEffectDist2Rd = 500 MaxDist2Road = 2000 ///////////////////////////////////////////////////////////////  // Tracking Year = 0 ForestSize = 0 THLBSize = 0 AreaBurned = 0    // Area burned in period AreaHarvested = 0 // Area (ha) harvested in period VolHarvested = 0  // Volume (m3) harvested in period VolumeLost = 0    // Volume lost due to decay VolumeSalvaged = 0 OFF SumOfAgeHarvested = 0 // Total sum of stand ages/ha harvested in period (divied by area to get mean) AnnualHarvestTarget = 0          // Target volume (m3) to harvest in period BlockSizeDist[NumBlockClasses] = 0 // block size dist BlockSizeArea[NumBlockClasses] = 0 AreaByAgeClass[2,MaxAgeClass+1] = 0 // age clas dist. THBL vs Non-contributing // For tracking roads built and stream crossings KmRoadsBuilt = 0 NumStreamCrossings[NumStreamClasses+1] = 0 ///////////////////////////////////////////////////////////////  // Frequency (in years) of outputting spatial layers ReportingInterval = 10 Output Frequency: 1 Output Model Frequency: StandAge Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative OldForest Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative // Spp Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative DistanceToRoads Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative VolHarvTotal Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative VolPerHa Freq: ReportingInterval Directory: ".\cell" Type: GRASS COMPRESSED Relative

Fire.lse
// Empirical fire model for the Labrador District 19a "study copy" LSEVENT: Fire DEFINITIONS GLOBAL CONSTANT: HaPerCell, MaxTimeSinceDisturbance GLOBAL VARIABLE: BaseTimestep // Base Time step (in years) GLOBAL VARIABLE: ReturnInterval, MeanFiresPerYear // fire return interval LAYER: StudyArea, StandAge, Spp LAYER: Elevation, Aspect    // elevation in metres; aspect 0 to 360 + flat (361) LAYER: Burnt  // Show areas that were recently burned LAYER: ForestStatus GLOBAL CONSTANT: Harvestable LAYER: SalvageableVolPerHa, VolPerHa GLOBAL VARIABLE: ForestSize, AreaBurned CLUSTER VARIABLE: FireExtent ENDDEF // Run after succession RETURNTIME RETURNTIME = IF Time EQ 0 THEN 0.5 ELSE BaseTimestep // Update time since disturbance information // Need to do this here, since Burnt may include non-forested cells OVER REGION WHOLE MAP DECISION Burnt > 0 Burnt = MAX(0,Burnt - BaseTimestep) ENDFN AreaBurned = 0 ENDRT // Only allow ignitions in forested cells EVENTLOCATION STATIC REGION WHOLE MAP DECISION (StudyArea > 0) AND (Spp > 0) ENDEL // Draw the number of openings from an exponential distribution NUMCLUSTERS meanClusters = MeanFiresPerYear * BaseTimestep nClusters = FLOOR(NEGEXP(meanClusters)) NUMCLUSTERS = nClusters ENDNC PROBINIT // Depends on aspect pAspect = 0.35*0.5*(SIN((Aspect-67.5))+1) PROBINIT = pAspect // For each opening, select an opening size meanExtent = (ForestSize/(MeanFiresPerYear*ReturnInterval)) FireExtent = ROUND(NEGEXP(meanExtent)) ENDPI TRANSITIONS // Continue if there is still extent to be burned // AND if the stand didn't burn during this event already TRANSITIONS = (FireExtent > 0) AND (Burnt < MaxTimeSinceDisturbance) IF Spp > 0 // Decrement the number of cells remaining to burn for this opening // Only do this for forested cells. NonForested cells will just pass on the fire FireExtent = FireExtent - HaPerCell AreaBurned = AreaBurned + HaPerCell ENDFN // Set the stand age StandAge = 0 Burnt = MaxTimeSinceDisturbance // Assume any unsalvaged yield is burned IF (ForestStatus EQ Harvestable) SalvageableVolPerHa = ROUND(0.8 * VolPerHa) ENDFN VolPerHa = 0 ENDTR // Spread timestep: time is irrelevant for this empirical model. SPREADTIMESTEP = 0 // Spread to the eight neighbours SPREADLOCATION REGION CENTRED(1, 1.5) DECISION (StudyArea > 0) ENDSL // Pick number of spread openings to be based on the // convolution index, but force at least one NUMRECIPIENTS // Mid-Complex shapes. This is the number of neighbours to // spread to from a burning cell. Lower mean gives more complex shapes. NUMRECIPIENTS = CLAMP(NORMAL(2,1), 1, 8) ENDNR // Prefer southwest aspects and uphill SPREADPROB pElevation = IF (SOURCE Elevation) < Elevation THEN 5 ELSE IF (SOURCE Elevation) > Elevation THEN 0.5 ELSE 1 pAspect = 0.35*0.5*(SIN((Aspect-67.5))+1) SPREADPROB = (pElevation^0.5) * (pAspect^0.5) ENDSP

GrowingStock.lse
// Landscape event definition for growing stock/cell availability in the Labrador District 19a "study copy" LSEVENT: GROWING_STOCK DEFINITIONS GLOBAL CONSTANT: HaPerCell /****************************************/  // Parameters GLOBAL VARIABLE: BaseTimestep // Base Time step (in years) // Length of time for greenup constraints (default: 15 years) GLOBAL VARIABLE: GreenupYears // The maximum distance from road for harvesting is assumed to be 2000m. Set higher (along with MinNoEffectDist2Rd) to remove access restrictions. GLOBAL VARIABLE: MaxDist2Road // Minimum yield that will be harvested in m3 per ha  GLOBAL CONSTANT: MinYieldGreenPerHa GLOBAL VARIABLE: MinYieldSalvagePerHa GLOBAL CONSTANT: SiteInfo[], rMHA // Volume curves GLOBAL CONSTANT: VolTable[], MaxDecade LAYER: StudyArea, StandAge LAYER: ForestStatus GLOBAL CONSTANT: Harvestable LAYER: SiteType, DistanceFromRoad LAYER: VolPerHa, SalvageableVolPerHa // in 100ths of a cubic metre GLOBAL CONSTANT: MaxVolPerHa // Treatment tracking LAYER: PotentialTreatmentType GLOBAL CONSTANT: Unavailable, GreenBlock, SalvageBlock // Limiting Constraints GLOBAL CONSTANT: ConstraintTypes[], NumConstraintTypes, lTotal, lTHLB, lMHA, lDist, lAdjGreenup, lAvailable LOCAL: Vol[NumConstraintTypes], Area[NumConstraintTypes] LOCAL: AreaLimited[NumConstraintTypes], AreaLimited2[NumConstraintTypes] OUTPUT VARIABLE: limitingConstraints = limitingConstraints.txt OUTPUT VARIABLE: growingStock = growingStock.txt ENDDEF RETURNTIME // Output report at the end of the event (which is captured just before the next instance is scheduled) IF (Time > 0) OVER INDEX SEQUENCE(0, NumConstraintTypes-1) OUTPUT RECORD(growingStock) Run Year: FLOOR(Time) CategoryId: Index Category: $ConstraintTypes[Index] Vol: Vol[Index] Area: Area[Index] ENDFN OUTPUT RECORD(limitingConstraints) Run Year: FLOOR(Time) CategoryId: Index Category: $ConstraintTypes[Index] NetArea: AreaLimited[Index] TotalArea: AreaLimited2[Index] ENDFN ENDFN ENDFN RETURNTIME = IF Time EQ 0 THEN 0.01 ELSE BaseTimestep // Clear tracking variables Vol [=] 0 Area [=] 0 AreaLimited [=] 0 AreaLimited2 [=] 0 ENDRT // Only process productive forest cells EVENTLOCATION STATIC REGION WHOLE MAP DECISION (StudyArea > 0) AND (ForestStatus > 0) ENDEL PROBINIT PROBINIT = 1 // Process all cells in productive forest. PotentialTreatmentType = Unavailable decade = CLAMP(StandAge/10,0,MaxDecade) l = FLOOR(decade) u = CEILING(decade) pLower = u - decade TotalVolPerHa = pLower * VolTable[SiteType, l] + (1-pLower) * VolTable[SiteType, u]   VolPerHa = CLAMP(ROUND(TotalVolPerHa*100),0,MaxVolPerHa*100) Yield = TotalVolPerHa * HaPerCell // Step through zones and stratifications to determine availability // areas and volumes Vol[lTotal] = Vol[lTotal] + Yield Area[lTotal] = Area[lTotal] + HaPerCell // Only apply remaining constraints in THLB IF (ForestStatus EQ Harvestable) lConstraint = lAvailable harvestArea = HaPerCell Vol[lTHLB] = Vol[lTHLB] + Yield Area[lTHLB] = Area[lTHLB] + harvestArea // First check minimum harvest age isMature = (StandAge >= SiteInfo[SiteType, rMHA]) IF isMature Vol[lMHA] = Vol[lMHA] + Yield Area[lMHA] = Area[lMHA] + harvestArea ELSE lConstraint = lMHA AreaLimited2[lMHA] = AreaLimited2[lMHA] + harvestArea ENDFN // Accessible by road isAccessible = (DistanceFromRoad <= MaxDist2Road) IF isAccessible Vol[lDist] = Vol[lDist] + Yield Area[lDist] = Area[lDist] + harvestArea ELSE lConstraint = IF (lConstraint NEQ lAvailable) THEN lConstraint ELSE lDist AreaLimited2[lDist] = AreaLimited2[lDist] + harvestArea * (lConstraint NEQ lMHA) ENDFN // Adjacency and greenup inBuffer = FALSE IF (GreenupYears > 0) OVER REGION CENTRED(1, 1.5, EUCLIDEAN) inBuffer = inBuffer OR (StandAge <= GreenupYears) ENDFN ENDFN IF !inBuffer Vol[lAdjGreenup] = Vol[lAdjGreenup] + Yield Area[lAdjGreenup] = Area[lAdjGreenup] + harvestArea ELSE lConstraint = IF (lConstraint NEQ lAvailable) THEN lConstraint ELSE lAdjGreenup AreaLimited2[lAdjGreenup] = AreaLimited2[lAdjGreenup] + harvestArea * (lConstraint NEQ lMHA) ENDFN // Finally: we know if this cell is available IF (lConstraint EQ lAvailable) Vol[lAvailable] = Vol[lAvailable] + Yield Area[lAvailable] = Area[lAvailable] + harvestArea ENDFN PotentialTreatmentType = IF (lConstraint EQ lAvailable) THEN GreenBlock ELSE IF ((SalvageableVolPerHa/100) >= MinYieldSalvagePerHa*HaPerCell) AND (DistanceFromRoad <= MaxDist2Road) THEN SalvageBlock ELSE Unavailable AreaLimited[lConstraint] = AreaLimited[lConstraint] + harvestArea ENDFN ENDPI

Logging.lse
// Landscape event definition for logging in the Labrador District 19a "study copy" // It is driven either by the input variable HarvestRotation, which specifies the area // to cut per year or by the input variable AAC, which specifies the volume of wood // to cut per year. Enough cutblocks are laid down to cut the required target LSEVENT: LOGGING DEFINITIONS GLOBAL CONSTANT: HaPerCell, CellWidth /****************************************/  // Parameters GLOBAL VARIABLE: BaseTimestep    // Base Time step (in years) // If TRUE use area vs. volume based harvest. The first is specified as a rotation, the second as a target in m3   GLOBAL VARIABLE: AreaBasedHarvest, HarvestRotation, AAC GLOBAL VARIABLE: SalvagePreference // Preference for salvage vs. green blocks // Length of time for greenup constraints (default: 15 years) GLOBAL VARIABLE: GreenupYears GLOBAL VARIABLE: BlockSizeMin, BlockSizeMax GLOBAL CONSTANT: SiteInfo[], rMHA // Range between start of deline in prob. // The maximum distance from road for harvesting is assumed to be 2000m. Set higher (along with MinNoEffectDist2Rd) to remove access restrictions. GLOBAL VARIABLE: MinNoEffectDist2Rd, MaxDist2Road, DistanceEffectSlope /****************************************/  LAYER: StudyArea, StandAge LAYER: ForestStatus GLOBAL CONSTANT: Harvestable LAYER: SiteType LAYER: Logged // Time since logging is MaxTimeSinceDisturbance - Logged GLOBAL CONSTANT: MaxTimeSinceDisturbance LAYER: DistanceFromRoad, NearestRoadLoc LAYER: VolPerHa, SalvageableVolPerHa, VolHarvTotal // in 100ths of a cubic metre LAYER: PotentialTreatmentType GLOBAL CONSTANT: Unavailable, GreenBlock, SalvageBlock GLOBAL VARIABLE: THLBSize, TotalHarvestTarget EVENT VARIABLE: EventExtent, HarvestTarget CLUSTER VARIABLE: BlockExtent, SalvageLogging CELL VARIABLE: NearestRoadLocation, NewDistanceFromRoad, StillLogging // Amount and area cut in year GLOBAL VARIABLE: VolHarvested, AreaHarvested, AnnualHarvestTarget, SumOfAgeHarvested, VolumeSalvaged GLOBAL VARIABLE: Year // Block size class distribution reporting GLOBAL CONSTANT: BlockClassSize, NumBlockClasses GLOBAL VARIABLE: BlockSizeDist[], BlockSizeArea[] CLUSTER VARIABLE: currBlockSize // For tracking stream crossings, roads built and stems removed LAYER: StreamClass GLOBAL VARIABLE: KmRoadsBuilt, NumStreamCrossings[] GLOBAL CONSTANT: NumStreamClasses ENDDEF INITIALSTATE DistanceEffectSlope = IF (MaxDist2Road > MinNoEffectDist2Rd) THEN 1 / (MinNoEffectDist2Rd - MaxDist2Road) ELSE 0 INITIALSTATE = 1 ENDIS // Process after fire RETURNTIME RETURNTIME = IF Time EQ 0 THEN 0.9 ELSE BaseTimestep Year = FLOOR(Time) // Reset annual harvest stats VolHarvested = 0 AreaHarvested = 0 SumOfAgeHarvested = 0 VolumeSalvaged = 0 KmRoadsBuilt = 0 NumStreamCrossings [=] 0 ENDRT // Only log in forested cells in the harvestable landbase EVENTLOCATION STATIC REGION WHOLE MAP DECISION ForestStatus EQ Harvestable ENDEL NUMCLUSTERS // Determine harvest target for entire time period TotalHarvestTarget = IF AreaBasedHarvest THEN (THLBSize/HarvestRotation) * BaseTimestep ELSE AAC * BaseTimestep AnnualHarvestTarget = TotalHarvestTarget HarvestTarget = TotalHarvestTarget // Keep logging until target is met or there are not more available cells NUMCLUSTERS = WHILE HarvestTarget > 0 // Check if this is salvage SalvageLogging = (PotentialTreatmentType EQ SalvageBlock) // Select from a uniform distribution BlockExtent = ROUND(UNIFORM(BlockSizeMin/HaPerCell, BlockSizeMax/HaPerCell)) currBlockSize = 0 ENDNC // Assumes that the growing stock sub-model has run to determine which cells // are available for harvest at the start of this time period PROBINIT ORDERED Logged = MAX(Logged-BaseTimestep, 0) MinHarvestAge = SiteInfo[SiteType, rMHA] // Probability based on distance pDistance = IF DistanceFromRoad <= MinNoEffectDist2Rd THEN 1 ELSE MAX(0,1 + (DistanceFromRoad- MinNoEffectDist2Rd)*DistanceEffectSlope) // Access is the only hard constraint for salvage pSalvageBlock = IF PotentialTreatmentType EQ SalvageBlock THEN SalvagePreference * (SalvageableVolPerHa/100) ELSE 0 // NOTE: These probabilities are computed at the start of the logging year and do not reflect // changes as harvesting progresses. Hence, to avoid violations, a check is made again in TRANSITIONS p = PRODUCT PotentialTreatmentType > Unavailable // If this is a naturally disturbed cell with salvageable timber... IF PotentialTreatmentType EQ SalvageBlock THEN pSalvageBlock // Otherwise, this is a green cell... // Relative oldest first ELSE (StandAge - MinHarvestAge + 1) pDistance ENDFN PROBINIT = p  NearestRoadLocation = NearestRoadLoc NewDistanceFromRoad = DistanceFromRoad StillLogging = (HarvestTarget > 0) ENDPI TRANSITIONS // FIRST: Update "still logging" flag // Check if we should still be logging. This must be done to avoid logging a cell that is in a  constraining zone. // This may occur during spread, and also for the first cell of a new block (in which case,we must also create an empty opening) IF StillLogging // There is still extent left in the block and area/volume to harvest isAvailable = (BlockExtent >= 1) AND (HarvestTarget > 0) IF SalvageLogging isAvailable = isAvailable AND (StandAge < 10) AND (PotentialTreatmentType EQ SalvageBlock) ELSE // Min. harvest age isAvailable = isAvailable AND (StandAge >= SiteInfo[SiteType, rMHA]) // Adjacency inBuffer = FALSE IF (GreenupYears > 0) OVER REGION CENTRED(1, 1.5, EUCLIDEAN) inBuffer = inBuffer OR (StandAge <= GreenupYears) ENDFN ENDFN isAvailable = isAvailable AND (!inBuffer) ENDFN // Finally: we know if this cell is available StillLogging = isAvailable ENDFN // Make a transition only if we are still logging or updating roads TRANSITIONS = OR                   StillLogging (currBlockSize > 0) AND (NewDistanceFromRoad < DistanceFromRoad) AND (NewDistanceFromRoad <= MaxDist2Road) //AND (DistanceFromRoad > MinNoEffectDist2Rd) ENDFN // If we are still logging, change logged cells IF StillLogging MinHarvestAge = SiteInfo[SiteType, rMHA] Yield = (SalvageableVolPerHa/100 + VolPerHa/100) * HaPerCell // Decrement the number of cells remaining to log for this block BlockExtent = BlockExtent - 1 currBlockSize = currBlockSize + HaPerCell IF AreaBasedHarvest HarvestTarget = HarvestTarget - HaPerCell TotalHarvestTarget = TotalHarvestTarget - HaPerCell ELSE HarvestTarget = HarvestTarget - Yield TotalHarvestTarget = TotalHarvestTarget - Yield ENDFN // Increment the yearly volume and areas harvested VolHarvested = VolHarvested + Yield AreaHarvested = AreaHarvested + HaPerCell VolHarvTotal = VolHarvTotal + ROUND(Yield * 100) SumOfAgeHarvested = SumOfAgeHarvested + StandAge * HaPerCell // Reset the stand StandAge = 0 Logged = MaxTimeSinceDisturbance // Update salvaged volume IF SalvageLogging VolumeSalvaged = VolumeSalvaged + (SalvageableVolPerHa/100)*HaPerCell ENDFN VolPerHa = 0 SalvageableVolPerHa = 0 // Reset the salvage volume // If this is an initiating cell: create a spur IF (currBlockSize EQ HaPerCell) AND (DistanceFromRoad > (MaxDist2Road/4)) KmRoadsBuilt = KmRoadsBuilt + DISTANCE(Location, NearestRoadLoc) * (CellWidth/1000) OnStreamCrossing = FALSE NearestRoadLocation = Location NewDistanceFromRoad = 0 OVER REGION VECTOR(Location, NearestRoadLoc) // Count stream crossings (only once per stream) IF (StreamClass > 0) IF (OnStreamCrossing EQ FALSE) NumStreamCrossings[StreamClass] = NumStreamCrossings[StreamClass] + 1 OnStreamCrossing = TRUE ENDFN ELSE OnStreamCrossing = FALSE ENDFN NearestRoadLoc = Location DistanceFromRoad = 0 ENDFN ENDFN ENDFN // Update the Distance from Road layer DistanceFromRoad = NewDistanceFromRoad NearestRoadLoc = NearestRoadLocation ENDTR // When a block is finished ENDCLUSTER // If at least one cell was harvested IF currBlockSize > 0 // Update block size distribution blockClass = MIN(FLOOR(currBlockSize / BlockClassSize), NumBlockClasses-1) BlockSizeDist[blockClass] = BlockSizeDist[blockClass] + 1 BlockSizeArea[blockClass] = BlockSizeArea[blockClass] + currBlockSize ENDFN ENDCLUSTER = TRUE ENDEC // Use immediate spread so that an entire block is processed before the next is initiated SPREADTIMESTEP = -1 // Spread to cardinal neighbours SPREADLOCATION REGION CENTRED(1,1) DECISION (StudyArea > 0) // AND (ForestStatus EQ Harvestable) ENDSL SPREADPROB // Keep spreading: transitions will stop when the cutblock is finished, and // the distance to road update is finished SPREADPROB = 1 // Transfer "still logging" status from source cell and we are in the THLB StillLogging = (SOURCE StillLogging) AND (ForestStatus EQ Harvestable) // If we are on a new road, refresh nearest road location variable IF (NearestRoadLoc EQ Location) //OR (Logged EQ MaxTimeSinceDisturbance) NearestRoadLocation = Location NewDistanceFromRoad = 0 ELSE NearestRoadLocation = SOURCE NearestRoadLocation NewDistanceFromRoad = CellWidth * DISTANCE(Location, NearestRoadLocation) ENDFN ENDSP

ReportResults.lse
// Landscape event definition for reporting in the Labrador FMD19a "study copy" LSEVENT: REPORT_RESULTS DEFINITIONS GLOBAL CONSTANT: HaPerCell /****************************************/  // Parameters GLOBAL VARIABLE: BaseTimestep // Base Time step (in years) GLOBAL VARIABLE: AreaBasedHarvest /****************************************/  // Stats by harvesting GLOBAL VARIABLE: VolHarvested, AreaHarvested, AreaBurned, VolumeLost GLOBAL VARIABLE: VolumeSalvaged, AnnualHarvestTarget, SumOfAgeHarvested GLOBAL VARIABLE: KmRoadsBuilt // For tracking stream crossings GLOBAL VARIABLE: NumStreamCrossings[] GLOBAL CONSTANT: NumStreamClasses // Block size class distribution reporting GLOBAL CONSTANT: BlockClassSize, NumBlockClasses GLOBAL VARIABLE: BlockSizeDist[], BlockSizeArea[] // Age class dist. reporting GLOBAL CONSTANT: MaxAgeClass, rContributing, rNonContributing GLOBAL VARIABLE: AreaByAgeClass[] MIN OUTPUT VARIABLE: ReportFile = "annualRecord.txt" MIN OUTPUT VARIABLE: BlockDist = BlockSizeDist.txt MIN OUTPUT VARIABLE: ageClassFile = ageClass.txt ENDDEF // Evaluate yearly (or by "BaseTimestep" years) RETURNTIME RETURNTIME = IF Time EQ 0 THEN 0.02 ELSE BaseTimestep MeanVolPerHa = VolHarvested/AreaHarvested OUTPUT RECORD(ReportFile) DECISION (Time > 1) Run Year: FLOOR(Time) AreaBurned: AreaBurned/BaseTimestep AreaHarvested: AreaHarvested/BaseTimestep VolHarvested: VolHarvested/BaseTimestep MeanVolPerHa PercentOfTarget: IF AreaBasedHarvest THEN ROUND(100*AreaHarvested/AnnualHarvestTarget) ELSE ROUND(100*VolHarvested/AnnualHarvestTarget) MeanHarvestAge: HaPerCell * SumOfAgeHarvested / AreaHarvested VolumeSalvaged:VolumeSalvaged/BaseTimestep VolumeLost: VolumeLost/BaseTimestep KmRoadsBuilt:KmRoadsBuilt/BaseTimestep StreamCrossings#1:NumStreamClasses: NumStreamCrossings[Index] ENDFN // Age class file OUTPUT RECORD(ageClassFile) DECISION (Time > 1) Run: Run Year: FLOOR(Time) C#0:MaxAgeClass: AreaByAgeClass[rContributing, Index] NC#0:MaxAgeClass: AreaByAgeClass[rNonContributing, Index] ENDFN // Block size dist IF (Time > 0) AND ((Time +BaseTimestep) > EndTime) n = 0 a = 0 OVER INDEX SEQUENCE(0,NumBlockClasses-1) n = n + BlockSizeDist[Index] a = a + BlockSizeArea[Index] ENDFN OVER INDEX SEQUENCE(0,NumBlockClasses-1) OUTPUT RECORD(BlockDist) Run: Run blockClass: Index blockSize: (Index + 0.5) * BlockClassSize NumBlocks: BlockSizeDist[Index] Area: BlockSizeArea[Index] pNumBlocks: BlockSizeDist[Index]/n pArea: BlockSizeArea[Index] / a       ENDFN ENDFN ENDFN ENDRT // Setting the number of clusters to zero will cause no cell initiations to occur NUMCLUSTERS = 0

Succession.lse
// Landscape event definition for succession in the Labrador District 19a "study copy" LSEVENT: SUCCESSION DEFINITIONS GLOBAL CONSTANT: MaxStandAge, HaPerCell /****************************************/  // Parameters GLOBAL VARIABLE: BaseTimestep  // Base Time step (in years) /****************************************/  GLOBAL VARIABLE: ForestSize, THLBSize LAYER: StudyArea, StandAge, OldForest LAYER: Spp LAYER: ForestStatus GLOBAL CONSTANT: Harvestable LAYER: SalvageableVolPerHa GLOBAL VARIABLE: VolumeLost // For age class info GLOBAL CONSTANT: AgeClassSize, MaxAgeClass, rContributing, rNonContributing GLOBAL VARIABLE: AreaByAgeClass[] ENDDEF INITIALSTATE // Set up the initial state OVER REGION WHOLE MAP DECISION (StudyArea > 0) AND (Spp > 0) // Old forest is over 150 years OldForest = (StandAge >= 150) ForestSize = ForestSize + HaPerCell IF ForestStatus EQ Harvestable THLBSize = THLBSize + HaPerCell ENDFN ENDFN INITIALSTATE = 1 ENDIS // Aging at start of year (before any other models), but not in first period RETURNTIME RETURNTIME = BaseTimestep VolumeLost = 0 AreaByAgeClass [=] 0 ENDRT // Only process forested cells EVENTLOCATION STATIC REGION WHOLE MAP DECISION (StudyArea > 0) AND (Spp > 0) ENDEL // In each forested cell TRANSITIONS // Increment the stand age in years StandAge = MIN(StandAge + BaseTimestep, MaxStandAge) // Old forest is over 150 years OldForest = (StandAge >= 150) // loss of salvageable wood (20%/year) IF SalvageableVolPerHa > 0 volLoss = CEIL(0.2 * SalvageableVolPerHa) SalvageableVolPerHa = SalvageableVolPerHa - volLoss VolumeLost = VolumeLost + volLoss ENDFN // Set up age class info ageClass = CLAMP(FLOOR(StandAge/AgeClassSize), 0, MaxAgeClass) status = IF (ForestStatus EQ Harvestable) THEN rContributing ELSE rNonContributing AreaByAgeClass[status, ageClass] = AreaByAgeClass[status, ageClass] + HaPerCell TRANSITIONS = TRUE // Markov chain succession... ENDTR

Suggested Experiments
To explore this cellular automata model further, try the following: