From 39f3e716dd8fbb203524d95e78a9129df0d8a057 Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Mon, 28 Nov 2022 11:07:51 +0000
Subject: [PATCH] Changes to timber yield calculation and optimisation.

---
 GAMS/IntExtOpt.gms                            |  70 ++++----
 src/ac/ed/lurg/ModelConfig.java               |  15 +-
 src/ac/ed/lurg/ModelMain.java                 |   3 +-
 src/ac/ed/lurg/Timestep.java                  |  14 --
 .../ed/lurg/country/AbstractCountryAgent.java |   2 +-
 src/ac/ed/lurg/country/CountryAgent.java      |  24 +--
 .../lurg/country/CountryBoundaryReader.java   |  15 +-
 .../lurg/country/gams/GamsCountryInput.java   |   6 +-
 .../lurg/country/gams/GamsLocationInput.java  |  11 +-
 .../country/gams/GamsLocationOptimiser.java   |  98 ++++++-----
 .../lurg/country/gams/GamsLocationOutput.java |  10 +-
 .../country/gams/GamsRasterOptimiser.java     | 162 +++++++++++++-----
 .../ed/lurg/demand/AbstractDemandManager.java |  38 ++--
 .../lurg/forestry/ForestryDataOutputer.java   |  16 +-
 src/ac/ed/lurg/forestry/ForestryPlanner.java  |   7 +
 src/ac/ed/lurg/forestry/WoodYieldData.java    |  38 ++--
 src/ac/ed/lurg/forestry/WoodYieldItem.java    | 159 +++++++----------
 .../ed/lurg/forestry/WoodYieldRasterSet.java  |  21 ++-
 src/ac/ed/lurg/forestry/WoodYieldReader.java  |  65 +++----
 .../ed/lurg/forestry/WoodYieldResponse.java   |  95 ++++++++++
 src/ac/ed/lurg/landuse/LandCoverItem.java     |  27 ++-
 src/ac/ed/lurg/landuse/LandUseItem.java       |  34 ++++
 src/ac/ed/lurg/types/LandCoverType.java       |  37 ++--
 23 files changed, 581 insertions(+), 386 deletions(-)
 create mode 100644 src/ac/ed/lurg/forestry/ForestryPlanner.java
 create mode 100644 src/ac/ed/lurg/forestry/WoodYieldResponse.java

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 730865bf..8a750b8c 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -14,10 +14,13 @@
  SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
  SET managed_forest(land_cover) / timberForest, carbonForest /;
  SET agriculture(land_cover) / cropland, pasture /;
+ SET wood_producing(land_cover) / timberForest, natural /;
  ALIAS (land_cover, land_cover_before);
  ALIAS (land_cover, land_cover_after);
 
  SET location;
+ SET year;
+
  PARAMETER suitableLandArea(location)            areas of land in Mha;
  PARAMETER previousCropArea(crop, location)      areas for previous timestep in Mha;
  PARAMETER previousFertIntensity(crop, location);
@@ -48,12 +51,11 @@
  PARAMETER subsidyRate(all_types)                 rates of subsidy compared to costs;
  
  PARAMETER previousLandCoverArea(land_cover, location)                      land cover area in Mha;
- PARAMETER carbonFluxRate(land_cover_before, land_cover_after, location)   carbon flux from LULUCF - tC-eq per ha;
  PARAMETER carbonCreditRate(land_cover_before, land_cover_after, location) potential carbon credits from LULUCF - tC-eq per ha;
- PARAMETER woodYieldRota(location)                                          wood yield from forest rotation - tC per ha;
- PARAMETER woodYieldLUC(land_cover_before, land_cover_after, location)      wood yield from land cover change;
- PARAMETER forestRotationPeriod(location)           timber forest rotation period - years;
+ PARAMETER woodYieldLUC(wood_producing, land_cover_after, location)      wood yield from land cover change;
  PARAMETER conversionCost(land_cover_before, land_cover_after, location)   cost of converting from one land cover to another - $1000 per ha;
+ PARAMETER woodYieldMax(location);
+ PARAMETER woodYieldParam(location);
 
  SCALAR meatEfficency                        efficiency of converting feed and pasture into animal products;
  SCALAR fertiliserUnitCost                   fert cost at max fert rate;
@@ -69,18 +71,17 @@
  SCALAR forestManagementCost    cost $1000 per ha;
  SCALAR vegClearingCostRate     cost of clearing vegetation $1000 per tC;
  
- SCALAR carbonHorizon           period over which carbon credits are calculated;
  SCALAR carbonForestMaxProportion maximum proportion of land cover as carbon forest;
 
 *$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
 $gdxin %gdxincname%
-$load location, suitableLandArea, demand
+$load location, year, suitableLandArea, demand
 $load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
 $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
 $load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxGrossLccRate, subsidyRate
 $load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
-$load previousLandCoverArea, carbonFluxRate, carbonCreditRate, woodYieldRota, woodYieldLUC, forestRotationPeriod, conversionCost
-$load forestManagementCost, vegClearingCostRate, carbonHorizon, carbonForestMaxProportion, maxFertChange, maxIrrigChange
+$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam
+$load forestManagementCost, vegClearingCostRate, carbonForestMaxProportion, maxFertChange, maxIrrigChange
 $gdxin
 
  SCALAR delta "use to smooth power function see 7.5 www.gams.com dd docs solversconopt.pdf" / 0.00000000001 /;
@@ -122,7 +123,7 @@ $gdxin
  PARAMETER animalProdCost(animal);
  animalProdCost('ruminants') = 0.28;
  animalProdCost('monogastrics') = 0.24;
- 
+
 * maxNetImport(import_types) = min(maxNetImport(import_types), demand(import_types));
 * minNetImport(import_types) = min(minNetImport(import_types), demand(import_types));
 
@@ -145,16 +146,15 @@ $gdxin
        landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
        totalConversionCost(location)      land cover conversion cost - $1000 per ha
 
-       woodSupplyRota
-       woodSupplyLUC
+       woodYieldRota(location)
+       rotationPeriod(location)
        forestryCost(location)
-       carbonFlux(location)               total carbon flux  - Mt C
        carbonCredits(location)
 *       A                                  "artificial variable for debugging https://www.gams.com/blog/2017/07/misbehaving-model-infeasible/"
        total_cost                         total cost of domestic supply including net imports;
 
- POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount, totalFeedDM,
-                   landCoverArea, landCoverChange, totalConversionCost, woodSupplyRota, woodSupplyLUC, animalProd;
+ POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount, totalFeedDM, 
+                   landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationPeriod;
 
 
 * POSITIVE VARIABLE A;
@@ -193,12 +193,12 @@ $gdxin
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
        CONVERSION_COST_CALC(location)                          cost of land cover conversion
 
-       WOOD_SUPPLY_ROTA_CALC
-       WOOD_SUPPLY_LUC_CALC
-       WOOD_SUPPLY_CONSTRAINT
+       WOOD_YIELD_CALC(location)
+       WOOD_DEMAND_CONSTRAINT
+       ROTATION_PERIOD_MIN_CONSTRAINT(location)
+       ROTATION_PERIOD_MAX_CONSTRAINT(location)
        FORESTRY_COST_CALC(location)
  
-       CARBON_FLUX_CALC(location)                         calc carbon flux
        CARBON_CREDIT_CALC(location)
        CARBON_CREDIT_CONSTRAINT
        CARBON_FOREST_CONSTRAINT
@@ -268,29 +268,31 @@ $gdxin
  PASTURE_LAND_COVER_CALC(location) .. landCoverArea('pasture', location) =E= cropArea('pasture', location);
 
  LAND_COVER_CHANGE_CALC(land_cover, location) .. landCoverArea(land_cover, location) =E= previousLandCoverArea(land_cover, location) +
-                                                                                         sum(land_cover_before, landCoverChange(land_cover_before, land_cover, location)) -
-                                                                                         sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location));
+                                                                                         sum((land_cover_before), landCoverChange(land_cover_before, land_cover, location)) -
+                                                                                         sum((land_cover_after), landCoverChange(land_cover, land_cover_after, location));
 
 
  LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) .. sum(land_cover_after, landCoverChange(land_cover, land_cover_after, location)) =E= previousLandCoverArea(land_cover, location);
 
  CONVERSION_COST_CALC(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after),
                                                                      landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after, location));
+                                                                     
+* MIN_LAND_COVER_AREA_CONSTRAINT(land_cover, location) .. landCoverArea(land_cover, location) =G= minLandCoverArea(land_cover, location);
+ 
 ************* Forestry ***********************************
-
- WOOD_SUPPLY_ROTA_CALC .. woodSupplyRota =E= sum(location, landCoverArea('timberForest', location) * woodYieldRota(location));
  
- WOOD_SUPPLY_LUC_CALC .. woodSupplyLUC =E= sum((land_cover_before, land_cover_after, location), landCoverChange(land_cover_before, land_cover_after, location) * woodYieldLUC(land_cover_before, land_cover_after, location));
+ WOOD_YIELD_CALC(location) .. woodYieldRota(location) =E= woodYieldMax(location) * (1 - exp(-woodYieldParam(location) * rotationPeriod(location))) / rotationPeriod(location);
  
- WOOD_SUPPLY_CONSTRAINT .. woodSupplyRota =G= demand('wood') + exportAmount('wood') - importAmount('wood');
+ ROTATION_PERIOD_MIN_CONSTRAINT(location) .. rotationPeriod(location) =G= 5;
+ ROTATION_PERIOD_MAX_CONSTRAINT(location) .. rotationPeriod(location) =L= 160;
  
- FORESTRY_COST_CALC(location) .. forestryCost(location) =E= landCoverArea('timberForest', location) * forestManagementCost / forestRotationPeriod(location) +
-                                                            landCoverArea('carbonForest', location) * forestManagementCost / carbonHorizon;
+ WOOD_DEMAND_CONSTRAINT .. demand('wood') + exportAmount('wood') - importAmount('wood') =L= sum(location, woodYieldRota(location) * landCoverArea('timberForest', location));
  
-*********** Carbon fluxes ***********************************
-
- CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
-                                                        
+ FORESTRY_COST_CALC(location) .. forestryCost(location) =E= sum(managed_forest, landCoverArea(managed_forest, location) * 0.05) +
+                                                            forestManagementCost / rotationPeriod(location);
+ 
+*********** Carbon ***********************************
+                                       
  CARBON_CREDIT_CALC(location) .. carbonCredits(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonCreditRate(land_cover_before, land_cover_after, location));
                                                         
  CARBON_CREDIT_CONSTRAINT .. sum(location, carbonCredits(location)) =G= demand('carbonCredits') + exportAmount('carbonCredits') - importAmount('carbonCredits');
@@ -324,6 +326,7 @@ $gdxin
  monogastricFeed.L(feed_crop) = previousMonogastricFeed(feed_crop);
  importAmount.L(all_types) = previousImportAmount(all_types);
  exportAmount.L(all_types) = previousExportAmount(all_types);
+ rotationPeriod.L(location) = 40;
 
 * LAND_USE.OptFile = 1;
 
@@ -346,8 +349,9 @@ $gdxin
  parameter netImportCost(all_types);
  parameter feedCostRate(feed_crop);
  parameter productionShock(all_types);
- scalar netCarbonFlux;
+ parameter carbonFlux(location);
  scalar netCarbonCredits;
+ scalar woodSupplyRota;
 
 * Production quantities based on smaller area (before unhandledCropArea adjustment applied)
  totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
@@ -367,8 +371,10 @@ $gdxin
 
  netImportAmount(import_types) = importAmount.l(import_types) - exportAmount.l(import_types);
  netImportCost(import_types) = importAmount.l(import_types) * importPrices(import_types) - exportAmount.l(import_types) * exportPrices(import_types);
- netCarbonFlux = SUM(location, carbonFlux.L(location));
+ 
  netCarbonCredits = SUM(location, carbonCredits.L(location));
+ 
+ woodSupplyRota = sum(location, woodYieldRota.L(location) * landCoverArea.L('timberForest', location));
 
  Scalar totalCostsLessLU;
 
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 5c5ba546..dec4e057 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -242,9 +242,9 @@ public class ModelConfig {
 	public static final double CELL_SIZE_Y = getDoubleProperty("CELL_SIZE_Y", CELL_SIZE_X);
 	public static final String SPATIAL_DIR_NAME = getProperty("SPATIAL_DIR_NAME", "halfdeg");
 	public static final String SPATIAL_DATA_DIR = getProperty("SPATIAL_DATA_DIR", DATA_DIR + File.separator + SPATIAL_DIR_NAME);
-	public static final String INITAL_LAND_COVER_FILENAME = getProperty("INITAL_LAND_COVER_FILENAME", "hurtt_2010.txt");
+	public static final String INITAL_LAND_COVER_FILENAME = getProperty("INITAL_LAND_COVER_FILENAME", "hilda_plus_2019.txt");
 	public static final String INITAL_LAND_COVER_FILE = SPATIAL_DATA_DIR + File.separator + INITAL_LAND_COVER_FILENAME;
-	public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.asc";
+	public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.csv";
 	public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost.asc";
 	public static final String IRRIGATION_CONSTRAINT_FILE = SPATIAL_DATA_DIR + File.separator + "blue_water_available_pseudoCRU_rcp8p5_2004_2013_grid_allhdyro_mm.txt";
 	public static final String FPU_BOUNDARIES_FILE = SPATIAL_DATA_DIR + File.separator + "FPU.asc";
@@ -270,7 +270,9 @@ public class ModelConfig {
 	public static final String WOOD_YIELD_NTRL_TO_PAST_FILENAME = getProperty("WOOD_YIELD_NTRL_TO_PAST_FILENAME", "landsymm_plutW_from_ntrl_to_past.dat");
 	public static final String WOOD_YIELD_NTRL_TO_CROP_FILENAME = getProperty("WOOD_YIELD_NTRL_TO_CROP_FILENAME", "landsymm_plutW_from_ntrl_to_crop.dat");
 	public static final String WOOD_YIELD_NTRL_TO_FORS_FILENAME = getProperty("WOOD_YIELD_NTRL_TO_FORS_FILENAME", "landsymm_plutW_from_ntrl_to_forC.dat");
-	public static final String WOOD_YIELD_ROTA_FILENAME = getProperty("WOOD_YIELD_ROTA_FILENAME", "landsymm_pcutW_sts_forC.dat");
+	public static final String WOOD_YIELD_ROTA_FORS_FILENAME = getProperty("WOOD_YIELD_ROTA_FORS_FILENAME", "landsymm_pcutW_sts_forC.dat");
+	public static final String WOOD_YIELD_ROTA_NTRL_FILENAME = getProperty("WOOD_YIELD_ROTA_NTRL_FILENAME", "landsymm_pcutW_sts_ntrl.dat");
+	public static final String WOOD_YIELD_PARAM_FILENAME = getProperty("WOOD_YIELD_PARAM_FILENAME", "pcutW_forC_param.out");
 	
 	public static final String C_FLUX_NTRL_TO_CROP_FILENAME = getProperty("C_FLUX_NTRL_TO_CROP_FILENAME", "landsymm_plutC_from_ntrl_to_crop.dat");
 	public static final String C_FLUX_NTRL_TO_PAST_FILENAME = getProperty("C_FLUX_NTRL_TO_PAST_FILENAME", "landsymm_plutC_from_ntrl_to_past.dat");
@@ -501,7 +503,7 @@ public class ModelConfig {
 	
 	// Forestry
 	public static final boolean IS_FORESTRY_ON = getBooleanProperty("IS_FORESTRY_ON", true);
-	public static final String WOOD_DEMAND_FILENAME = getProperty("WOOD_DEMAND_FILENAME", "wood_demand.csv");
+	public static final String WOOD_DEMAND_FILENAME = getProperty("WOOD_DEMAND_FILENAME", "wood_demand_adj.csv");
 	public static final String WOOD_DEMAND_FILE = getProperty("WOOD_DEMAND_FILE", DATA_DIR + File.separator + WOOD_DEMAND_FILENAME);
 	public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 1500.0); // MtC-eq
 	public static final double IND_ROUNDWOOD_DEMAND_ELASTICITY = getDoubleProperty("IND_ROUNDWOOD_DEMAND_ELASTICITY",  0.3123881); // fitted to FAO data
@@ -509,12 +511,13 @@ public class ModelConfig {
 	public static final double WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 3e-7); // m3 to MtC-eq p.16 [https://doi.org/10.5194/gmd-13-5425-2020] 
 	public static final double DISCOUNT_RATE = getDoubleProperty("DISCOUNT_RATE", 0.05); // for calculating rotation period
 	public static final double VEGETATION_CLEARING_COST = getDoubleProperty("VEGETATION_CLEARING_COST", 0.001); //$1000/tC
-	public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 5.0) : 0.0; // establishment, management etc. $1000/ha
+	public static final double FOREST_MANAGEMENT_COST = IS_FORESTRY_ON ? getDoubleProperty("FOREST_ESTABLISHMENT_COST", 3.0) : 0.0; // establishment, management etc. $1000/ha
 	public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/tC
-	public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.22) : 0.0; // $1000/tC-eq
+	public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.1) : 0.0; // $1000/tC-eq
 	public static final boolean ENABLE_VEGETATION_CLEARANCE_COST = getBooleanProperty("ENABLE_VEGETATION_CLEARANCE_COST", false); // calculates cost of clearing vegetation even if forestry off
 	public static final int MIN_FOREST_ROTATION_PERIOD = getIntProperty("MIN_FOREST_ROTATION_PERIOD", 5); // Minimum forest rotation period; Bauhus et al., 2010. Ecosystem goods and services from plantation forests.
 	public static final int MAX_FOREST_ROTATION_PERIOD = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years
+	public static final int LAND_AGE_GROUP_SIZE = getIntProperty("LAND_AGE_GROUP_SIZE", 5); // grouping size for land cover age
 	
 	// Carbon
 	public static final boolean IS_CARBON_ON = getBooleanProperty("IS_CARBON_ON", false);
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 83762f5d..c529dc40 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -158,10 +158,11 @@ public class ModelMain {
 		
 		getWoodYieldData(timestep);
 		getCarbonFluxData(timestep);
+		
+		LogWriter.println("Memory usage 1: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
 
 		ConversionCostReader conCostReader = new ConversionCostReader(timestep);
 		Map<LccKey, Double> conversionCosts = conCostReader.getConversionCosts();
-
 		
 		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, carbonFluxData, woodYieldData, conversionCosts);
 		
diff --git a/src/ac/ed/lurg/Timestep.java b/src/ac/ed/lurg/Timestep.java
index da4ed902..c92e3e41 100644
--- a/src/ac/ed/lurg/Timestep.java
+++ b/src/ac/ed/lurg/Timestep.java
@@ -104,18 +104,4 @@ public class Timestep implements Serializable {
 		return rootDir + File.separator + (endYear-ModelConfig.LPJG_TIMESTEP_SIZE+1) + "-" + endYear;
 	}
 	
-	public String getWoodAndCarbonYearSubDir(String rootDir) {
-		int startYear;
-		
-		if (ModelConfig.CHANGE_YIELD_DATA_YEAR) {
-			double currentYear = (float) getYear();
-			double timestepSize = (float) ModelConfig.CARBON_DATA_TIMESTEP_SIZE;
-			startYear = (int) (Math.floor(currentYear / timestepSize) * timestepSize) + 10;
-		}
-		else {
-			startYear = ModelConfig.BASE_YEAR;
-		}
-		
-		return rootDir + File.separator + startYear;
-	}
 }
diff --git a/src/ac/ed/lurg/country/AbstractCountryAgent.java b/src/ac/ed/lurg/country/AbstractCountryAgent.java
index 08d16b6c..d308ca22 100644
--- a/src/ac/ed/lurg/country/AbstractCountryAgent.java
+++ b/src/ac/ed/lurg/country/AbstractCountryAgent.java
@@ -31,7 +31,7 @@ public abstract class AbstractCountryAgent {
 	protected Timestep currentTimestep;
 	protected Map<CommodityType, Map<CropType, Double>> currentDemandFract;
 	protected double currentGen2EcDemand;
-	protected Map<WoodType, Double> currentWoodDemand;
+	protected Map<Integer, Map<WoodType, Double>> currentWoodDemand;
 	protected double currentCarbonDemand;
 	protected double exportTaxRate;
 	
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 5147bb42..59f94d17 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -147,8 +147,6 @@ public class CountryAgent extends AbstractCountryAgent {
 			
 			if (!ModelConfig.IS_CALIBRATION_RUN)
 				incrementLandCoverAge();
-			
-			doForestRotation(woodYieldData);
 
 			return result;
 		}
@@ -250,11 +248,12 @@ public class CountryAgent extends AbstractCountryAgent {
 			
 			if ((ModelConfig.IS_CALIBRATION_RUN || ModelConfig.USE_INITIAL_CROP_USAGE_DATA) && currentTimestep.isInitialTimestep()) {
 				double prodFract = woodUsage.getHarvest() / totalWoodProd; 
-				double potentialMaxYield = WoodYieldRasterSet.getMaxWoodHarvest(woodYieldData, previousGamsRasterOutput.getLandUses());
-				double potentialAdjYield = potentialMaxYield * prodFract; // adjust based on demand and cap at max 50% realised
-				double currentDemand = currentWoodDemand.get(woodType);
+				double potentialMaxYield = WoodYieldRasterSet.getPotentialTimberProduction(woodYieldData, previousGamsRasterOutput.getLandUses()) * 0.9;
+				double potentialAdjYield = potentialMaxYield * prodFract;
+				double currentDemand = currentWoodDemand.get(0).get(woodType);
 				double importsNeeded = currentDemand - potentialAdjYield;
-				baseTrade = importsNeeded > 0 ? importsNeeded : importsNeeded;
+				if (baseTrade > 0 & importsNeeded > baseTrade)
+					baseTrade = importsNeeded;
 			}
 			
 			double changeUp;
@@ -263,7 +262,7 @@ public class CountryAgent extends AbstractCountryAgent {
 			double maxOfProdOrSupply = woodUsage.getHarvest() + Math.max(baseTrade, 0);
 			
 			if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
-				changeUp = changeDown = 0;
+				changeUp = changeDown = 0.0;
 			} else {
 				changeUp = maxOfProdOrSupply * globalWoodPrice.getMaxImportChange();
 				changeDown = maxOfProdOrSupply * globalWoodPrice.getMaxExportChange();
@@ -343,17 +342,6 @@ public class CountryAgent extends AbstractCountryAgent {
 		}
 	}
 	
-	private void doForestRotation(RasterSet<WoodYieldItem> woodYieldData) {
-		for (Map.Entry<RasterKey, LandUseItem> entry : previousGamsRasterOutput.getLandUses().entrySet()) {
-			RasterKey key = entry.getKey();
-			LandUseItem luItem = entry.getValue();
-			WoodYieldItem wyItem = woodYieldData.get(key);
-			if (wyItem == null) 
-				continue;
-			luItem.doForestRotation(wyItem.getOptimalRotation());
-		}
-	}
-	
 	public CarbonUsageData getCarbonUsageData() {
 		return previousGamsRasterOutput.getCarbonUsageData();
 	}
diff --git a/src/ac/ed/lurg/country/CountryBoundaryReader.java b/src/ac/ed/lurg/country/CountryBoundaryReader.java
index e1eb38c5..22c8d9f6 100644
--- a/src/ac/ed/lurg/country/CountryBoundaryReader.java
+++ b/src/ac/ed/lurg/country/CountryBoundaryReader.java
@@ -1,17 +1,22 @@
 package ac.ed.lurg.country;
 
-import ac.sac.raster.AbstractRasterReader;
+import java.util.Map;
+
+import ac.sac.raster.AbstractTabularRasterReader;
+import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterSet;
 
-public class CountryBoundaryReader extends AbstractRasterReader<CountryBoundaryItem> {	
+public class CountryBoundaryReader extends AbstractTabularRasterReader<CountryBoundaryItem> {	
+	
+	private static int MIN_COLS = 2;
 	
 	public CountryBoundaryReader (RasterSet<CountryBoundaryItem> dataset) {
-		super(dataset);
+		super(",", MIN_COLS, dataset);
 	}
 
 	@Override
-	public void setData(CountryBoundaryItem item, String token) {
-		int isoCode = Integer.parseInt(token);
+	public void setData(RasterKey key, CountryBoundaryItem item, Map<String, Double> rowValues) {
+		int isoCode = (int) getValueForCol(rowValues, "FaoCode");
 		SingleCountry c = CountryManager.getForCode(isoCode);
 		item.setCountry(c);
 	}
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 4b5971c3..f621e4fa 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -26,7 +26,7 @@ public class GamsCountryInput {
 	private double carbonDemand;
 	private CountryPrice carbonPrice;
 	private TradeConstraint carbonTradeConstraint;
-	private Map<WoodType, Double> woodDemand;
+	private Map<Integer, Map<WoodType, Double>> woodDemand;
 	private CountryPrice woodPrice;
 	private Map<WoodType, TradeConstraint> woodTradeConstraints;
 	private Map<WoodType, WoodUsageData> previousWoodUsageData;
@@ -34,7 +34,7 @@ public class GamsCountryInput {
 	public GamsCountryInput(CompositeCountry country, Map<CommodityType, Double> projectedDemand, Map<CropType, CountryPrice> countryPrices, 
 			Map<CropType, TradeConstraint> importConstraints, Map<CropType, CropUsageData> previousCropUsageData, 
 			Map<CommodityType, Map<CropType, Double>> minDemandFracts, Map<CropType, Double> subsidyRates, double secondGenBioenergyDemand,
-			double carbonDemand, CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint, Map<WoodType, Double> woodDemand, 
+			double carbonDemand, CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint, Map<Integer, Map<WoodType, Double>> woodDemand, 
 			CountryPrice woodPrice, Map<WoodType, TradeConstraint> woodTradeConstraints, Map<WoodType, WoodUsageData> previousWoodUsageData) {
 		super();
 		this.country = country;
@@ -120,7 +120,7 @@ public class GamsCountryInput {
 		return carbonTradeConstraint;
 	}
 
-	public Map<WoodType, Double> getWoodDemand() {
+	public Map<Integer, Map<WoodType, Double>> getWoodDemand() {
 		return woodDemand;
 	}
 
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index 211b4974..8dfe79e3 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -20,12 +20,11 @@ public class GamsLocationInput {
 	private Map<Integer, ? extends WoodYieldData> woodYields;
 	private Map<LccKey, Double> conversionCosts;
 	private GamsCountryInput countryInput;
-	private Map<Integer, Double> rotationPeriods;
 	
 	public GamsLocationInput(Timestep timestep, Map<Integer, ? extends YieldResponsesItem> yields, Map<Integer, ? extends LandUseItem> previousLandUse,
 			Map<Integer, ? extends IrrigationItem> irrigationCosts, Map<Integer, ? extends CarbonFluxItem> carbonFluxes, 
 			Map<Integer, ? extends WoodYieldData> woodYields,
-			Map<LccKey, Double> conversionCosts, GamsCountryInput countryInput, Map<Integer, Double> rotationPeriods) {
+			Map<LccKey, Double> conversionCosts, GamsCountryInput countryInput) {
 		super();
 		this.timestep = timestep;
 		this.yields = yields;
@@ -34,8 +33,7 @@ public class GamsLocationInput {
 		this.carbonFluxes = carbonFluxes;
 		this.woodYields = woodYields;
 		this.conversionCosts = conversionCosts;
-		this.countryInput = countryInput;	
-		this.rotationPeriods = rotationPeriods;
+		this.countryInput = countryInput;
 	}
 		
 	public Map<Integer, ? extends YieldResponsesItem> getYields() {
@@ -65,12 +63,9 @@ public class GamsLocationInput {
 	public GamsCountryInput getCountryInput() {
 		return countryInput;
 	}
-	
-	public Map<Integer, Double> getRotationPeriods() {
-		return rotationPeriods;
-	}
 
 	public Timestep getTimestep() {
 		return timestep;
 	}
+
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index a2920a6b..c86ce111 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -117,6 +117,10 @@ public class GamsLocationOptimiser {
 			//if (DEBUG) LogWriter.println("     " + locId);
 			locationSet.addRecord(locId.toString());
 		}
+		GAMSSet yearSet = inDB.addSet("year", 1);
+		for (Integer year = 0; year <= 160; year++) {
+			yearSet.addRecord(year.toString());
+		}
 
 		if (DEBUG) LogWriter.println("\nPrevious crop and land areas");	
 		GAMSParameter prevCropP = inDB.addParameter("previousCropArea", 2);
@@ -191,11 +195,13 @@ public class GamsLocationOptimiser {
 			
 			// Previous land covers
 			for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
-				Vector<String> v = new Vector<String>();
-				v.add(lc.getName());
-				v.add(locString);
-				double area = landUseItem.getLandCoverArea(lc, LandProtectionType.CONVERTIBLE);
-				setGamsParamValueTruncate(prevLandCoverP.addRecord(v), area, 6);
+					Vector<String> v = new Vector<String>();
+					v.add(lc.getName());
+					v.add(locString);
+					double area = landUseItem.getLandCoverArea(lc, LandProtectionType.CONVERTIBLE);
+					if (area == 0)
+						continue;
+					setGamsParamValueTruncate(prevLandCoverP.addRecord(v), area, 6);			
 			}				
 		}
 		
@@ -356,9 +362,6 @@ public class GamsLocationOptimiser {
 		addScalar(inDB, "maxIrrigChange", ModelConfig.MAX_IRRIG_CHANGE, 5);
 		
 		// Forestry
-		Map<WoodType, Double> woodDemandMap = countryInput.getWoodDemand();
-		double totalWoodDemand = woodDemandMap.values().stream().reduce(0.0, Double::sum);	
-		setGamsParamValue(demandP.addRecord("wood"), totalWoodDemand, 5);
 		
 		// Not simulating production of different wood products so need to aggregate
 		double woodMinNetImport = 0;
@@ -380,38 +383,52 @@ public class GamsLocationOptimiser {
 		// Yield from timber forest rotation
 		Map<Integer, ? extends WoodYieldData> woodYieldData = inputData.getWoodYields();
 		
-		GAMSParameter woodYieldRotaP = inDB.addParameter("woodYieldRota", 1);
-		GAMSParameter woodYieldLuc = inDB.addParameter("woodYieldLUC", 3);
+		GAMSParameter woodYieldLucP = inDB.addParameter("woodYieldLUC", 3);
+		GAMSParameter woodRespMaxYieldP = inDB.addParameter("woodYieldMax", 1);
+		GAMSParameter woodRespSlopeP = inDB.addParameter("woodYieldParam", 1);
+		
+		Map<Integer, ? extends LandUseItem> prevLandUses = inputData.getPreviousLandUse();
 		
 		for (Entry<Integer, ? extends WoodYieldData> entry : woodYieldData.entrySet()) {
 			Integer locationId = entry.getKey();
 			String locString = Integer.toString(locationId);
 			WoodYieldData wYield = entry.getValue();
 
-			Vector<String> v = new Vector<String>();
-			v.add(locString);
-			setGamsParamValue(woodYieldRotaP.addRecord(v), wYield.getYieldRota(), 5);
-			
 			for (LandCoverType fromLc : LandCoverType.getForestedTypes()) {
-				for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
-					if (fromLc.equals(toLc))
-						continue;
-					Vector<String> w = new Vector<String>();
+				for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {	
+				    Vector<String> w = new Vector<String>();
 					w.add(fromLc.getName());
 					w.add(toLc.getName());
 					w.add(locString);
-					setGamsParamValue(woodYieldLuc.addRecord(w), wYield.getYieldLuc(fromLc, toLc), 5);
+					double stock = wYield.getYieldLuc(fromLc, toLc);
+					double area = prevLandUses.get(locationId).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
+					double yield = area == 0 ? 0.0: stock / area;
+					setGamsParamValue(woodYieldLucP.addRecord(w), yield, 5);
 				}
 			}
+
+			{
+				Vector<String> w = new Vector<String>();
+				w.add(locString);
+				setGamsParamValue(woodRespMaxYieldP.addRecord(w), wYield.getMaxYieldParam(), 5);
+				setGamsParamValue(woodRespSlopeP.addRecord(w), wYield.getSlopeParam(), 5);
+			}			
+			
 		}
 		
-		GAMSParameter forestRotationP = inDB.addParameter("forestRotationPeriod", 1);
-		for (Map.Entry<Integer, Double> entry : inputData.getRotationPeriods().entrySet()) {
-			int locationId = entry.getKey();
-			double rotPeriod = entry.getValue();
-			String locString = Integer.toString(locationId);
-			setGamsParamValue(forestRotationP.addRecord(locString), rotPeriod, 0);
+		Map<Integer, Map<WoodType, Double>> woodDemandMap = countryInput.getWoodDemand();
+		GAMSParameter woodDemandP = inDB.addParameter("expectedWoodSupply", 1);
+		double woodNetImportsExpected = countryInput.getPreviousWoodUsageData().values().stream().mapToDouble(o -> o.getNetImport()).sum();
+		double woodDemandExpected = countryInput.getWoodDemand().get(0).values().stream().reduce(0.0, Double::sum);
+		double woodNetImportsPropOfDemand = woodDemandExpected > 0 ? woodNetImportsExpected / woodDemandExpected : 0;
+		for (Integer i : woodDemandMap.keySet()) {
+			double demand = woodDemandMap.get(i).values().stream().reduce(0.0, Double::sum);
+			double productionExpected = Math.max(0.0, demand - demand * woodNetImportsPropOfDemand);			
+			setGamsParamValue(woodDemandP.addRecord(i.toString()), productionExpected, 5);
 		}
+		
+		double currentWoodDemand = woodDemandMap.get(0).values().stream().reduce(0.0, Double::sum);
+		setGamsParamValue(demandP.addRecord("wood"), currentWoodDemand, -1);
 
 		// Carbon 
 		setGamsParamValue(demandP.addRecord("carbonCredits"), countryInput.getCarbonDemand(), 5);
@@ -592,9 +609,9 @@ public class GamsLocationOptimiser {
 		Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges = new HashMap<Integer, ArrayList<LandCoverChangeItem>>();
 		
 		for (GAMSVariableRecord rec : varLandCoverChange) {
-			String fromLcStr = rec.getKeys()[0];
-			String toLcStr = rec.getKeys()[1];
-			String locationName = rec.getKeys()[2];
+			String fromLcStr = rec.getKey(0);
+			String toLcStr = rec.getKey(1);
+			String locationName = rec.getKey(2);
 			int locId = Integer.parseInt(locationName);
 			double change = rec.getLevel();
 			
@@ -608,14 +625,9 @@ public class GamsLocationOptimiser {
 			ArrayList<LandCoverChangeItem> changesList = landCoverChanges.computeIfAbsent(locId, k -> new ArrayList<LandCoverChangeItem>());
 			changesList.add(new LandCoverChangeItem(fromLc, toLc, change));
 		}
-
 		
 		// Timber harvest
-		// Ignore harvest from LUC during calibration as we want all wood demand to be met by rotation harvest
-		//double woodSupplyLuc = (ModelConfig.IS_CALIBRATION_RUN) ? 0.0 : outDB.getVariable("woodSupplyLuc").getFirstRecord().getLevel();
-		double woodSupplyLuc = 0.0;
-		double woodSupplyRota = outDB.getVariable("woodSupplyRota").getFirstRecord().getLevel();
-		double totalWoodHarvest = woodSupplyLuc + woodSupplyRota;
+		double woodSupplyRota = outDB.getParameter("woodSupplyRota").getFirstRecord().getValue();
 		double netWoodImport = getParmValue(parmNetImports, "wood");
 		
 		// Need to disaggregate wood harvest and import by type. This is just for future proofing. Eventually, could have separate harvest for each wood type.
@@ -631,17 +643,25 @@ public class GamsLocationOptimiser {
 			double previousHarvest = previousWoodUsageMap.get(wt).getHarvest();
 			
 			// assuming equal split if no previous harvest
-			double newHarvest = (previousWoodHarvestSum > 0) ? totalWoodHarvest * (previousHarvest / previousWoodHarvestSum) : totalWoodHarvest / numWoodTypes;
-			newWoodUsageMap.put(wt, new WoodUsageData(newHarvest, newNetImport, woodSupplyLuc / numWoodTypes));
+			double newHarvest = (previousWoodHarvestSum > 0) ? woodSupplyRota * (previousHarvest / previousWoodHarvestSum) : woodSupplyRota / numWoodTypes;
+			newWoodUsageMap.put(wt, new WoodUsageData(newHarvest, newNetImport, 0));
+		}
+		
+		Map<Integer, Double> rotationPeriods = new HashMap<Integer, Double>();
+		GAMSVariable varRotaP = outDB.getVariable("rotationPeriod");
+		for (GAMSVariableRecord rec : varRotaP) {
+			int loc = Integer.parseInt(rec.getKey(0));
+			double period = rec.getLevel();
+			rotationPeriods.put(loc, period);
 		}
 		
 		// Carbon
 		double carbonSequestered = outDB.getParameter("netCarbonCredits").getFirstRecord().getValue();
 		double netCarbonImport = getParmValue(parmNetImports, "carbonCredits");
-		double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue();
-		CarbonUsageData carbonUsageData = new CarbonUsageData(carbonSequestered, netCarbonImport, netCarbonFlux);
+		CarbonUsageData carbonUsageData = new CarbonUsageData(carbonSequestered, netCarbonImport, 0.0);
 
-		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, carbonUsageData, newWoodUsageMap);
+		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, carbonUsageData, 
+				newWoodUsageMap, rotationPeriods);
 		return results;
 	}
 
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index a5b82623..309b8eee 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -21,12 +21,14 @@ public class GamsLocationOutput {
 	Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges;
 	private CarbonUsageData carbonUsageData;
 	private Map<WoodType, WoodUsageData> woodUsageData;
+	private Map<Integer, Double> rotationPeriods;
 	
 	public GamsLocationOutput(ModelStat status, 
 			Map<Integer, LandUseItem> landUses, 
 			Map<CropType, CropUsageData> cropUsageData,
 			Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges,
-			CarbonUsageData carbonUsageData, Map<WoodType, WoodUsageData> woodUsageData) {
+			CarbonUsageData carbonUsageData, Map<WoodType, WoodUsageData> woodUsageData,
+			Map<Integer, Double> rotationPeriods) {
 		super();
 		this.status = status;
 		this.landUses = landUses;
@@ -34,6 +36,7 @@ public class GamsLocationOutput {
 		this.landCoverChanges = landCoverChanges;
 		this.carbonUsageData = carbonUsageData;
 		this.woodUsageData = woodUsageData;
+		this.rotationPeriods = rotationPeriods;
 	}
 	
 	public ModelStat getStatus() {
@@ -58,4 +61,9 @@ public class GamsLocationOutput {
 	public Map<WoodType, WoodUsageData> getWoodUsageData() {
 		return woodUsageData;
 	}
+
+	public Map<Integer, Double> getRotationPeriods() {
+		return rotationPeriods;
+	}
+	
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 78d0f007..492cd248 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -1,18 +1,23 @@
 package ac.ed.lurg.country.gams;
 
 import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.HashMap;
 import java.util.HashSet;
+import java.util.List;
 import java.util.Map;
 import java.util.Map.Entry;
 import java.util.Set;
+
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.carbon.CarbonFluxItem;
 import ac.ed.lurg.country.LandCoverChangeItem;
 import ac.ed.lurg.forestry.WoodYieldData;
 import ac.ed.lurg.forestry.WoodYieldItem;
+import ac.ed.lurg.forestry.WoodYieldRasterSet;
 import ac.ed.lurg.landuse.Intensity;
 import ac.ed.lurg.landuse.IrrigationItem;
+import ac.ed.lurg.landuse.LandCoverTile;
 import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.CropType;
@@ -58,6 +63,26 @@ public class GamsRasterOptimiser {
 
 	private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {		
 		RasterSet<LandUseItem> newIntensityRaster = allocAreas(gamsInput.getPreviousLandUse(), gamsOutput, gamsInput.getTimestep().getYear());
+		Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges =gamsOutput.getLandCoverChanges(); 
+		Map<Integer, ? extends WoodYieldData> woodYieldDataMap = gamsInput.getWoodYields();
+		double woodSupplyLuc = 0;
+		for (int loc : landCoverChanges.keySet()) {
+			ArrayList<LandCoverChangeItem> lccArr = landCoverChanges.get(loc);
+			WoodYieldData wyData = woodYieldDataMap.get(loc);
+			for (LandCoverChangeItem lccItem : lccArr) {
+				double area = lccItem.getArea();
+				double yield = wyData.getYieldLuc(lccItem.getFromLandCover(), lccItem.getToLandCover());
+				woodSupplyLuc += area * yield;
+			}
+		}
+
+		if (!ModelConfig.IS_CALIBRATION_RUN) {
+			double woodNetImports = gamsOutput.getWoodUsageData().values().stream().mapToDouble(o -> o.getNetImport()).sum();
+			double woodDemand = gamsInput.getCountryInput().getWoodDemand().get(0).values().stream().reduce(0.0, Double::sum);
+			double residual = allocateTimberHarvest(rasterInputData.getPreviousLandUses(), gamsOutput.getRotationPeriods(), woodDemand, woodNetImports, woodSupplyLuc);
+			if (residual > 0.0)
+				LogWriter.printlnWarning("Not enough wood to harvest");
+		}
 
 		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(),
 				gamsOutput.getCarbonUsageData(), gamsOutput.getWoodUsageData());
@@ -127,7 +152,7 @@ public class GamsRasterOptimiser {
 					LandCoverType fromLC = lccItem.getFromLandCover();
 					LandCoverType toLC = lccItem.getToLandCover();
 					double change = lccItem.getArea();
-					allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId);
+					allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId, false);
 					
 				}
 			}
@@ -148,7 +173,8 @@ public class GamsRasterOptimiser {
 		return newLandUseRaster;
 	}
 	
-	private void allocAllLandCropsTop(RasterSet<LandUseItem> newLandUseRaster, Set<RasterKey> keys, LandCoverType toLC, LandCoverType fromLC, double change, Integer locId) {
+	private void allocAllLandCropsTop(RasterSet<LandUseItem> newLandUseRaster, Set<RasterKey> keys, LandCoverType toLC, 
+			LandCoverType fromLC, double change, Integer locId, boolean natToTimb) {
 		LandCoverType toType = toLC;
 		LandCoverType fromType = fromLC;
 
@@ -159,14 +185,15 @@ public class GamsRasterOptimiser {
 			fromType = toLC;
 		}
 
-		double shortfall = allocAllLandCrops(newLandUseRaster, keys, toType, fromType, change);
+		double shortfall = allocAllLandCrops(newLandUseRaster, keys, toType, fromType, change, natToTimb);
 
 		if (shortfall > 0.00001) {
 			LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes: from " + fromLC + " to " + toLC + " " + locId + ": " + shortfall);
 		}
 	}
 
-	private double allocAllLandCrops(RasterSet<LandUseItem> newLandUseRaster, Set<RasterKey> keys, LandCoverType toType, LandCoverType fromType, double change) {
+	private double allocAllLandCrops(RasterSet<LandUseItem> newLandUseRaster, Set<RasterKey> keys, LandCoverType toType, 
+			LandCoverType fromType, double change, boolean natToTimb) {
 		if (DEBUG) LogWriter.println("allocAllLandCrops: from " + fromType + " to " + toType + " change " + change, 3);
 		if (change < 0.00001)
 			return 0;
@@ -186,7 +213,12 @@ public class GamsRasterOptimiser {
 
 			if (newLandUseItem!=null) {
 				double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getLandCoverArea(fromType, LandProtectionType.CONVERTIBLE)/totalUnprotectedLC : change / keys.size();
-				double shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
+				double shortfall = 0;
+				if (natToTimb) {
+					shortfall = newLandUseItem.convertNaturalToTimber(cellChange);
+				} else {
+					shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
+				}
 				if (shortfall == 0)
 					keysWithSpace.add(key);
 				else
@@ -200,6 +232,65 @@ public class GamsRasterOptimiser {
 		return totalShortfall;
 	}
 
+	private double allocateTimberHarvest(RasterSet<LandUseItem> landUseRaster, Map<Integer, Double> rotationPeriods,
+			double demand, double net_imports, double woodSupplyLuc) {
+
+		double stockToRemove = demand - net_imports - woodSupplyLuc;
+		if (stockToRemove <= 0)
+			return 0.0;
+		
+		Map<LandCoverType, Map<Integer, Double>> stock = new HashMap<LandCoverType, Map<Integer, Double>>();
+		for (RasterKey key : landUseRaster.keySet()) {
+			LandUseItem luItem = landUseRaster.get(key);
+			WoodYieldItem wyItem = rasterInputData.getWoodYields().get(key);
+			if (luItem == null | wyItem == null) 
+				continue;
+			IntegerRasterItem locId = mapping.get(key);
+			double rotaPeriod = locId == null ? 160 : rotationPeriods.get(locId.getInt());
+			for (LandCoverType lcType : LandCoverType.getWoodProducers()) {
+				Map<Integer, Double> mapping = stock.computeIfAbsent(lcType, k -> new HashMap<Integer, Double>());
+				LandCoverTile tiles = luItem.getLandCoverTiles().get(lcType);
+				for (int age : tiles.getAgeKeys()) {
+					if (age >= rotaPeriod) {
+						double s = tiles.getArea(LandProtectionType.CONVERTIBLE, age) * wyItem.getYieldRota(lcType, age);
+						mapping.merge(age, s, (a, b) -> a + b);
+					}
+				}
+			}
+		}
+
+		List<LandCoverType> lcOrder = new ArrayList<LandCoverType>(Arrays.asList(LandCoverType.TIMBER_FOREST, LandCoverType.NATURAL));
+		for (LandCoverType lcType : lcOrder) {
+			Map<Integer, Double> stockMapping = stock.get(lcType);
+			for (int age = 251; age >= 0; age--) {
+				if (stockMapping.containsKey(age)) {
+					double s = stockMapping.get(age);
+					double toRemoveThisTime = Math.min(stockToRemove, s);
+					double factor = 1 - toRemoveThisTime / s;
+					for (RasterKey key : landUseRaster.keySet()) {
+						LandUseItem luItem = landUseRaster.get(key);
+						if (luItem == null)
+							continue;
+						LandCoverTile tiles = luItem.getLandCoverTiles().get(lcType);
+						IntegerRasterItem locId = mapping.get(key);
+						double rotaPeriod = locId == null ? 160 : rotationPeriods.get(locId.getInt());
+						if (age >= rotaPeriod) {
+							double previousArea = tiles.getArea(LandProtectionType.CONVERTIBLE, age);
+							if (previousArea == 0.0) 
+								continue;
+							double newArea = previousArea * factor;
+							tiles.setArea(LandProtectionType.CONVERTIBLE, age, newArea);
+						}
+					}
+					stockToRemove -= toRemoveThisTime;
+					if (stockToRemove <= 1e-10) 
+						break;
+				}
+			}
+		}
+		return stockToRemove;
+	}
+
 	@SuppressWarnings("serial")
 	private GamsLocationInput convertFromRaster(GamsRasterInput rasterInputData) {
 		// as a first attempt only going to look at pasture and cereal yields, assuming other yields will be correlated to one or the other.
@@ -239,8 +330,7 @@ public class GamsRasterOptimiser {
 		};
 		LazyTreeMap<Integer, WoodYieldData> aggregatedWoodYields = new LazyTreeMap<Integer, WoodYieldData>() { 
 			protected WoodYieldData createValue() { return new WoodYieldData(); }
-		};		
-		Map<Integer, Double> aggregatedForestRotation = new HashMap<Integer, Double>();
+		};
 
 		int yRespFound = 0, yRespMissing = 0;
 		int mappingsFound = 0, mappingsMissing = 0;
@@ -312,29 +402,32 @@ public class GamsRasterOptimiser {
 			}			
 			
 			if (woodYieldItem != null) {
-				for (Map.Entry<LccKey, Double> wyEntry : woodYieldItem.getYieldMap().entrySet()) {
+				// wood from LUC
+				for (Map.Entry<LccKey, Map<Integer, Double>> wyEntry : woodYieldItem.getYieldLucMap().entrySet()) {
 					LccKey lccKey = wyEntry.getKey();
-					double wYieldThisTime = wyEntry.getValue();
+					Map<Integer, Double> yieldMap = wyEntry.getValue();
+					double stock = 0;
+					for (int age : yieldMap.keySet()) {
+						double y = yieldMap.get(age);
+						double a = landUseItem.getLandCoverArea(lccKey.getFromLc(), LandProtectionType.CONVERTIBLE, age);
+						stock += y * a;
+					}
+					double wYieldThisTime = stock;
 					double wYieldSoFar = aggWYield.getYieldLuc(lccKey);
-					double areaThisTime = landUseItem.getLandCoverArea(lccKey.getFromLc());
-					double areaSoFar = aggLandUse.getLandCoverArea(lccKey.getFromLc());
-					double wYieldAgg = aggregateMean(wYieldSoFar, areaSoFar, wYieldThisTime, areaThisTime);
-					aggWYield.setYieldLuc(lccKey, wYieldAgg);		
+					double wYieldAgg = wYieldSoFar + wYieldThisTime;
+					aggWYield.setYieldLuc(lccKey, wYieldAgg);
 				}
-
-				{
-					// Rotation wood yield
-					double wYieldThisTime = woodYieldItem.getYieldAtRotation();
-					double wYieldSoFar = aggWYield.getYieldRota();
-					double wYieldAgg = aggregateMean(wYieldSoFar, suitableAreaSoFar, wYieldThisTime, suitableAreaThisTime);
-					aggWYield.setYieldRota(wYieldAgg);	
+				
+				for (Map.Entry<Integer, Double> entry : woodYieldItem.getYieldResponseMap().entrySet()) {
+					double areaThisTime = landUseItem.getLandCoverArea(LandCoverType.TIMBER_FOREST);
+					double areaSoFar = aggLandUse.getLandCoverArea(LandCoverType.TIMBER_FOREST);
+					int age = entry.getKey();
+					double yieldThisTime = entry.getValue();
+					double yieldSoFar = aggWYield.getYieldRota(age);
+					double yieldAgg = aggregateMean(yieldSoFar, areaSoFar, yieldThisTime, areaThisTime);
+					aggWYield.setYieldRota(age, yieldAgg);
 				}
 
-				// Rotation period
-				double forestRotaThisTime = woodYieldItem.getOptimalRotation();
-				double forestRotaSoFar = (aggregatedForestRotation.containsKey(clusterId)) ? aggregatedForestRotation.get(clusterId) : forestRotaThisTime;
-				double forestRotaAgg = aggregateMean(forestRotaSoFar, suitableAreaSoFar, forestRotaThisTime, suitableAreaThisTime);
-				aggregatedForestRotation.put(clusterId, forestRotaAgg);				
 			}
 
 			// Crops yields and area fractions
@@ -391,7 +484,6 @@ public class GamsRasterOptimiser {
 				double protAreaSoFar = aggLandUse.getLandCoverArea(lcType, LandProtectionType.PROTECTED);
 				aggLandUse.setLandCoverArea(lcType, LandProtectionType.PROTECTED, protAreaSoFar + protAreaThisTime);
 			}			
-
 		}	
 
 		LogWriter.println("Mapping: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + mappingsFound + ", countMissing=" + mappingsMissing, 3);
@@ -399,25 +491,9 @@ public class GamsRasterOptimiser {
 
 		checkedTotalAreas(landUseRaster, "before");
 		checkedTotalAreas(aggregatedAreas, "after");
-		
-		
-		if (!ModelConfig.IS_FORESTRY_ON) {
-			for (Integer locId : aggregatedYields.keySet()) {
-				aggregatedForestRotation.put(locId, 1.0); // shouldn't matter what this is as long as it's >0
-			}
-		}
-		
-		// TODO shouldn't have missing data in the first place
-		double meanRotation = aggregatedForestRotation.values().stream().mapToDouble(o -> o.doubleValue()).sum() / aggregatedForestRotation.keySet().size();
-		for (Integer locId : aggregatedYields.keySet()) {
-			if (!aggregatedForestRotation.containsKey(locId)) {
-				aggregatedForestRotation.put(locId, meanRotation);
-				LogWriter.printlnWarning("GamsRasterOptimister.convertFromRaster missing timber rotation; substituting average");
-			}
-		}
 
 		return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, 
-				aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput(), aggregatedForestRotation);
+				aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
 	}
 
 	private void logWarningWithCoord(String message, RasterKey key, YieldRaster yieldRaster, CropType crop) {
diff --git a/src/ac/ed/lurg/demand/AbstractDemandManager.java b/src/ac/ed/lurg/demand/AbstractDemandManager.java
index 127c3066..d72c4b85 100644
--- a/src/ac/ed/lurg/demand/AbstractDemandManager.java
+++ b/src/ac/ed/lurg/demand/AbstractDemandManager.java
@@ -114,27 +114,31 @@ public abstract class AbstractDemandManager implements Serializable {
 
 	protected abstract Map<WoodType, Double> getWoodDemand(SingleCountry country, int year);
 	
-	public Map<WoodType, Double> getWoodDemandComposite(CompositeCountry cc, int year) {
+	public Map<Integer, Map<WoodType, Double>> getWoodDemandComposite(CompositeCountry cc, int year) {
 		if (!ModelConfig.CHANGE_DEMAND_YEAR)
 			year = ModelConfig.BASE_YEAR;
-
-		Map<WoodType, Double> compositeDemandMap = new HashMap<WoodType, Double>();
-		Map<WoodType, Double> singleDemandMap;
-
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
-			singleDemandMap = getWoodDemand(c, year);
-			if (singleDemandMap.isEmpty()) {
-				LogWriter.printlnWarning("Wood demand not found for: " + c.getCountryName());
-				continue;
-			}
-			for (WoodType w : WoodType.values()) {
-				double totalDemand = (compositeDemandMap.containsKey(w)) ? compositeDemandMap.get(w) : 0.0;
-				double demand = singleDemandMap.get(w);
-				totalDemand += demand;
-				compositeDemandMap.put(w, totalDemand);				
+		Map<Integer, Map<WoodType, Double>>  annualDemandMap = new HashMap<Integer, Map<WoodType, Double>>();
+		for (int i = 0; i <= 2100 - year; i++) {
+			Map<WoodType, Double> compositeDemandMap = new HashMap<WoodType, Double>();
+			Map<WoodType, Double> singleDemandMap;
+
+			for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+				singleDemandMap = getWoodDemand(c, year + i);
+				if (singleDemandMap.isEmpty()) {
+					LogWriter.printlnWarning("Wood demand not found for: " + c.getCountryName());
+					continue;
+				}
+				for (WoodType w : WoodType.values()) {
+					double totalDemand = (compositeDemandMap.containsKey(w)) ? compositeDemandMap.get(w) : 0.0;
+					double demand = singleDemandMap.get(w);
+					totalDemand += demand;
+					compositeDemandMap.put(w, totalDemand);				
+				}
 			}
+			annualDemandMap.put(i, compositeDemandMap);
 		}
-		return compositeDemandMap;
+
+		return annualDemandMap;
 	}
 	
 	public double getWorldCarbonDemand(int year) {
diff --git a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
index a5760f6c..41760185 100644
--- a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
+++ b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
@@ -1,16 +1,7 @@
 package ac.ed.lurg.forestry;
 
-import java.io.BufferedWriter;
-import java.io.File;
-import java.io.FileWriter;
-import java.io.IOException;
-import java.util.Map.Entry;
-
 import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.output.AbstractLandUseOutputer;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.utils.LogWriter;
-import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterSet;
 
 public class ForestryDataOutputer extends AbstractLandUseOutputer {
@@ -23,10 +14,11 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
 
 	@Override
 	public void writeOutput() {
+		/*
 		File outputDir = getOutputDir(year);
 		String outputFileName = outputDir.getPath() + File.separator + "Forestry.txt";
 		BufferedWriter outputWriter = null;
-		try {
+		try {	
 			outputWriter = new BufferedWriter(new FileWriter(outputFileName, false));
 			outputWriter.write("Lon Lat RotationAge YieldAtRotation RotationCutArea StandingStockTimber StandingStockCarbon StandingStockNatural");
 			outputWriter.newLine();
@@ -40,7 +32,6 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
 				if (wyItem == null)
 					continue;
 
-				
 				double lat = landUseRaster.getXCoordin(key);
 				double lon = landUseRaster.getYCoordin(key);
 				int rotaAge = wyItem.getOptimalRotation();
@@ -53,6 +44,7 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
 				outputWriter.write(String.format("%.2f %.2f %d %.14f %.14f %.14f %.14f %.14f", lat, lon, 
 						rotaAge, rotaYield, rotaArea, stockTimber, stockCarbon, stockNatural));
 				outputWriter.newLine();
+				
 			}
 		}
 		catch (IOException e) {
@@ -69,6 +61,6 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
 			}
 		}
 		
-		
+		*/
 	}
 }
diff --git a/src/ac/ed/lurg/forestry/ForestryPlanner.java b/src/ac/ed/lurg/forestry/ForestryPlanner.java
new file mode 100644
index 00000000..7de2c87c
--- /dev/null
+++ b/src/ac/ed/lurg/forestry/ForestryPlanner.java
@@ -0,0 +1,7 @@
+package ac.ed.lurg.forestry;
+
+public class ForestryPlanner {
+	
+
+
+}
diff --git a/src/ac/ed/lurg/forestry/WoodYieldData.java b/src/ac/ed/lurg/forestry/WoodYieldData.java
index ab8b8b00..e3a504c1 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldData.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldData.java
@@ -7,37 +7,39 @@ import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
 
 public class WoodYieldData {
-	private double yieldRota; // wood yield at current rotation length
-	private Map<LccKey, Double> yieldLuc = new HashMap<LccKey, Double>(); // wood yield from clearing current forest
-	
-	public void setYieldRota(double yield) {
-		yieldRota = yield;
-	}
+	private Map<LccKey, Double> yieldLuc = new HashMap<LccKey, Double>();
+	private WoodYieldResponse yieldResponse = new WoodYieldResponse();
 	
 	public void setYieldLuc(LccKey key, double yield) {
 		yieldLuc.put(key, yield);
 	}
 	
-	public double getYieldRota() {
-		return yieldRota;
-	}
-	
 	public double getYieldLuc(LandCoverType fromLc, LandCoverType toLc) {
 		LccKey key = new LccKey(fromLc, toLc);
-		if (yieldLuc.containsKey(key)) 
-			return yieldLuc.get(key);
-		
-		return Double.NaN;
+		return getYieldLuc(key);
 	}
 	
 	public double getYieldLuc(LccKey key) {
-		if (yieldLuc.containsKey(key)) 
-			return yieldLuc.get(key);
-		
-		return Double.NaN;
+		return yieldLuc.getOrDefault(key, 0.0);
 	}
 	
 	public Map<LccKey, Double> getYieldLucMap() {
 		return yieldLuc;
 	}
+
+	public double getYieldRota(int age) {
+		return yieldResponse.getYield(age);
+	}
+	
+	public void setYieldRota(int age, double yield) {
+		yieldResponse.setYield(age, yield);
+	}
+	
+	public double getMaxYieldParam() {
+		return yieldResponse.getMaxYield(); 
+	}
+
+	public double getSlopeParam() {
+		return yieldResponse.getSlope();
+	}
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index d27c6e4a..aa0aea7c 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -1,138 +1,97 @@
 package ac.ed.lurg.forestry;
 
-import java.util.ArrayList;
-import java.util.Collections;
 import java.util.HashMap;
-import java.util.List;
 import java.util.Map;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.landuse.LandCoverTile;
-import ac.ed.lurg.landuse.LandKey;
 import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.types.LandProtectionType;
-import ac.ed.lurg.utils.Interpolator;
 import ac.sac.raster.RasterItem;
 
 public class WoodYieldItem implements RasterItem {	
-	private Map<LccKey, Double> yield = new HashMap<LccKey, Double>();
-	private int optimalRotation; // Only applies to TIMBER_FOREST
-	private double yieldAtRotation; // Only applies to TIMBER_FOREST
-	private double rotationArea = 0; // current area clear-cut
+	private Map<LccKey, Map<Integer, Double>> yieldLuc = new HashMap<LccKey, Map<Integer, Double>>();
+	private Map<LandCoverType, Map<Integer, Double>> yieldRota = new HashMap<LandCoverType, Map<Integer, Double>>();
+	private WoodYieldResponse yieldResponse = new WoodYieldResponse();
 	private double maxYield = 0; // Yield at maximum age. Used for clustering
+	private double potentialYield = 0;
 	
 	public WoodYieldItem() {}
 	
-	public void calcYieldData(LccKey key, Double[] yields, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
-		// Mean wood yield for grid cell
+	public void setYieldLuc(LccKey key, Double[] yields, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
 		LandCoverTile tiles = landUseTiles.get(key.getFromLc());
-
-		double totalArea = tiles.getTotalArea(LandProtectionType.CONVERTIBLE); // Assuming no harvest from protected areas
-		if (totalArea < 1e-6) { // Avoid floating point errors due to small areas
-			this.yield.put(key, 0.0) ;
-		} else {				
-			double totalYield = 0;
-			for (int age : tiles.getAgeKeys()) {
-				int ageCapped = Math.min(age, maxAge);
-				totalYield += yields[ageCapped] * tiles.getArea(LandProtectionType.CONVERTIBLE, age);
-			}
-			double meanYield = totalYield / totalArea;
-			this.yield.put(key, meanYield);
+		Map<Integer, Double> yieldMap = yieldLuc.computeIfAbsent(key, k -> new HashMap<Integer, Double>());
+		for (Integer age : tiles.getAgeKeys()) {
+			int ageCapped = Math.min(age, maxAge);
+			yieldMap.put(ageCapped, yields[ageCapped]);
+		}
+	}
+	
+	public void setYieldRota(LandCoverType lcType, Double[] yields, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
+		LandCoverTile tiles = landUseTiles.get(lcType);
+		Map<Integer, Double> yieldMap = yieldRota.computeIfAbsent(lcType, k -> new HashMap<Integer, Double>());
+		for (Integer age : tiles.getAgeKeys()) {
+			int ageCapped = Math.min(age, maxAge);
+			yieldMap.put(ageCapped, yields[ageCapped]);
 		}
 
 		// For clustering
-		if (key.getFromLc().equals(LandCoverType.NATURAL) && key.getToLc().equals(LandCoverType.PASTURE))
+		if (lcType.equals(LandCoverType.TIMBER_FOREST)) {
 			maxYield = yields[ModelConfig.MAX_FOREST_ROTATION_PERIOD];
+			potentialYield = yields[20] / 20;
+		}
 	}
-
 	
-	public void calcRotationData(Map<Integer, Double> woodYields, Map<LandCoverType, LandCoverTile> landUseTiles, double woodPrice) {
-		
-		// This is in case ModelConfig.ENABLE_VEGETATION_CLEARANCE_COST is true. We want conversion yields but don't care about rotation.
-		if (!ModelConfig.IS_FORESTRY_ON) {
-			optimalRotation = 1;
-			yieldAtRotation = 0;
-			return;
-		}
-		
-		// We only have yields for every CARBON_DATA_TIMESTEP_SIZE years so need to interpolate.
-		Double[] yieldsInterp = new Double[ModelConfig.MAX_FOREST_ROTATION_PERIOD + 1];
-		woodYields.putIfAbsent(0, 0.0);
-		List<Integer> ages = woodYields.keySet().stream().toList();
-		List<Integer> agesSorted = new ArrayList<Integer>(ages);
-		Collections.sort(agesSorted);
-		for (int i = 0; i < agesSorted.size() - 1; i++) {
-			int fromAge = agesSorted.get(i);
-			int toAge = agesSorted.get(i + 1);
-			double fromYield = woodYields.get(fromAge);
-			double toYield = woodYields.get(toAge);
-			for (int age = fromAge; age < toAge; age++) {
-				if (age > ModelConfig.MAX_FOREST_ROTATION_PERIOD)
-					break;
-				double factor = ((double) age - fromAge) / (toAge - fromAge);
-				double yield = Interpolator.interpolate(fromYield, toYield, factor);
-				yieldsInterp[age] = yield;
-			}
-		}
-		
-		// Optimal rotation
-		Map<Integer, Double> levArr = new HashMap<Integer, Double>();
-		for (int age = ModelConfig.MIN_FOREST_ROTATION_PERIOD; age < ModelConfig.MAX_FOREST_ROTATION_PERIOD; age++) {
-			double lev = (woodPrice * yieldsInterp[age] * Math.exp(-ModelConfig.DISCOUNT_RATE * age) - ModelConfig.FOREST_MANAGEMENT_COST) / (1 - Math.exp(-ModelConfig.DISCOUNT_RATE * age));
-			levArr.put(age, lev);
-		}
-		
-		double maxLev = Collections.max(levArr.values());
-		List<Integer> candidates = new ArrayList<Integer>(); // There may be several equal maximum values
-		for (int age : levArr.keySet()) {
-			if (levArr.get(age) == maxLev) {
-				candidates.add(age);
-			}
-		}
-		
-		optimalRotation = Collections.min(candidates); // Choose shortest rotation		
-		yieldAtRotation = yieldsInterp[optimalRotation] / optimalRotation;
-		
-		// Area that is currently being clear-cut
-		LandCoverTile forestTile = landUseTiles.get(LandCoverType.TIMBER_FOREST);
-		for (int age : forestTile.getAgeKeys()) {
-			if (age >= optimalRotation) {
-				LandKey key = new LandKey(LandProtectionType.CONVERTIBLE, age);
-				double area = forestTile.getArea(key);
-				this.rotationArea += area;
-			}
+	public void setYieldResponse(Double[] yields, int modelYear) {
+		int dataYear = Math.floorDiv(modelYear - 10, ModelConfig.CARBON_DATA_TIMESTEP_SIZE) * ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 10;
+		for (int age = modelYear - dataYear; age <= 160; age += ModelConfig.CARBON_DATA_TIMESTEP_SIZE) {
+			yieldResponse.setYield(age, yields[age]);
 		}
+		yieldResponse.setYield(0, 0); // in case missing
 	}
 	
-	public void setDefaultForMissingData() {
-		
+	public WoodYieldResponse getYieldResponse() {
+		return yieldResponse;
+	}
+
+
+	public double getYieldLuc(LandCoverType fromLc, LandCoverType toLc, int age) {
+		return getYieldLuc(new LccKey(fromLc, toLc), age);
 	}
 	
-	public double getYield(LandCoverType fromLc, LandCoverType toLc) {
-		return yield.get(new LccKey(fromLc, toLc));
+	public double getYieldLuc(LccKey key, int age) {
+		if (yieldLuc.containsKey(key)) {
+			return yieldLuc.get(key).getOrDefault(age, 0.0);
+		} else {
+			return 0.0;
+		}
 	}
 	
-	public double getYield(LccKey key) {
-		return yield.get(key);
+	public double getYieldRota(LandCoverType lcType, int age) {
+		if (yieldRota.containsKey(lcType)) {
+			return yieldRota.get(lcType).getOrDefault(age, 0.0);
+		} else {
+			return 0.0;
+		}
 	}
 	
-	public Map<LccKey, Double> getYieldMap() {
-		return yield;
+	public double getMaxYield() {
+		return maxYield;
 	}
-
-	public int getOptimalRotation() {
-		return optimalRotation;
+	
+	public Map<LccKey, Map<Integer, Double>> getYieldLucMap() {
+		return yieldLuc;
 	}
-
-	public double getYieldAtRotation() {
-		return yieldAtRotation;
+	
+	public Map<LandCoverType, Map<Integer, Double>> getYieldRotaMap() {
+		return yieldRota;
 	}
-
-	public double getMaxYield() {
-		return maxYield;
+	
+	public Map<Integer, Double> getYieldResponseMap() {
+		return yieldResponse.getYieldMap();
 	}
 	
-	public double getRotationArea() {
-		return rotationArea;
+	public double getPotentialYield() {
+		return potentialYield;
 	}
+	
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java b/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
index 635e1d79..ae4ab078 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
@@ -3,6 +3,7 @@ package ac.ed.lurg.forestry;
 import java.util.Collection;
 import java.util.Map;
 
+import ac.ed.lurg.landuse.LandCoverTile;
 import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.LandProtectionType;
@@ -29,17 +30,19 @@ public class WoodYieldRasterSet extends RasterSet<WoodYieldItem> {
 		return subsetForestGrowthRaster;
 	}
 	
-	public static double getMaxWoodHarvest(RasterSet<WoodYieldItem> woodYields, RasterSet<LandUseItem> landUses) {
+	public static double getPotentialTimberProduction(RasterSet<WoodYieldItem> woodYieldData, RasterSet<LandUseItem> landUseData) {
 		double total = 0;
-		for (Map.Entry<RasterKey, WoodYieldItem> entry : woodYields.entrySet()) {
-			RasterKey key = entry.getKey();
-			WoodYieldItem wyItem = entry.getValue();
-			LandUseItem luItem = landUses.get(key);
-			if (wyItem == null)
+		for (RasterKey key : landUseData.keySet()) {
+			LandUseItem luItem = landUseData.get(key);
+			WoodYieldItem wyItem  = woodYieldData.get(key);
+			if (luItem == null | wyItem == null)
 				continue;
-			double area = luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE);
-			double rotaYield = wyItem.getYieldAtRotation();
-			total += area * rotaYield;
+			double area = 0;
+			for (LandCoverType lcType : LandCoverType.getForestedTypes()) {
+				area += luItem.getLandCoverArea(lcType, LandProtectionType.CONVERTIBLE);
+			}
+			double yield = wyItem.getPotentialYield(); // assuming rotation period of 20 years
+			total += area * yield;
 		}
 		return total;
 	}
diff --git a/src/ac/ed/lurg/forestry/WoodYieldReader.java b/src/ac/ed/lurg/forestry/WoodYieldReader.java
index af8ff6fd..6a9a40e1 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldReader.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java
@@ -18,6 +18,8 @@ import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.yield.YieldResponsesItem;
+import ac.sac.raster.AbstractTabularRasterReader;
 import ac.sac.raster.RasterHeaderDetails;
 import ac.sac.raster.RasterKey;
 import ac.sac.raster.RasterSet;
@@ -27,6 +29,7 @@ public class WoodYieldReader {
 	private RasterHeaderDetails rasterProj;
 	private int modelYear;
 	private Map<Integer, RasterKey> keyMap;
+	double woodPrice;
 	
 	public WoodYieldReader(RasterHeaderDetails rasterProj) {
 		this.rasterProj = rasterProj;
@@ -36,9 +39,14 @@ public class WoodYieldReader {
 		WoodYieldRasterSet woodYieldData = new WoodYieldRasterSet(rasterProj);
 		keyMap = getRasterKeys(woodYieldData);
 		modelYear = ModelConfig.IS_CALIBRATION_RUN ? ModelConfig.BASE_YEAR : timestep.getYear();
+		this.woodPrice = woodPrice;
 		long startTime = System.currentTimeMillis();
 			// Wood yield from rotation
-			readRotationData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_FILENAME, woodPrice); 
+			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_FORS_FILENAME, 
+					Arrays.asList(new LccKey(LandCoverType.TIMBER_FOREST, LandCoverType.TIMBER_FOREST))); 
+			
+			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_NTRL_FILENAME, 
+					Arrays.asList(new LccKey(LandCoverType.NATURAL, LandCoverType.NATURAL))); 
 			
 			// Wood yield from LCC
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_PAST_FILENAME, 
@@ -102,47 +110,30 @@ public class WoodYieldReader {
 				WoodYieldItem wyItem = woodYieldData.get(key.getCol(), key.getRow());
 				LandUseItem luItem = landUseRaster.get(key);
 				for (LccKey lccKey : lccMappings) {
-					wyItem.calcYieldData(lccKey, map.get(key), luItem.getLandCoverTiles(), maxAge);
+					if (lccKey.getFromLc().equals(lccKey.getToLc())) {
+						wyItem.setYieldRota(lccKey.getFromLc(), map.get(key), luItem.getLandCoverTiles(), maxAge);
+						if (lccKey.getFromLc().equals(LandCoverType.TIMBER_FOREST)) {
+							wyItem.setYieldResponse(map.get(key), modelYear);
+						}
+					} else {
+						wyItem.setYieldLuc(lccKey, map.get(key), luItem.getLandCoverTiles(), maxAge);
+					}
 				}
 			}
 		}
 		LogWriter.println("Reading " + fileName + " took " + (System.currentTimeMillis() - startTime) + " ms");
 	}
 	
-	public void readRotationData(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, String fileName, double carbonPrice) {
-		long startTime = System.currentTimeMillis();
-		int maxYear = Math.floorDiv(modelYear - 10, ModelConfig.CARBON_DATA_TIMESTEP_SIZE) * ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 10;
-		Map<RasterKey, Map<Integer, Double>> map = new HashMap<RasterKey, Map<Integer, Double>>();
-		for (int dataYear = maxYear; dataYear >= ModelConfig.CARBON_DATA_MIN_YEAR; dataYear -= ModelConfig.CARBON_DATA_TIMESTEP_SIZE) {
-			int nCol = keyMap.size();
-			int age = modelYear - dataYear;
-			String filePath = getDataDir(fileName, dataYear);
-			try {
-				RandomAccessFile reader = new RandomAccessFile(filePath, "r");
-				int startPos = nCol * age * Double.BYTES;
-				int nRead = nCol * Double.BYTES;
-				reader.seek(startPos);
-				byte[] byteData = new byte[nRead];
-				reader.read(byteData);
-				Double[] data = byteArrayToDouble(byteData);
-				for (int i = 0; i < nCol; i++) {
-					Map<Integer, Double> locData = map.computeIfAbsent(keyMap.get(i), o -> new HashMap<Integer, Double>());
-					locData.put(age, data[i] * CONVERSION_FACTOR);
-				}
-				reader.close();
-			} catch (IOException e) {
-				e.printStackTrace();
-			}
-		}
-		for (RasterKey key : map.keySet()) {
-			if (landUseRaster.containsKey(key)) {
-				WoodYieldItem wyItem = woodYieldData.get(key.getCol(), key.getRow());
-				LandUseItem luItem = landUseRaster.get(key);
-				wyItem.calcRotationData(map.get(key), luItem.getLandCoverTiles(), carbonPrice);
-			}
-		}
-		LogWriter.println("Reading " + fileName + " took " + (System.currentTimeMillis() - startTime) + " ms");
-	}
+//	public void readYieldResponses(WoodYieldRasterSet woodYieldRaster) {
+//		AbstractTabularRasterReader<WoodYieldItem> reader = new AbstractTabularRasterReader<WoodYieldItem>("\\s+", 4, woodYieldRaster) {
+//			protected void setData(RasterKey key, WoodYieldItem item, Map<String, Double> rowValues) {
+//				item.setYieldResponse(getValueForCol(rowValues, "a") * CONVERSION_FACTOR, getValueForCol(rowValues, "b"));
+//			}
+//		};
+//		int dataYear = Math.floorDiv(modelYear - 10, ModelConfig.CARBON_DATA_TIMESTEP_SIZE) * ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 10;
+//		String fileName = getDataDir(ModelConfig.WOOD_YIELD_PARAM_FILENAME, dataYear);
+//		reader.getRasterDataFromFile(fileName);
+//	}
 	
 	public Map<Integer, RasterKey> getRasterKeys(WoodYieldRasterSet woodYieldData) {
 		Map<Integer, RasterKey> keyMap = new HashMap<Integer, RasterKey>();
@@ -192,5 +183,5 @@ public class WoodYieldReader {
 	private String getDataDir(String filename, int dataYear) {
 		return ModelConfig.FORESTRY_DIR + File.separator + dataYear + File.separator + filename;
 	}
-	
+
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldResponse.java b/src/ac/ed/lurg/forestry/WoodYieldResponse.java
new file mode 100644
index 00000000..e68902f0
--- /dev/null
+++ b/src/ac/ed/lurg/forestry/WoodYieldResponse.java
@@ -0,0 +1,95 @@
+package ac.ed.lurg.forestry;
+
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.Vector;
+
+public class WoodYieldResponse {
+	private Map<Integer, Double> yields = new HashMap<Integer, Double>();
+	private Double maxYield;
+	private Double slope;
+
+	public void setYield(int age, double yield) {
+		yields.put(age, yield);
+	}
+	
+	public double getYield(int age) {
+		return yields.getOrDefault(age, Double.NaN);
+	}
+	
+	public void calcParams() {
+		if (yields.size() == 0) {
+			maxYield = 0.0;
+			slope = 0.0;
+			return;
+		}
+		
+		maxYield = yields.values().stream().reduce(0.0, Double::max) + 1e-6;
+		Vector<Double> kEst = new Vector<Double>();
+		for (int a : yields.keySet()) {
+			if (a == 0)
+				continue;
+			double y = yields.get(a);
+			kEst.add(Math.log(1 - y / maxYield) / (-a));
+		}
+		Collections.sort(kEst);
+		slope = kEst.get((int) kEst.size() / 2); 
+
+		double minA = maxYield * 0.8;
+		double stepA = (maxYield * 1.2 - minA) / 100;
+		double minB = slope * 0.8;
+		double stepB = (slope * 1.2 - minB) / 100;
+		
+		double bestA = maxYield;
+		double bestB = slope;
+		double bestError = calcResidualSqr(bestA, bestB);
+		
+		for (int i = 0; i < 100; i++) {
+			double a = minA + stepA * i;
+			for (int j = 0; j < 100; j++) {
+				double b = minB + stepB * j;
+				double error = calcResidualSqr(a, b);
+				if (error < bestError) {
+					bestA = a;
+					bestB = b;
+					bestError = error;
+				}
+			}
+		}
+		
+		maxYield = bestA;
+		slope = bestB;	
+	}
+	
+	public double getMaxYield() {
+		if (maxYield == null)
+			calcParams();
+		return maxYield;
+	}
+	
+	public double getSlope() {
+		if (slope == null)
+			calcParams();
+		return slope;
+	}
+	
+	private double yieldFunction(double a, double b, int t) {
+		return a * (1 - Math.exp(-b * t));
+	}
+	
+	private double calcResidualSqr(double maxYield, double slope) {
+		double sumOfSquares = 0;
+		for (int age : yields.keySet()) {
+			double y = yields.get(age);
+			double fit = yieldFunction(maxYield, slope, age);
+			double res2 = Math.pow(y - fit, 2);
+			sumOfSquares += res2;
+		}
+		return sumOfSquares;
+	}
+	
+	public Map<Integer, Double> getYieldMap() {
+		return yields;
+	}
+}
diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java
index e2a62e53..cd292e71 100644
--- a/src/ac/ed/lurg/landuse/LandCoverItem.java
+++ b/src/ac/ed/lurg/landuse/LandCoverItem.java
@@ -81,6 +81,7 @@ public class LandCoverItem implements RasterItem {
 	public Map<LandCoverType, LandCoverTile> getLandCoverTiles() {
 		createLandCoverTiles();
 		setUnavailableArea();
+		cleanUp();
 		return landCoverAreas;
 	}
 	
@@ -120,6 +121,10 @@ public class LandCoverItem implements RasterItem {
 			}
 
 			Map<Integer, Double> ageMap = getLandUseAgeDist().get(lcType);
+			if (ageMap == null) {
+				ageMap = new HashMap<Integer, Double>();
+				ageMap.put(15, 1.0);
+			}
 			double totalProp = ageMap.values().stream().reduce(0.0,  Double::sum);
 			if (totalProp == 0.0) { // If no distribution given then assume land is in oldest age group
 				int maxAgeGroup = ageMap.keySet().stream().max(Integer::compare).get();
@@ -132,17 +137,23 @@ public class LandCoverItem implements RasterItem {
 				}
 			}
 			
-			// Assuming land cover area all in middle of age group
 			for (Map.Entry<Integer, Double> ageEntry : ageMap.entrySet()) {
 				double prop = ageEntry.getValue();
 				if (prop == 0.0)
 					continue;
 				int ageGroup = ageEntry.getKey();
 				int minAge = (int) (ageGroup - 1) * ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE + 1;
-				int midAge = minAge + ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE / 2;
-				
-				landCoverAreas.get(lcType).setArea(LandProtectionType.CONVERTIBLE, midAge, unprotectedArea * prop);
-				landCoverAreas.get(lcType).setArea(LandProtectionType.PROTECTED, midAge, protectedArea * prop);				
+				int maxAge = minAge + ModelConfig.LAND_COVER_INIT_AGE_GROUP_SIZE - 1;
+				double groupProp = 0.2;
+				// Assuming even distribution within each age group
+				for (int a = minAge; a <= maxAge; a+=2) {
+					double unprotectedA = unprotectedArea * prop * groupProp;
+					if (unprotectedA > 0)
+						landCoverAreas.get(lcType).setArea(LandProtectionType.CONVERTIBLE, a, unprotectedA);
+					double protectedA = protectedArea * prop * groupProp;
+					if (protectedA > 0)
+						landCoverAreas.get(lcType).setArea(LandProtectionType.PROTECTED, a, protectedA);		
+				}	
 			}
 		}
 		
@@ -197,4 +208,10 @@ public class LandCoverItem implements RasterItem {
 		}
 	}
 	
+	private void cleanUp() {
+		for (LandCoverTile tile : landCoverAreas.values()) {
+			tile.cleanUp();
+		}
+	}
+	
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index fd3b53c9..4ba22ac7 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -74,6 +74,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return landCoverAreas.get(lcType).getTotalArea(lpType);
 	}
 	
+	public double getLandCoverArea(LandCoverType lcType, LandProtectionType lpType, int age) {
+		return landCoverAreas.get(lcType).getArea(lpType, age);
+	}
+	
 	public double getLandCoverArea(LandCoverType lcType) {
 		Double d = landCoverAreas.get(lcType).getTotalArea();
 		return d == null ? 0.0 : d;
@@ -139,6 +143,27 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
         return area - areaMoved;
 	}
 	
+	public double convertNaturalToTimber(double area) {
+		if (area == 0)
+			return 0;
+        double prevFrom = getLandCoverArea(LandCoverType.NATURAL, LandProtectionType.CONVERTIBLE);
+        double areaMoved;
+        areaMoved = Math.min(area, prevFrom);
+        if (Math.abs(area - areaMoved) / area > 1e-6) {
+        	throw new RuntimeException("Move unprotected area shortfall. Should not happen due to GAMS constraint.");
+        }
+		LandCoverTile toTiles = landCoverAreas.get(LandCoverType.TIMBER_FOREST);
+		LandCoverTile fromTiles = landCoverAreas.get(LandCoverType.NATURAL);
+		double prop = areaMoved / prevFrom;
+		for (int age : fromTiles.getAgeKeys()) {
+			double a = fromTiles.getArea(LandProtectionType.CONVERTIBLE, age);
+			double propA = a * prop;
+			fromTiles.removeArea(LandProtectionType.CONVERTIBLE, age, propA);
+			toTiles.addArea(LandProtectionType.CONVERTIBLE, age, propA);
+		}
+		return area - areaMoved;
+	}
+	
 	public void resetLccAreas() {
 		landCoverChange.clear();
 	}
@@ -156,6 +181,15 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		landCoverAreas.get(c).removeAllArea(p);
 		landCoverAreas.get(c).addArea(p, 0, landCover);
 	}
+	
+	public void setLandCoverArea(LandCoverType c, LandProtectionType p, int age, double d) {
+		if (Double.isNaN(d) || Double.isInfinite(d))
+			throw new RuntimeException("AreasItem for " + c + " is " + d);
+
+		double landCover = (d < 0.0) ? 0.0 : d;
+		
+		landCoverAreas.get(c).setArea(p, age, landCover);
+	}
 
 	public Map<LandCoverType, Double> getLandCoverAreaMap() {
 		Map<LandCoverType, Double> areaMap = new HashMap<LandCoverType, Double>();
diff --git a/src/ac/ed/lurg/types/LandCoverType.java b/src/ac/ed/lurg/types/LandCoverType.java
index dad53765..45318c61 100644
--- a/src/ac/ed/lurg/types/LandCoverType.java
+++ b/src/ac/ed/lurg/types/LandCoverType.java
@@ -9,26 +9,29 @@ import ac.ed.lurg.utils.LogWriter;
 
 public enum LandCoverType {
 	
-	TIMBER_FOREST("timberForest", false, true, true, true),
-	CARBON_FOREST("carbonForest", false, true, true, true),
-	NATURAL ("natural", true, true, true, false),
-	CROPLAND ("cropland", false, true, false, false),
-	PASTURE ("pasture", false, true, false, false),
-	BARREN ("barren", false, false, false, false),
-	URBAN("urban", true, false, false, false);
+	TIMBER_FOREST("timberForest", false, true, true, true, true),
+	CARBON_FOREST("carbonForest", false, true, true, true, false),
+	NATURAL ("natural", true, true, true, false, true),
+	CROPLAND ("cropland", false, true, false, false, false),
+	PASTURE ("pasture", false, true, false, false, false),
+	BARREN ("barren", false, false, false, false, false),
+	URBAN("urban", true, false, false, false, false);
 
 	private String name;
 	private boolean isProtectable;
 	private boolean isConvertible;
 	private boolean isNatural;
 	private boolean isManagedForest;
+	private boolean isWoodProducer;
 	
-	LandCoverType(String name, boolean isProtectable, boolean isConvertible, boolean isNatural, boolean isManagedForest) {
+	LandCoverType(String name, boolean isProtectable, boolean isConvertible, boolean isNatural, 
+			boolean isManagedForest, boolean isWoodProducer) {
 		this.name = name;
 		this.isProtectable = isProtectable;
 		this.isConvertible = isConvertible;
 		this.isNatural = isNatural;
 		this.isManagedForest = isManagedForest;
+		this.isWoodProducer = isWoodProducer;
 	}
 	
 	public String getName() {
@@ -94,39 +97,39 @@ public enum LandCoverType {
 	}
 	
 	public static Collection<LandCoverType> getManagedForestTypes() {
-
 		Collection<LandCoverType> managedForestTypes = new HashSet<LandCoverType>();
-
 		for (LandCoverType c : values())
 			if (c.isManagedForest)
 				managedForestTypes.add(c);
 
 		return managedForestTypes;
-
 	}
 	
 	public static Collection<LandCoverType> getForestedTypes() {
-
 		Collection<LandCoverType> naturalTypes = new HashSet<LandCoverType>();
-
 		for (LandCoverType c : values())
 			if (c.isNatural)
 				naturalTypes.add(c);
-
+		
 		return naturalTypes;
-
 	}
 	
 	public static Collection<LandCoverType> getAgriculturalTypes() {
-
 		Collection<LandCoverType> agriTypes = new HashSet<LandCoverType>();
-
 		for (LandCoverType c : values())
 			if (!c.isNatural && c.isConvertible)
 				agriTypes.add(c);
 
 		return agriTypes;
+	}
+	
+	public static Collection<LandCoverType> getWoodProducers() {
+		Collection<LandCoverType> woodProducers = new HashSet<LandCoverType>();
+		for (LandCoverType c : values())
+			if (c.isWoodProducer)
+				woodProducers.add(c);
 
+		return woodProducers;
 	}
 }
 
-- 
GitLab