diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index d92b5acb12e1bd4f80ae11d001edfa35294c766a..caecc0dd63659880278a2e5e08f9a44627ed31d6 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -12,6 +12,8 @@
  SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
  SET forest(land_cover) / timberForest, carbonForest, natural /;
  SET non_forest(land_cover) / cropland, pasture /;
+ SET non_timber_forest_before(land_cover) / cropland, pasture, carbonForest, natural /;
+ ALIAS(non_timber_forest_before, non_timber_forest_after);
  ALIAS (land_cover, land_cover_before);
  ALIAS (land_cover, land_cover_after);
 
@@ -112,6 +114,7 @@ $gdxin
  otherIntCost('pasture') = 0.8 + otherICost;
 
  SCALAR forestExpansionMaxRate / 1 /;
+ SCALAR restrictedDeforestationCost / 10000000 /;
 
  VARIABLES
        cropArea(crop, location)             total area for crops - Mha
@@ -143,7 +146,8 @@ $gdxin
 
  POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
                    totalFeedDM,
-                   landCoverArea, landCoverChange, woodHarvest;
+                   deforestedArea, reservedAreaDeforested, prematureDeforestationCost,
+                   landCoverArea, landCoverChange, woodHarvest, timberHarvest, timberForwardContract, totalConversionCost;
 
 * POSITIVE VARIABLE A;
 
@@ -264,7 +268,7 @@ $gdxin
                                                                                           sqrt(sqr(deforestedArea(forest, location) - [previousLandCoverArea(forest, location) - minimumLandCoverArea(forest, location)] + delta)) +
                                                                                                sqr(delta)] / 2;
 
- PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location) * woodYield('pasture', forest, location) * woodPrice);
+ PREMATURE_DEFORESTATION_COST_CALC .. prematureDeforestationCost =E= sum((forest, location), reservedAreaDeforested(forest, location)) * restrictedDeforestationCost;
 
 * UNMANAGED_FOREST_CONSTRAINT(location) .. sum((land_cover)$[not sameAs(land_cover, 'unmanagedForest')], landCoverChange(land_cover, 'unmanagedForest', location)) =E= 0;
 
@@ -272,13 +276,14 @@ $gdxin
 *                                                                         sum(land_cover_before, landCoverChange(land_cover_before, 'timberForest', location));
 
 * Timber from land cover change
- WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((land_cover_before, land_cover_after)$[not sameAs(land_cover_before, 'timberForest')], landCoverChange(land_cover_before, land_cover_after, location) * woodYield(land_cover_before, land_cover_after, location));
+ WOOD_HARVEST_CALC(location) .. woodHarvest(location) =E= sum((non_timber_forest_before, non_timber_forest_after), landCoverChange(non_timber_forest_before, non_timber_forest_after, location) * woodYield(non_timber_forest_before, non_timber_forest_after, location));
 
 * Timber from timberForest rotation
  TIMBER_HARVEST_CALC(location) .. timberHarvest(location) =E=  (sum(land_cover_after, landCoverChange('timberForest', land_cover_after, location)) - minimumLandCoverArea('timberForest', location)) * timberForestYield(location);
 
 * Future timber from newly planted timberForest
- TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(land_cover_before, landCoverChange(land_cover_before, 'timberForest', location) * woodYield(land_cover_before, 'timberForest', location));
+ TIMBER_FORWARD_CONTRACT_CALC(location) .. timberForwardContract(location) =E= sum(non_timber_forest_before, landCoverChange(non_timber_forest_before, 'timberForest', location) * woodYield(non_timber_forest_before, 'timberForest', location)) +
+                                                                               (landCoverChange('timberForest', 'timberForest', location) - minimumLandCoverArea('timberForest', location)) * woodYield('timberForest', 'timberForest', location);
 
  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));
 
@@ -292,7 +297,7 @@ $gdxin
          (
               (   SUM((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) ) * domesticPriceMarkup +
 
-               SUM(location, carbonFlux(location)) * carbonPrice - SUM(location, (woodHarvest(location) + timberForwardContract(location))) * woodPrice +
+               SUM(location, carbonFlux(location)) * carbonPrice - SUM(location, woodHarvest(location) + timberForwardContract(location)) * woodPrice +
 
                SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
 
@@ -325,6 +330,7 @@ $gdxin
  parameter feedCostRate(feed_crop);
  parameter productionShock(all_types);
  scalar netCarbonFlux;
+ scalar totalTimberProduction;
 
 * Production quantities based on smaller area (before unhandledCropArea adjustment applied)
  totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
@@ -346,6 +352,7 @@ $gdxin
  netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
 
  netCarbonFlux = SUM(location, carbonFlux.L(location));
+ totalTimberProduction = SUM(location, woodHarvest.L(location) + timberHarvest.L(location));
 
  Scalar totalCostsLessLU;
 
diff --git a/data/timber_demand.csv b/data/timber_demand.csv
new file mode 100644
index 0000000000000000000000000000000000000000..b352093be853a2cdb6915d715071b2e5d9baa2c6
--- /dev/null
+++ b/data/timber_demand.csv
@@ -0,0 +1,22 @@
+year,value
+2000,2458
+2005,2584
+2010,2715
+2015,2854
+2020,3000
+2025,3153
+2030,3313
+2035,3482
+2040,3660
+2045,3847
+2050,4043
+2055,4249
+2060,4466
+2065,4694
+2070,4933
+2075,5185
+2080,5450
+2085,5728
+2090,6020
+2095,6327
+2100,6650
diff --git a/debug_config.properties b/debug_config.properties
index 4c8cfb1366466445293bbcaf3f874f7e1464d55f..0bdf160eee80b1c21862ff49216f454ef389098c 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -9,7 +9,7 @@ START_TIMESTEP=0
 TIMESTEP_SIZE=1
 END_TIMESTEP=10
 
-IS_CALIBRATION_RUN=true
+IS_CALIBRATION_RUN=false
 END_FIRST_STAGE_CALIBRATION=10
 
 GENERATE_NEW_YIELD_CLUSTERS=false
@@ -20,6 +20,6 @@ DEBUG_LIMIT_COUNTRIES=true
 DEBUG_COUNTRY_NAME=Central America
 GAMS_COUNTRY_TO_SAVE=Central America
 
-WOOD_PRICE=200
+WOOD_PRICE=0.2
 CARBON_PRICE=0.02
-FOREST_LOCKIN_PERIOD=3
\ No newline at end of file
+FOREST_LOCKIN_PERIOD=30
\ No newline at end of file
diff --git a/src/ac/ed/lurg/CarbonMarket.java b/src/ac/ed/lurg/CarbonMarket.java
deleted file mode 100644
index 735a277fded5f0bfa7db4b2741de753fd842ad05..0000000000000000000000000000000000000000
--- a/src/ac/ed/lurg/CarbonMarket.java
+++ /dev/null
@@ -1,91 +0,0 @@
-package ac.ed.lurg;
-
-import java.io.FileInputStream;
-import java.io.FileOutputStream;
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import java.io.ObjectOutputStream;
-import java.util.Collection;
-import java.util.Map;
-import java.util.Map.Entry;
-
-import ac.ed.lurg.country.AbstractCountryAgent;
-import ac.ed.lurg.country.CountryAgent;
-import ac.ed.lurg.country.GlobalPrice;
-import ac.ed.lurg.landuse.CropUsageData;
-import ac.ed.lurg.types.CropToDoubleMap;
-import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.utils.LogWriter;
-
-public class CarbonMarket {
-	private GlobalPrice carbonPrice;
-	
-	public CarbonMarket() {
-		if (ModelConfig.IS_CALIBRATION_RUN) {
-			carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
-		} else {
-			carbonPrice = deserializeCarbonPrice();
-		}
-	}
-	
-	public GlobalPrice getCarbonPrice() {
-		return carbonPrice;
-	}
-	
-	/*
-	public void updateCarbonPrice(Collection<CountryAgent> countryAgents, double carbonDemand, Timestep timestep) {
-		double globalNetCarbonFlux;
-		for (CountryAgent ca : countryAgents) {
-			double netCarbonFlux = ca.getNetCarbonFlux();
-			globalNetCarbonFlux += netCarbonFlux;
-		}
-
-			GlobalPrice prevPrice = carbonPrice;
-			double imports = carbonDemand;
-			double exportsBeforeTransportLosses = globalNetCarbonFlux;
-			LogWriter.println(timestep.getYear() + " Updating carbon prices");
-			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exportsBeforeTransportLosses, timestep, globalNetCarbonFlux);
-			LogWriter.println( String.format("Price for carbon updated from %s to %s \n", prevPrice, adjustedPrice));
-			if (adjustedPrice.getStockLevel() < 0)
-				LogWriter.println("Global stocks are below zero carbon, " + timestep.getYear());
-
-			carbonPrice = adjustedPrice;
-	}
-	*/
-	
-	public void serializeGlobalPrices() {
-		try {
-			LogWriter.println("Starting serializing CarbonPrice to " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			ObjectOutputStream out = new ObjectOutputStream(fileOut);
-			out.writeObject(carbonPrice);
-			out.close();
-			fileOut.close();
-			LogWriter.println("Serialized carbon price is saved");
-		} catch (IOException i) {
-			i.printStackTrace();
-		}
-	}
-	
-	@SuppressWarnings("unchecked")
-	private GlobalPrice deserializeCarbonPrice() {
-		try {
-			GlobalPrice initCarbonPrice;
-			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			ObjectInputStream in = new ObjectInputStream(fileIn);
-			initCarbonPrice = (GlobalPrice) in.readObject();
-			in.close();
-			fileIn.close();
-			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			return initCarbonPrice;
-		} catch (IOException i) {
-			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			LogWriter.print(i);
-			return null;
-		} catch (ClassNotFoundException c) {
-			LogWriter.printlnError("GlobalPrice class not found");
-			c.printStackTrace();
-			return null;
-		}
-	}
-}
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index c069e5178955dab12f8028706ee698d607734b0f..028b0a233f746b352b12c88ff0a2e6f57a209a60 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -42,11 +42,13 @@ public class InternationalMarket {
 			}
 			
 			carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
+			timberPrice = GlobalPrice.createInitial(ModelConfig.WOOD_PRICE, 0);
 		}
 		else {
 			List<Object> deserializedPrices = deserializeGlobalPrice();
 			worldPrices = (Map<CropType, GlobalPrice>) deserializedPrices.get(0);
 			carbonPrice = (GlobalPrice) deserializedPrices.get(1);
+			timberPrice = (GlobalPrice) deserializedPrices.get(2);
 		}
 		priceShockManager = new PriceShockManager();
 	}
@@ -58,6 +60,10 @@ public class InternationalMarket {
 	public GlobalPrice getCarbonPrice() {
 		return carbonPrice;
 	}
+	
+	public GlobalPrice getTimberPrice() {
+		return timberPrice;
+	}
 
 	private Map<CropType, Double> getInitialStockLevels() {
 		Map<CropType, Double> initialStockLevels = new StockReader().read();
@@ -65,7 +71,7 @@ public class InternationalMarket {
 		return initialStockLevels;
 	}
 
-	void determineInternationalTrade(Collection<AbstractCountryAgent> countryAgents, double gen2EcDDemand, double carbonDemand, Timestep timestep) {
+	void determineInternationalTrade(Collection<AbstractCountryAgent> countryAgents, double gen2EcDDemand, double carbonDemand, double timberDemand, Timestep timestep) {
 		CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
 		CropToDoubleMap totalExportCommodities = new CropToDoubleMap();
 		CropToDoubleMap totalProduction = new CropToDoubleMap();
@@ -117,15 +123,28 @@ public class InternationalMarket {
 				totalCarbonSequestered -= netCarbonFlux;
 		}
 		
-		totalCarbonSequestered = Math.max(totalCarbonSequestered, 0.0000001);
-		GlobalPrice prevPrice = carbonPrice;
-				LogWriter.println(timestep.getYear() + " Updating carbon prices");
-		GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(carbonDemand, totalCarbonSequestered, timestep, totalCarbonSequestered, false);
-		LogWriter.println( String.format("Price for carbon updated from %s to %s \n", prevPrice, adjustedPrice));
-		if (adjustedPrice.getStockLevel() < 0)
+		totalCarbonSequestered = Math.max(totalCarbonSequestered, 0.0000001); // avoid division by 0
+		GlobalPrice prevCPrice = carbonPrice;
+		LogWriter.println(timestep.getYear() + " Updating carbon price");
+		GlobalPrice adjustedCPrice = prevCPrice.createWithUpdatedMarketPrices(carbonDemand, totalCarbonSequestered, timestep, totalCarbonSequestered, false);
+		LogWriter.println( String.format("Price for carbon updated from %s to %s \n", prevCPrice, adjustedCPrice));
+		if (adjustedCPrice.getStockLevel() < 0)
 			LogWriter.println("Global stocks are below zero carbon, " + timestep.getYear());
-
-		carbonPrice = adjustedPrice;
+		carbonPrice = adjustedCPrice;
+		
+		// Update timber price
+		double totalTimberProduced = 0;
+		for (AbstractCountryAgent ca : countryAgents) {
+			totalTimberProduced += ca.getTimberProduction();
+		}
+		totalTimberProduced = Math.max(totalTimberProduced, 0.0000001); // avoid division by 0
+		GlobalPrice prevTPrice = timberPrice;
+		LogWriter.println(timestep.getYear() + " Updating timber price");
+		GlobalPrice adjustedTPrice = prevTPrice.createWithUpdatedMarketPrices(timberDemand, totalTimberProduced, timestep, totalTimberProduced, false);
+		LogWriter.println( String.format("Price for timber updated from %s to %s \n", prevTPrice, adjustedTPrice));
+		if (adjustedTPrice.getStockLevel() < 0)
+			LogWriter.println("Global stocks are below zero timber, " + timestep.getYear());
+		timberPrice = adjustedTPrice;
 	}
 
 	void writeGlobalMarketFile(Timestep timestep, BufferedWriter outputFile) throws IOException {
@@ -154,6 +173,7 @@ public class InternationalMarket {
 		List<Object> pricesToSerialize = new ArrayList<Object>();
 		pricesToSerialize.add(worldPrices);
 		pricesToSerialize.add(carbonPrice);
+		pricesToSerialize.add(timberPrice);
 		
 		try {
 			LogWriter.println("Starting serializing GlobalPrice to " + ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE);
@@ -168,21 +188,6 @@ public class InternationalMarket {
 		}
 	}
 	
-/*
-	public void serializeCarbonPrices() {
-		try {
-			LogWriter.println("Starting serializing CarbonPrice to " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
-			ObjectOutputStream out = new ObjectOutputStream(fileOut);
-			out.writeObject(carbonPrice);
-			out.close();
-			fileOut.close();
-			LogWriter.println("Serialized carbon price is saved");
-		} catch (IOException i) {
-			i.printStackTrace();
-		}
-	}
-*/
 	@SuppressWarnings("unchecked")
 	private List<Object> deserializeGlobalPrice() {
 		try {
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 583b7be0f835c0680c86ce95eb178c62e52eda70..7009cb271a6cb227cd3f37d7e30ac0ea63785dc8 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -19,8 +19,7 @@ public class ModelConfig {
 	
 	private Properties configFile;
 	private static ModelConfig modelConfig;
-	//public static final String CONFIG_FILE = System.getProperty("CONFIG_FILE");
-	public static final String CONFIG_FILE = "C:\\Users\\Bart\\git\\plumv2\\debug_config.properties";
+	public static final String CONFIG_FILE = System.getProperty("CONFIG_FILE");
 	
 	private static StringTabularReader shocksReader;
 	
@@ -429,7 +428,9 @@ public class ModelConfig {
 	public static final String CARBON_DEMAND_FILENAME = getProperty("CARBON_DEMAND_FILENAME", "carbon_demand.csv");
 	public static final String CARBON_DEMAND_FILE = getProperty("CARBON_DEMAND_FILE", DATA_DIR + File.separator + CARBON_DEMAND_FILENAME);
 	public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 0.02); // $1000/tC-eq
-	public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.05); // $1000/tC-eq
+	public static final String TIMBER_DEMAND_FILENAME = getProperty("TIMBER_DEMAND_FILENAME", "timber_demand.csv");
+	public static final String TIMBER_DEMAND_FILE = getProperty("TIMBER_DEMAND_FILE", DATA_DIR + File.separator + TIMBER_DEMAND_FILENAME);
+	public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.2); // $1000/tC-eq
 	public static final int FOREST_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
 	
 }
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 04a7ca553d3a2b8f01614efc118a8be03792a1cf..ca92fa0f200af50fe050088e80322cd97c18f573 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -74,7 +74,6 @@ public class ModelMain {
 	private RasterHeaderDetails desiredProjection;
 
 	private InternationalMarket internationalMarket;
-	//private CarbonMarket carbonMarket;
 	private IrrigationRasterSet currentIrrigationData;
 	private RasterSet<LandUseItem> globalLandUseRaster;
 	private RasterSet<IntegerRasterItem> clusterIdRaster;
@@ -111,7 +110,6 @@ public class ModelMain {
 
 		globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
 		internationalMarket = new InternationalMarket();
-		//carbonMarket = new CarbonMarket();
 		createCountryAgents(compositeCountryManager.getAll());
 	}
 
@@ -148,20 +146,25 @@ public class ModelMain {
 		double carbonDemand = demandManager.getCarbonDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
 		double carbonDemandIncrease = (carbonDemand > previousCarbonDemand) ? carbonDemand - previousCarbonDemand: 0;
 		
+		double previousTimberDemand = (timestep.isInitialTimestep() || ModelConfig.IS_CALIBRATION_RUN ) ? 0: demandManager.getTimberDemand(timestep.getPreviousTimestep());
+		double timberDemand = demandManager.getTimberDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
+		double timberDemandIncrease = (timberDemand > previousTimberDemand) ? timberDemand - previousTimberDemand: 0;
+		
 		CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
 		WoodYieldRasterSet currentWoodYieldData = getWoodYieldData(timestep);
 		
 		DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new ConversionCostReader().read();
 		
-		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData, conversionCosts, carbonDemandIncrease);
-		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep);
+		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, 
+				currentWoodYieldData, conversionCosts, carbonDemandIncrease, timberDemandIncrease);
+		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timberDemand, timestep);
 		
 		// loop for iterations.  Could check within a tolerance using internationalMarket.findMaxPriceDiff, not doing so as doesn't find a solution due to inelastic consumption
 		for (int i=0; i < ModelConfig.DEMAND_RECALC_MAX_ITERATIONS; i++) {
 			LogWriter.println("\n++ Re-estimating prices and demand " + i);
 			countryAgents.recalculateDemandForAll(); // recalculate demand from new prices
 			countryAgents.updateNetImportsForAll(); // calculate imports and exports
-			internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep); // calculate prices
+			internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timberDemand, timestep); // calculate prices
 		}
 		internationalMarket.applyPriceShocks(timestep);
 		
@@ -534,10 +537,11 @@ public class ModelMain {
 	}
 
 	public void createCountryAgents(Collection<CompositeCountry> countryGrouping) {
-		countryAgents = new CountryAgentManager(compositeCountryManager, demandManager, countryBoundaryRaster, internationalMarket, clusterIdRaster, globalLandUseRaster);
+		WoodYieldRasterSet initWoodYieldData = getWoodYieldData(new Timestep(0));
+		countryAgents = new CountryAgentManager(compositeCountryManager, demandManager, countryBoundaryRaster, internationalMarket, clusterIdRaster, globalLandUseRaster, initWoodYieldData);
 		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = getInitialCropUsageData();
 		RasterSet<LandUseItem> initLU = getInitialLandUse();
-
+		
 		for (CompositeCountry cc : countryGrouping) {
 			countryAgents.addForCountry(cc, cropUsageDataMap, initLU);
 		}
diff --git a/src/ac/ed/lurg/country/AbstractCountryAgent.java b/src/ac/ed/lurg/country/AbstractCountryAgent.java
index 879373879f462668ce74a2ecccd2dd486726509e..8aa5d9c87cc0256fe372206449cb3ce96b08bf14 100644
--- a/src/ac/ed/lurg/country/AbstractCountryAgent.java
+++ b/src/ac/ed/lurg/country/AbstractCountryAgent.java
@@ -19,6 +19,7 @@ public abstract class AbstractCountryAgent {
 	protected Map<CommodityType, Double> currentProjectedDemand;
 	protected Map<CropType, CountryPrice> currentCountryPrices;
 	protected CountryPrice currentCarbonPrice;
+	protected CountryPrice currentTimberPrice;
 	protected Map<CommodityType, Map<CropType, Double>> baseDemandFact;
 	protected Timestep currentTimestep;
 	protected Map<CommodityType, Map<CropType, Double>> currentMinDemandFract;
@@ -40,7 +41,7 @@ public abstract class AbstractCountryAgent {
 		currentMinDemandFract = getMinDemandFraction();
 	}
 
-	private void calculateCountryPrices(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice){
+	private void calculateCountryPrices(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice){
 		Map<CropType, CountryPrice> countryPrices = new HashMap <CropType, CountryPrice>();
 
 		for (CropType c : CropType.getImportedTypes()) {
@@ -53,6 +54,7 @@ public abstract class AbstractCountryAgent {
 
 		currentCountryPrices = countryPrices;
 		currentCarbonPrice = new CountryPrice(0.0, carbonPrice.getExportPrice());
+		currentTimberPrice = new CountryPrice(0.0, timberPrice.getExportPrice());
 	}
 	
 	private void calculateDemand(boolean outputGamsDemand) {
@@ -60,8 +62,8 @@ public abstract class AbstractCountryAgent {
 		currentProjectedDemand = demandManager.getDemand(country, currentTimestep.getYear(), producerPrices, outputGamsDemand);
 	}
 	
-	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, boolean outputGamsDemand) {
-		calculateCountryPrices(worldPrices, carbonPrice);
+	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice, boolean outputGamsDemand) {
+		calculateCountryPrices(worldPrices, carbonPrice, timberPrice);
 		calculateDemand(outputGamsDemand);
 	}
 	
@@ -174,4 +176,6 @@ public abstract class AbstractCountryAgent {
 	
 	abstract public double getNetCarbonFlux();
 	
+	abstract public double getTimberProduction();
+	
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 5baf3d7eb7525b9e3521f0bfc0afce394fa899b0..d1015b72be9666a8bbc5f6ea351c4c7e246b010c 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -25,6 +25,7 @@ import ac.ed.lurg.landuse.CropUsageData;
 import ac.ed.lurg.landuse.IrrigationItem;
 import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.landuse.WoodYieldItem;
+import ac.ed.lurg.landuse.WoodYieldRasterSet;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
@@ -52,7 +53,7 @@ public class CountryAgent extends AbstractCountryAgent {
 
 	public CountryAgent(AbstractDemandManager demandManager,CompositeCountry country, RasterSet<LandUseItem> cropAreaRaster,
 			Map<CropType, CropUsageData> cropUsageData, Map<CropType, Double> tradeBarriers, RasterSet<IntegerRasterItem> yieldClusters,
-			Map<CropType, Double> subsidyRates) {
+			Map<CropType, Double> subsidyRates, WoodYieldRasterSet initWoodYieldData) {
 
 		super(demandManager, country, tradeBarriers);
 		this.yieldClusters = yieldClusters;
@@ -63,7 +64,8 @@ public class CountryAgent extends AbstractCountryAgent {
 		
 		saveGamsGdxFiles = (ModelConfig.GAMS_COUNTRY_TO_SAVE != null && country.getName().equals(ModelConfig.GAMS_COUNTRY_TO_SAVE));
 		
-		initialiseForestManager(yieldClusters);
+		if (yieldClusters != null)
+			initialiseForestManager(initWoodYieldData);
 		
 	}
 
@@ -117,10 +119,10 @@ public class CountryAgent extends AbstractCountryAgent {
 	public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData, 
 			Map<CropType, GlobalPrice> worldPrices, double globalGen2EcIncrease, RasterSet<CarbonFluxItem> carbonFluxData, 
 			RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
-			double carbonDemandIncrease, GlobalPrice carbonPrice) {
+			double carbonDemandIncrease, GlobalPrice carbonPrice, double timberDemandIncrease, GlobalPrice timberPrice, WoodYieldRasterSet initWoodYieldData) {
 
 		// project demand
-		calculateCountryPricesAndDemand(worldPrices, carbonPrice, false);
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, timberPrice, false);
 		
 		if (saveGamsGdxFiles && ModelConfig.PRICE_ELASTIC_DEMAND)
 			saveGDXFile("demand");
@@ -134,11 +136,14 @@ public class CountryAgent extends AbstractCountryAgent {
 		else if (ModelConfig.DEBUG_JUST_DEMAND_OUTPUT) { // if this debug flag is set we don't do the optimisation
 		}
 		else {			
-			if (yieldClusters==null)
+			if (yieldClusters==null) {
 				yieldClusters = calcYieldClusters(irrigData, countryYieldSurfaces);  // this should only be on the first timestep
+				initialiseForestManager(initWoodYieldData);
+			}
 			
 			// optimize areas and intensity 
-			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, woodYieldData, conversionCosts, carbonDemandIncrease);
+			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, woodYieldData, conversionCosts, 
+					carbonDemandIncrease, timberDemandIncrease);
 			GamsRasterOptimiser opti = new GamsRasterOptimiser(input, yieldClusters);
 			LogWriter.println("Running " + country.getName() + ", currentTimestep " + currentTimestep);
 
@@ -177,7 +182,7 @@ public class CountryAgent extends AbstractCountryAgent {
 
 	private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces, double gen2EcIncrease,
 			RasterSet<CarbonFluxItem> carbonFluxData, RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
-			double carbonDemandIncrease) {
+			double carbonDemandIncrease, double timberDemandIncrease) {
 		double allowedImportChange;
 
 		if (currentTimestep.isInitialTimestep() || (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION)) {  // initialisation time-step
@@ -222,22 +227,37 @@ public class CountryAgent extends AbstractCountryAgent {
 		}
 		
 		// Carbon import/export constraints TODO not used
-		double baseTrade = getNetCarbonFlux();
-		double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
-		double changeUp = 0.0;
-		if (allowedImportChange > 0.0)
-			changeUp = Math.max(baseTrade * allowedImportChange, carbonDemandIncrease * countryArea / 1500 * 2);
-			
-		TradeConstraint carbonTradeConstraint = new TradeConstraint(0, baseTrade + changeUp);
+		TradeConstraint carbonTradeConstraint;
+		{
+			double baseTrade = getNetCarbonFlux();
+			double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
+			double change = 0.0;
+			if (allowedImportChange > 0.0)
+				change = Math.max(baseTrade * allowedImportChange, carbonDemandIncrease * countryArea / 1500 * 2);
+			carbonTradeConstraint = new TradeConstraint(baseTrade - change, baseTrade + change);
+		}
+
+		// Timber import/export constraints
+		TradeConstraint timberTradeConstraint;
+		{
+			double baseTrade = getTimberProduction();
+			double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
+			double change = 0.0;
+			if (allowedImportChange > 0.0)
+				change = Math.max(baseTrade * allowedImportChange, timberDemandIncrease * countryArea / 1500 * 2);
+
+			timberTradeConstraint = new TradeConstraint(baseTrade - change, baseTrade + change);
+		}
 		
 		// Aggregate minimum land cover for current timestep
 		DoubleMap<Integer, LandCoverType, Double> minimumLandCoverForTimestep = forestManager.getReservedAreaForTimestep(currentTimestep); //<Location, LandCoverType, Area>
 		DoubleMap<Integer, LandCoverType, Double> timberYieldForTimestep = forestManager.getTimberYieldForTimestep(currentTimestep); //<Location, LandCoverType, Area>
+		Map<Integer, Double> naturalYieldFactorForTimestep = forestManager.getNaturalYieldFactorForTimestep(currentTimestep);
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
-				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates, currentCarbonPrice, carbonTradeConstraint);	
+				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates, currentCarbonPrice, carbonTradeConstraint, currentTimberPrice, timberTradeConstraint);	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, 
-				carbonFluxData, woodYieldData, minimumLandCoverForTimestep, timberYieldForTimestep, conversionCosts, countryLevelInputs);
+				carbonFluxData, woodYieldData, minimumLandCoverForTimestep, timberYieldForTimestep, naturalYieldFactorForTimestep, conversionCosts, countryLevelInputs);
 
 		return input;
 	}
@@ -324,57 +344,18 @@ public class CountryAgent extends AbstractCountryAgent {
 		}
 		
 	}
-	
-	@SuppressWarnings("serial")
-	private void initialiseForestManager(RasterSet<IntegerRasterItem> yieldClusters) {
+
+	private void initialiseForestManager(WoodYieldRasterSet initWoodYieldData) {
 		forestManager = new ForestManager();
-		
-		RasterSet<LandUseItem> initialLandUse = previousGamsRasterOutput.getLandUses();
-		
-		LazyTreeMap<Integer, Double> unmanagedForestArea = new LazyTreeMap<Integer, Double>() {
-			protected Double createValue() {return 0.0;}
-		};
-		LazyTreeMap<Integer, Double> otherNaturalArea = new LazyTreeMap<Integer, Double>() {
-			protected Double createValue() {return 0.0;}
-		};
-		
-		for (Entry<RasterKey, LandUseItem> entry : initialLandUse.entrySet()) {
-			RasterKey key = entry.getKey();
-			LandUseItem luItem = entry.getValue();
-			int clusterId = yieldClusters.get(key).getInt();
-			double unmanagedForestAreaSoFar = unmanagedForestArea.lazyGet(clusterId);
-			double otherNaturalAreaSoFar = otherNaturalArea.lazyGet(clusterId);
-			double forestarea = luItem.getInitialForestedNaturalArea();
-			double natarea = luItem.getLandCoverArea(LandCoverType.NATURAL) ;
-			unmanagedForestArea.put(clusterId, unmanagedForestAreaSoFar + luItem.getInitialForestedNaturalArea());
-			otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getLandCoverArea(LandCoverType.NATURAL) - luItem.getInitialForestedNaturalArea());
-		}
-		
-		int otherNaturalPseudoYear = ModelConfig.BASE_YEAR;
-		int unmanagedForestPseudoYear = ModelConfig.BASE_YEAR - ModelConfig.FOREST_LOCKIN_PERIOD;
-		
-		for (Entry<Integer, Double> entry : unmanagedForestArea.entrySet()) {
-			forestManager.addForestRecord(new ForestRecord(unmanagedForestPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
-		}
-		
-		for (Entry<Integer, Double> entry : otherNaturalArea.entrySet()) {
-			forestManager.addForestRecord(new ForestRecord(otherNaturalPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
-		}
-		
-		/*
-		int maxYear = Integer.MAX_VALUE;
-		for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
-			
-			for (IntegerRasterItem iri : yieldClusters.values()) {
-				// inefficient, better to store somewhere the total number of clusters
-				if (iri != null) // TODO Why?				
-					forestManager.put(maxYear, iri.getInt(), lc, 0.0);
-			}
-		}
-		*/
+		// Add initial natural areas and managed forest. Assume that unmanaged forest is a mature natural area (climax community) and other natural is a new natural area
+		forestManager.initialiseForestHistory(yieldClusters, previousGamsRasterOutput.getLandUses(), initWoodYieldData);
 	}
 	
 	public double getNetCarbonFlux() {
 		return previousGamsRasterOutput.getNetCarbonFlux();
 	}
+	
+	public double getTimberProduction() {
+		return previousGamsRasterOutput.getTimberProduction();
+	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index 3ca85261846c5cedd72f431c4266475f833a5b05..539e7692ee2220be3528f6504cfd55086ccedf60 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -9,7 +9,6 @@ import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
 
-import ac.ed.lurg.CarbonMarket;
 import ac.ed.lurg.InternationalMarket;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.Timestep;
@@ -41,6 +40,7 @@ public class CountryAgentManager {
 	private CountryBoundaryRaster countryBoundaryRaster;
 	private RasterSet<IntegerRasterItem> clusterIdRaster;
 	private RasterSet<LandUseItem> globalLandUseRaster;
+	private WoodYieldRasterSet initWoodYieldData;
 	
 	private Collection<AbstractCountryAgent> countryAgents = new ArrayList<AbstractCountryAgent>();
 	private Collection<CountryAgent> gamsCountryAgents = new ArrayList<CountryAgent>();
@@ -49,12 +49,13 @@ public class CountryAgentManager {
 	
 	public CountryAgentManager (CompositeCountryManager compositeCountryManager, AbstractDemandManager demandManager, 
 			CountryBoundaryRaster countryBoundaryRaster, InternationalMarket internationalMarket,
-			RasterSet<IntegerRasterItem> clusterIdRaster, RasterSet<LandUseItem> globalLandUseRaster) {
+			RasterSet<IntegerRasterItem> clusterIdRaster, RasterSet<LandUseItem> globalLandUseRaster, WoodYieldRasterSet initWoodYieldData) {
 		this.demandManager = demandManager;
 		this.countryBoundaryRaster = countryBoundaryRaster;
 		this.clusterIdRaster = clusterIdRaster;
 		this.internationalMarket = internationalMarket;
 		this.globalLandUseRaster = globalLandUseRaster;
+		this.initWoodYieldData = initWoodYieldData;
 		tradeBarrierManager = new TradeManager(compositeCountryManager);
 		subsidyRateManager = new SubsidyRateManager(compositeCountryManager);
 		craftyManager = new CraftyProdManager();
@@ -91,7 +92,7 @@ public class CountryAgentManager {
 				else {
 					Map<CropType, Double> subsidyRates = subsidyRateManager.getSubsidyRates(cc);
 					
-					CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, tradeBarriers, yieldClusters, subsidyRates);
+					CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, tradeBarriers, yieldClusters, subsidyRates, initWoodYieldData);
 					gamsCountryAgents.add(ca);
 					countryAgents.add(ca);
 				}
@@ -105,7 +106,7 @@ public class CountryAgentManager {
 
 	public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData, double gen2Increase,
 			CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
-			double carbonDemandIncrease) {
+			double carbonDemandIncrease, double timberDemandIncrease) {
 		
 		for (AbstractCountryAgent aca : countryAgents) {		
 			aca.setCurrentTimestep(timestep);
@@ -121,7 +122,7 @@ public class CountryAgentManager {
 			
 			try {
 				ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase, carbonFluxData, woodYieldData, 
-						conversionCosts, carbonDemandIncrease, internationalMarket.getCarbonPrice());
+						conversionCosts, carbonDemandIncrease, internationalMarket.getCarbonPrice(), timberDemandIncrease, internationalMarket.getTimberPrice(), initWoodYieldData);
 				
 				// update global rasters
 				globalLandUseRaster.putAll(ca.getLandUses());
@@ -137,7 +138,8 @@ public class CountryAgentManager {
 		}
 		
 		if (craftyCountryAgents.size() > 0) {
-			craftyManager.updateWithCraftyData(craftyCountryAgents, timestep, internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice());  // this will wait for the marker file from CRAFTY
+			craftyManager.updateWithCraftyData(craftyCountryAgents, timestep, internationalMarket.getWorldPrices(), 
+					internationalMarket.getCarbonPrice(), internationalMarket.getTimberPrice());  // this will wait for the marker file from CRAFTY
 		}
 	}
 	
@@ -150,7 +152,7 @@ public class CountryAgentManager {
 	public void recalculateDemandForAll() {
 		for (CountryAgent ca : gamsCountryAgents) {
 			ca.shockGDP();
-			ca.calculateCountryPricesAndDemand(internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice(), true);
+			ca.calculateCountryPricesAndDemand(internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice(), internationalMarket.getTimberPrice(), true);
 		}
 	}
 	
diff --git a/src/ac/ed/lurg/country/ForestManager.java b/src/ac/ed/lurg/country/ForestManager.java
index fc2c8439437927a5a53a9e9b774b2d66f479d2f4..49e264d6e940b6d0a13ed97c90371aea4ac90fc5 100644
--- a/src/ac/ed/lurg/country/ForestManager.java
+++ b/src/ac/ed/lurg/country/ForestManager.java
@@ -10,9 +10,16 @@ import java.util.stream.Collectors;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.Timestep;
+import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.landuse.WoodYieldRasterSet;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.utils.DoubleMap;
+import ac.ed.lurg.utils.LazyTreeMap;
 import ac.ed.lurg.utils.LogWriter;
+import ac.sac.raster.IntegerRasterItem;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+import ac.ed.lurg.landuse.WoodYieldItem;
 
 public class ForestManager {
 	private List<ForestRecord> forestHistory;
@@ -25,6 +32,78 @@ public class ForestManager {
 		forestHistory.add(new ForestRecord(record.getYearPlanted(), record.getLocation(), record.getForestType(), record.getArea(), record.getPotentialYield()));	
 	}
 	
+	public void initialiseForestHistory(RasterSet<IntegerRasterItem> yieldClusters, RasterSet<LandUseItem> initialLandUse, WoodYieldRasterSet initWoodYieldData) {
+		LazyTreeMap<Integer, Double> unmanagedForestArea = new LazyTreeMap<Integer, Double>() {
+			private static final long serialVersionUID = 1L;
+			protected Double createValue() {return 0.0;}
+		};
+		LazyTreeMap<Integer, Double> otherNaturalArea = new LazyTreeMap<Integer, Double>() {
+			private static final long serialVersionUID = 1L;
+			protected Double createValue() {return 0.0;}
+		};
+		LazyTreeMap<Integer, Double> timberForestArea = new LazyTreeMap<Integer, Double>() {
+			private static final long serialVersionUID = 1L;
+			protected Double createValue() {return 0.0;}
+		};
+		LazyTreeMap<Integer, Double> timberYieldTimesArea = new LazyTreeMap<Integer, Double>() {
+			private static final long serialVersionUID = 1L;
+			protected Double createValue() {return 0.0;}
+		};
+		
+		for (Entry<RasterKey, LandUseItem> entry : initialLandUse.entrySet()) {
+			RasterKey key = entry.getKey();
+			LandUseItem luItem = entry.getValue();
+			if (luItem == null) {
+				LogWriter.printlnWarning("CountryAgent.initialiseForestManager: null luItem");
+				continue;
+			}
+			try {
+				int clusterId = yieldClusters.get(key).getInt();
+
+				double unmanagedForestAreaSoFar = unmanagedForestArea.lazyGet(clusterId);
+				double otherNaturalAreaSoFar = otherNaturalArea.lazyGet(clusterId);
+				double timberForestAreaSoFar = timberForestArea.lazyGet(clusterId);
+				double timberYieldSoFar = timberYieldTimesArea.lazyGet(clusterId);
+				//LogWriter.println(key.toString());
+				WoodYieldItem wYItem = initWoodYieldData.get(key);
+				double timberYieldThisTime;
+				timberYieldThisTime = (wYItem != null) ? initWoodYieldData.get(key).getWoodYield(LandCoverType.PASTURE, LandCoverType.TIMBER_FOREST) : 0.0;
+		
+				unmanagedForestArea.put(clusterId, unmanagedForestAreaSoFar + luItem.getInitialForestedNaturalArea());
+				otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getLandCoverArea(LandCoverType.NATURAL) - luItem.getInitialForestedNaturalArea()); // Natural is other natural + unmanaged forest so need to subtract the latter
+				timberForestArea.put(clusterId, timberForestAreaSoFar + luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST));
+				timberYieldTimesArea.put(clusterId, timberYieldSoFar + luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST) * timberYieldThisTime);
+				
+				/*
+				if (clusterId == 3) {
+					LogWriter.println(""+timberYieldThisTime+","+timberYieldSoFar+","+luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST)+","+timberForestAreaSoFar);
+				}*/
+
+			} catch (Exception e) {
+				LogWriter.print(e);
+			}
+		}
+		
+		// Assume unmanaged forest is fully mature (max wood yield) and other natural is new (no wood yield).
+		int otherNaturalPseudoYear = ModelConfig.BASE_YEAR;
+		int forestPseudoYear = ModelConfig.BASE_YEAR - ModelConfig.FOREST_LOCKIN_PERIOD;
+		
+		for (Entry<Integer, Double> entry : unmanagedForestArea.entrySet()) {
+			addForestRecord(new ForestRecord(forestPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
+		}
+		
+		for (Entry<Integer, Double> entry : otherNaturalArea.entrySet()) {
+			addForestRecord(new ForestRecord(otherNaturalPseudoYear, entry.getKey(), LandCoverType.NATURAL, entry.getValue(), 0));
+		}
+		
+		for (Entry<Integer, Double> entry : timberForestArea.entrySet()) {
+			int clusterId = entry.getKey();
+			double area = entry.getValue();
+			double timberYield = (area > 0) ? timberYieldTimesArea.get(clusterId) / area : 0.0;
+			addForestRecord(new ForestRecord(forestPseudoYear, entry.getKey(), LandCoverType.TIMBER_FOREST, area, timberYield));
+		}
+	}
+	
 	public DoubleMap<Integer, LandCoverType, Double> getReservedAreaForTimestep(Timestep timestep) {
 		DoubleMap<Integer, LandCoverType, Double> minimumLandCoverForTimestep = new DoubleMap<Integer, LandCoverType, Double>();
 		double minYear = timestep.getYear() - ModelConfig.FOREST_LOCKIN_PERIOD;
@@ -38,31 +117,33 @@ public class ForestManager {
 	}
 	
 	public DoubleMap<Integer, LandCoverType, Double> getTimberYieldForTimestep(Timestep timestep) {
-		DoubleMap<Integer, LandCoverType, Double> timberYield = new DoubleMap<Integer, LandCoverType, Double>();
+		DoubleMap<Integer, LandCoverType, Double> yieldTimesArea = new DoubleMap<Integer, LandCoverType, Double>();
 		DoubleMap<Integer, LandCoverType, Double> area = new DoubleMap<Integer, LandCoverType, Double>();
 		double yearPlanted = timestep.getYear() - ModelConfig.FOREST_LOCKIN_PERIOD;
 		for (ForestRecord record : forestHistory) {
 			if (record.getYearPlanted() == yearPlanted && record.getForestType().isManagedForest()) {
-				timberYield.addTo(record.getLocation(), record.getForestType(), record.getPotentialYield() * record.getArea());
+				yieldTimesArea.addTo(record.getLocation(), record.getForestType(), record.getPotentialYield() * record.getArea());
 				area.addTo(record.getLocation(), record.getForestType(), record.getArea());
 			}
 		}
 		
+		DoubleMap<Integer, LandCoverType, Double> yield = new DoubleMap<Integer, LandCoverType, Double>();
 		// Weighted average
-		for (Entry<Integer, Map<LandCoverType, Double>> locMap : timberYield.getMap().entrySet()) {
+		for (Entry<Integer, Map<LandCoverType, Double>> locMap : yieldTimesArea.getMap().entrySet()) {
 			for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
-				double y = timberYield.get(locMap.getKey(), coverMap.getKey());
+				double ya = yieldTimesArea.get(locMap.getKey(), coverMap.getKey());
 				double a = area.get(locMap.getKey(), coverMap.getKey());
-				timberYield.put(locMap.getKey(), coverMap.getKey(), y / a);
+				double y = (a > 0) ? ya / a : a;
+				yield.put(locMap.getKey(), coverMap.getKey(), y);
 			}
 		}
 		
-		return timberYield;
+		return yield;
 	}
 	
 	// Calculates wood yield multiplier for natural areas based on age
 	double naturalWoodYieldFunction(double age) {
-		return Math.max(age, ModelConfig.FOREST_LOCKIN_PERIOD) / ModelConfig.FOREST_LOCKIN_PERIOD;
+		return Math.min(age, ModelConfig.FOREST_LOCKIN_PERIOD) / ModelConfig.FOREST_LOCKIN_PERIOD;
 	}
 	
 	public Map<Integer, Double> getNaturalYieldFactorForTimestep(Timestep timestep) {
@@ -82,7 +163,12 @@ public class ForestManager {
 			for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
 				double y = ageTimesArea.get(locMap.getKey(), coverMap.getKey());
 				double a = area.get(locMap.getKey(), coverMap.getKey());
-				yieldFactor.put(locMap.getKey(), naturalWoodYieldFunction(y / a));
+				if (a != 0) {
+					yieldFactor.put(locMap.getKey(), naturalWoodYieldFunction(y / a));
+				} else {
+					yieldFactor.put(locMap.getKey(), 0.0);
+				}
+				
 			}
 		}
 
diff --git a/src/ac/ed/lurg/country/GlobalPrice.java b/src/ac/ed/lurg/country/GlobalPrice.java
index 163eb9dadd79061430c65f18d00a98ee1cbe4cd4..c71ff9b6b09426f8c8d49035179919e4b4926634 100644
--- a/src/ac/ed/lurg/country/GlobalPrice.java
+++ b/src/ac/ed/lurg/country/GlobalPrice.java
@@ -104,6 +104,10 @@ public class GlobalPrice implements Serializable {
 				else // stock decreasing, so increase price.  stockChange is negative so adjustment > 1
 					adjustment = 1 - ModelConfig.MARKET_LAMBA * stockChange/Math.max(0, stockLevel);
 				
+				if (Double.isNaN(adjustment)) {
+					int foo = 1;
+				}
+				
 				LogWriter.print(String.format("     initally adjustment= %.4f", adjustment));
 			}
 			
diff --git a/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java b/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
index 8879fd0d1293db6b66a5dc39ff4c8d44c17abf1b..dc5d375319f53186bbc92df4b924cb0738090200 100644
--- a/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
+++ b/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
@@ -28,8 +28,8 @@ public class CraftyCountryAgent extends AbstractCountryAgent {
 		return cropUsageData;
 	}
 
-	public void updateProduction(Map<CropType, Double> cropProduction, Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice) {
-		calculateCountryPricesAndDemand(worldPrices, carbonPrice, false);
+	public void updateProduction(Map<CropType, Double> cropProduction, Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice) {
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, timberPrice, false);
 
 		cropUsageData = new HashMap<CropType, CropUsageData>();
 		for (Entry<CropType, Double> entry : cropProduction.entrySet()) {
@@ -54,4 +54,8 @@ public class CraftyCountryAgent extends AbstractCountryAgent {
 	public double getNetCarbonFlux() {
 		return 0;
 	}
+	
+	public double getTimberProduction() {
+		return 0;
+	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/crafty/CraftyProdManager.java b/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
index 7b67dccf625f129bb423ee309b458bd60d60fd31..21a9988ac39629afa43b103b4a6100032faa98ac 100644
--- a/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
+++ b/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
@@ -32,7 +32,8 @@ public class CraftyProdManager {
 		return craftyCountries;
 	}
 
-	public void updateWithCraftyData(Collection<CraftyCountryAgent> craftyCountryAgents, Timestep timestep, Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice) {
+	public void updateWithCraftyData(Collection<CraftyCountryAgent> craftyCountryAgents, Timestep timestep, Map<CropType, GlobalPrice> worldPrices, 
+			GlobalPrice carbonPrice, GlobalPrice timberPrice) {
 		String rootDir = ModelConfig.CRAFTY_PRODUCTION_DIR + File.separator + timestep.getYear();
 		long startTime = System.currentTimeMillis();
 		
@@ -66,7 +67,7 @@ public class CraftyProdManager {
 				if (cropProduction.size() < CropType.getImportedTypes().size()) {  // Don't need setaside or pasture, which aren't imported either
 					LogWriter.printlnError("Not all crops present in Crafty production for country: " + cca.getCountry() + " only " + cropProduction.size());
 				}
-				cca.updateProduction(cropProduction, worldPrices, carbonPrice);
+				cca.updateProduction(cropProduction, worldPrices, carbonPrice, timberPrice);
 			}
 			catch (Exception e) {
 				LogWriter.println("Problem getting Crafty data for: " + cca.getCountry());
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index f0fec4b8a08e670155076130ce97ed4d5cf7a732..fd47fa0c042bd85f38c04ad1e034c369421ce4ea 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -22,11 +22,13 @@ public class GamsCountryInput {
 	private Map<CropType, Double> subsidyRates;
 	private CountryPrice carbonPrice;
 	private TradeConstraint carbonTradeConstraint;
+	private CountryPrice timberPrice;
+	private TradeConstraint timberTradeConstraint;
 
 	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,
-			CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint) {
+			CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint, CountryPrice timberPrice, TradeConstraint timberTradeConstraint) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
@@ -37,6 +39,8 @@ public class GamsCountryInput {
 		this.subsidyRates = subsidyRates;
 		this.carbonPrice = carbonPrice;
 		this.carbonTradeConstraint = carbonTradeConstraint;
+		this.timberPrice= timberPrice;
+		this.timberTradeConstraint = timberTradeConstraint;
 	}
 		
 	public CompositeCountry getCountry() { 
@@ -96,4 +100,16 @@ public class GamsCountryInput {
 	public TradeConstraint getCarbonTradeConstraint() {
 		return carbonTradeConstraint;
 	}
+
+	public CountryPrice getTimberPrice() {
+		return timberPrice;
+	}
+
+	public TradeConstraint getTimberTradeConstraint() {
+		return timberTradeConstraint;
+	}
+
+	public void setTradeConstraints(Map<CropType, TradeConstraint> tradeConstraints) {
+		this.tradeConstraints = tradeConstraints;
+	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationInput.java b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
index e9f9f8d42dcc1fa9f57432702361e0f44f294cc8..669fd0f2dca089a2f2a367965c1cab18058e1f72 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationInput.java
@@ -21,13 +21,14 @@ public class GamsLocationInput {
 	private Map<Integer, ? extends WoodYieldItem> woodYields;
 	private DoubleMap<Integer, LandCoverType, Double> minimumLandCover;
 	private DoubleMap<Integer, LandCoverType, Double> timberYield;
+	private Map<Integer, Double> naturalWoodYieldFactor;
 	private DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts;
 	private GamsCountryInput countryInput;
 	
 	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 WoodYieldItem> woodYields, DoubleMap<Integer, LandCoverType, Double> minimumLandCover, 
-			DoubleMap<Integer, LandCoverType, Double> timberYield,
+			DoubleMap<Integer, LandCoverType, Double> timberYield, Map<Integer, Double> naturalWoodYieldFactor,
 			DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts, GamsCountryInput countryInput) {
 		super();
 		this.timestep = timestep;
@@ -38,6 +39,7 @@ public class GamsLocationInput {
 		this.woodYields = woodYields;
 		this.minimumLandCover = minimumLandCover;
 		this.timberYield = timberYield;
+		this.naturalWoodYieldFactor = naturalWoodYieldFactor;
 		this.conversionCosts = conversionCosts;
 		this.countryInput = countryInput;	
 	}
@@ -70,6 +72,10 @@ public class GamsLocationInput {
 		return timberYield;
 	}
 	
+	public Map<Integer, Double> getNaturalWoodYieldFactor() {
+		return naturalWoodYieldFactor;
+	}
+	
 	public DoubleMap<LandCoverType, LandCoverType, Double> getConversionCosts() {
 		return conversionCosts;
 	}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index aaef758990a6349c0c9ad159ab679402eb4c1b09..c426ebd7ce0527eb859f47ab6989f6b947159373 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -340,7 +340,7 @@ public class GamsLocationOptimiser {
 		addScalar(inDB, "maxLandExpansionRate", maxExpansion, 3);
 		
 		addScalar(inDB, "carbonPrice", countryInput.getCarbonPrice().getExportPrice(), 3);
-		addScalar(inDB, "woodPrice", ModelConfig.WOOD_PRICE, 3);
+		addScalar(inDB, "woodPrice", countryInput.getTimberPrice().getExportPrice(), 3);
 		
 		// Carbon fluxes
 		GAMSParameter cFluxRateP = inDB.addParameter("carbonFluxRate", 3);
@@ -365,17 +365,29 @@ public class GamsLocationOptimiser {
 		}
 		
 		// Wood yields
+		Map<Integer, ? extends WoodYieldItem> woodYieldData = inputData.getWoodYields();
+		Map<Integer, Double> naturalWoodYieldFactor = inputData.getNaturalWoodYieldFactor();
+		
+		// Scale yields for natural areas
+		for (Entry<Integer, Double> entry : naturalWoodYieldFactor.entrySet()) {
+			int locId = entry.getKey();
+			double factor = entry.getValue();
+			for (LandCoverType toLC : LandCoverType.getConvertibleTypes()) {
+				WoodYieldItem item = woodYieldData.get(locId);
+				double maxYield = item.getWoodYield(LandCoverType.NATURAL, toLC);
+				item.setWoodYield(LandCoverType.NATURAL, toLC, maxYield * factor);
+			}
+		}
+		
 		GAMSParameter woodYieldP = inDB.addParameter("woodYield", 3);
 		
-		for (Entry<Integer, ? extends WoodYieldItem> entry : inputData.getWoodYields().entrySet()) {
+		for (Entry<Integer, ? extends WoodYieldItem> entry : woodYieldData.entrySet()) {
 			Integer locationId = entry.getKey();
 			String locString = Integer.toString(locationId);
 			WoodYieldItem wYield = entry.getValue();
 
 			for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
 				for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
-					if (newLC.equals(LandCoverType.TIMBER_FOREST))
-						continue;
 					Vector<String> v = new Vector<String>();
 					v.add(prevLC.getName());
 					v.add(newLC.getName());
@@ -629,8 +641,11 @@ public class GamsLocationOptimiser {
 		
 		// Carbon flux
 		double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue();
+		
+		// Timber harvest
+		double totalTimberProduction = outDB.getParameter("totalTimberProduction").getFirstRecord().getValue();
 
-		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, newForest, reservedDeforested, netCarbonFlux);
+		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, newForest, reservedDeforested, netCarbonFlux, totalTimberProduction);
 		return results;
 	}
 	
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 5cd4ad93cfc54a3c2411d25bae172125218e9aba..f0b45e884670bf8654c84ab786a4b366eca36109 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -18,10 +18,11 @@ public class GamsLocationOutput {
 	
 	Map<Integer, LandUseItem> landUses;  // data mapped from id (not raster)
 	private Map<CropType, CropUsageData> cropUsageData;
-	TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
-	List<ForestRecord> newForest;
-	DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
-	double netCarbonFlux;
+	private TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
+	private List<ForestRecord> newForest;
+	private DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
+	private double netCarbonFlux;
+	private double timberProduction;
 	
 	public GamsLocationOutput(ModelStat status, 
 			Map<Integer, LandUseItem> landUses, 
@@ -29,7 +30,7 @@ public class GamsLocationOutput {
 			TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange,
 			List<ForestRecord> newForest,
 			DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
-			double netCarbonFlux) {
+			double netCarbonFlux, double timberProduction) {
 		super();
 		this.status = status;
 		this.landUses = landUses;
@@ -38,6 +39,7 @@ public class GamsLocationOutput {
 		this.newForest = newForest;
 		this.reservedDeforested = reservedDeforested;
 		this.netCarbonFlux = netCarbonFlux;
+		this.timberProduction = timberProduction;
 	}
 	
 	public ModelStat getStatus() {
@@ -66,4 +68,8 @@ public class GamsLocationOutput {
 	public double getNetCarbonFlux() {
 		return netCarbonFlux;
 	}
+	
+	public double getTimberProduction() {
+		return timberProduction;
+	}
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterInput.java b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
index cb06c7644319eef6e1690606e93ce82825e16673..f14c6c8ded6fe48a9e5f037c042e96dbaa9a2bd7 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterInput.java
@@ -1,5 +1,7 @@
 package ac.ed.lurg.country.gams;
 
+import java.util.Map;
+
 import ac.ed.lurg.Timestep;
 import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.landuse.WoodYieldItem;
@@ -20,12 +22,14 @@ public class GamsRasterInput {
 	private RasterSet<WoodYieldItem> woodYields;
 	private DoubleMap<Integer, LandCoverType, Double> minimumLandCover;
 	private DoubleMap<Integer, LandCoverType, Double> timberYield;
+	private Map<Integer, Double> naturalWoodYieldFactor;
 	private DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts;
 	private GamsCountryInput countryInput;
 
 	public GamsRasterInput(Timestep timestep, YieldRaster yields, RasterSet<LandUseItem> previousLandsUses, RasterSet<IrrigationItem> irrigationCost, 
 			RasterSet<CarbonFluxItem> carbonFluxes, RasterSet<WoodYieldItem> woodYields, DoubleMap<Integer, LandCoverType, Double> minimumLandCover,
-			DoubleMap<Integer, LandCoverType, Double> timberYield, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts, GamsCountryInput countryInput) {
+			DoubleMap<Integer, LandCoverType, Double> timberYield, Map<Integer, Double> naturalWoodYieldFactor, 
+			DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts, GamsCountryInput countryInput) {
 		super();
 		this.timestep = timestep;
 		this.yields = yields;
@@ -35,6 +39,7 @@ public class GamsRasterInput {
 		this.woodYields = woodYields;
 		this.minimumLandCover = minimumLandCover;
 		this.timberYield = timberYield;
+		this.naturalWoodYieldFactor = naturalWoodYieldFactor;
 		this.conversionCosts = conversionCosts;
 		this.countryInput = countryInput;
 	}
@@ -71,6 +76,10 @@ public class GamsRasterInput {
 		return timberYield;
 	}
 	
+	public Map<Integer, Double> getNaturalWoodYieldFactor() {
+		return naturalWoodYieldFactor;
+	}
+
 	public GamsCountryInput getCountryInput() {
 		return countryInput;
 	}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index fb99a88c3bea8f08212f01a2d4c255cc65bb89eb..2eb23fe87e4cc89845da9affcc20be061219a70b 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -53,7 +53,8 @@ public class GamsRasterOptimiser {
 	private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {		
 		RasterSet<LandUseItem> newIntensityRaster = allocAreas(gamsInput.getPreviousLandUse(), gamsOutput, gamsInput.getTimestep().getYear());
 
-		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getNewForest(), gamsOutput.getReservedDeforested(), gamsOutput.netCarbonFlux);
+		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getNewForest(), 
+				gamsOutput.getReservedDeforested(), gamsOutput.getNetCarbonFlux(), gamsOutput.getTimberProduction());
 	}
 
 	private RasterSet<LandUseItem> createWithSameLandCovers(RasterSet<LandUseItem> toCopy) {
@@ -293,27 +294,44 @@ public class GamsRasterOptimiser {
 				}
 				*/
 				
+				/*
+				if (clusterId == 16) {
+					int foo =1;
+				} else {
+					continue;
+				}
+				*/
+				
 				// Aggregate carbon fluxes and wood yields
 				for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
 					for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
 						// Carbon flux
+						double carbonFlux;
 						if (carbonFluxItem != null) {
-							if (!aggCFlux.checkForKeys(prevLC, newLC)) { // if not seen yet
-								aggCFlux.setCarbonFlux(prevLC, newLC, carbonFluxItem.getCarbonFlux(prevLC, newLC));
-							} else {
-								aggCFlux.setCarbonFlux(prevLC, newLC, aggregateMean(aggCFlux.getCarbonFlux(prevLC, newLC), suitableAreaSoFar, 
-										carbonFluxItem.getCarbonFlux(prevLC, newLC), suitableAreaThisTime));
-							}
+							carbonFlux = carbonFluxItem.getCarbonFlux(prevLC, newLC);
+						} else {
+							carbonFlux = 0.0; // if missing data assume 0 TODO is this right?
 						}
 
+						if (!aggCFlux.checkForKeys(prevLC, newLC)) { // if not seen yet
+							aggCFlux.setCarbonFlux(prevLC, newLC, carbonFlux);
+						} else {
+							aggCFlux.setCarbonFlux(prevLC, newLC, aggregateMean(aggCFlux.getCarbonFlux(prevLC, newLC), suitableAreaSoFar, 
+									carbonFlux, suitableAreaThisTime));
+						}
 						// Wood yield
+						double woodYield;
 						if (woodYieldItem != null) {
-							if (!aggWYield.checkForKeys(prevLC, newLC)) { // if not seen yet
-								aggWYield.setWoodYield(prevLC, newLC, woodYieldItem.getWoodYield(prevLC, newLC));
-							} else {
-								aggWYield.setWoodYield(prevLC, newLC, aggregateMean(aggWYield.getWoodYield(prevLC, newLC), suitableAreaSoFar, 
-										woodYieldItem.getWoodYield(prevLC, newLC), suitableAreaThisTime));
-							}
+							woodYield = woodYieldItem.getWoodYield(prevLC, newLC);
+						} else {
+							woodYield = 0.0; // if missing data assume 0
+						}
+
+						if (!aggWYield.checkForKeys(prevLC, newLC)) { // if not seen yet
+							aggWYield.setWoodYield(prevLC, newLC, woodYield);
+						} else {
+							aggWYield.setWoodYield(prevLC, newLC, aggregateMean(aggWYield.getWoodYield(prevLC, newLC), suitableAreaSoFar, 
+									woodYield, suitableAreaThisTime));
 						}
 					}
 				}
@@ -395,7 +413,8 @@ public class GamsRasterOptimiser {
 		checkedTotalAreas(aggregatedAreas, "after");
 		
 		return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, 
-				aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getMinimumLandCover(), rasterInputData.getTimberYield(), rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
+				aggregatedCarbonFluxes, aggregatedWoodYields, rasterInputData.getMinimumLandCover(), rasterInputData.getTimberYield(), rasterInputData.getNaturalWoodYieldFactor(),
+				rasterInputData.getConversionCosts(), rasterInputData.getCountryInput());
 	}
 
 	private void logWarningWithCoord(String message, RasterKey key, YieldRaster yieldRaster, CropType crop) {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index 110bb241e17f1f03c41dc65da5cee986f306263f..c6875491caa2ee6e1f779ccac611223ffa934eba 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -21,6 +21,7 @@ public class GamsRasterOutput {
 	private List<ForestRecord> newForest;
 	DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
 	double netCarbonFlux;
+	double timberProduction;
 
 	public GamsRasterOutput(RasterSet<LandUseItem> landUses, Map<CropType, CropUsageData> cropUsageData) {
 		super();
@@ -30,12 +31,13 @@ public class GamsRasterOutput {
 
 	public GamsRasterOutput(ModelStat status, RasterSet<LandUseItem> intensityRaster, Map<CropType, CropUsageData> cropUsageData,
 			List<ForestRecord> newForest, DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
-			double netCarbonFlux) {
+			double netCarbonFlux, double timberProduction) {
 		this(intensityRaster, cropUsageData);
 		this.status = status;
 		this.newForest = newForest;
 		this.reservedDeforested = reservedDeforested;
 		this.netCarbonFlux = netCarbonFlux;
+		this.timberProduction = timberProduction;
 	}
 	
 	public ModelStat getStatus() {
@@ -61,4 +63,9 @@ public class GamsRasterOutput {
 	public double getNetCarbonFlux() {
 		return netCarbonFlux;
 	}
+
+	public double getTimberProduction() {
+		return timberProduction;
+	}
+	
 }
diff --git a/src/ac/ed/lurg/demand/AbstractDemandManager.java b/src/ac/ed/lurg/demand/AbstractDemandManager.java
index f97074d441eb7ee5be93fa235b2f9e1d1b7e3fcd..e981c1fedbe9dd8ebed766ef5cb355b233deb13d 100644
--- a/src/ac/ed/lurg/demand/AbstractDemandManager.java
+++ b/src/ac/ed/lurg/demand/AbstractDemandManager.java
@@ -20,6 +20,7 @@ public abstract class AbstractDemandManager {
 	protected BioenergyDemandManager bioenergyDemandManager;
 	protected CerealFractionsManager cerealFractionsManager;
 	protected CarbonDemandManager carbonDemandManager;
+	protected TimberDemandManager timberDemandManager;
 
 	public AbstractDemandManager(CompositeCountryManager compositeCountryManager, CalorieManager calorieManager) {
 		this.compositeCountryManager = compositeCountryManager;
@@ -27,6 +28,7 @@ public abstract class AbstractDemandManager {
 		bioenergyDemandManager = new BioenergyDemandManager();
 		cerealFractionsManager = new CerealFractionsManager(compositeCountryManager);
 		carbonDemandManager = new CarbonDemandManager();
+		timberDemandManager = new TimberDemandManager();
 	}
 
 	public Map<CommodityType, Double> getDemand(CompositeCountry cc, int year, Map<CommodityType, Double> prices, boolean outputGamsDemand) {
@@ -136,5 +138,9 @@ public abstract class AbstractDemandManager {
 	public double getCarbonDemand(Timestep timestep) {
 		return carbonDemandManager.getGlobalCarbonDemand(timestep.getYear());
 	}
+	
+	public double getTimberDemand(Timestep timestep) {
+		return timberDemandManager.getGlobalTimberDemand(timestep.getYear());
+	}
 
 }
diff --git a/src/ac/ed/lurg/demand/TimberDemandManager.java b/src/ac/ed/lurg/demand/TimberDemandManager.java
new file mode 100644
index 0000000000000000000000000000000000000000..560ea3bc94548a78cf9774c33a67738fc4c363e9
--- /dev/null
+++ b/src/ac/ed/lurg/demand/TimberDemandManager.java
@@ -0,0 +1,64 @@
+package ac.ed.lurg.demand;
+
+import java.io.BufferedReader;
+import java.io.FileReader;
+import java.io.IOException;
+import java.util.HashMap;
+import java.util.Map;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.utils.Interpolator;
+import ac.ed.lurg.utils.LogWriter;
+
+public class TimberDemandManager {
+	private Map<Integer, Double> globalTimberDemand; // global demand for timber
+	private static final int YEAR_COL = 0;
+	private static final int DEMAND_COL = 1; 
+	
+	public TimberDemandManager() {
+		readTimberDemandData();
+	}
+	
+	public void readTimberDemandData() {
+		
+		Map<Integer, Double> data = new HashMap<Integer, Double>();
+		
+		String filename = ModelConfig.TIMBER_DEMAND_FILE;
+		try {
+			BufferedReader reader = new BufferedReader(new FileReader(filename)); 
+			String line;
+			Integer year;
+			Double demand;
+			reader.readLine(); // read header
+
+			while ((line=reader.readLine()) != null) {
+				String[] tokens = line.split(",");
+				
+				if (tokens.length < 2)
+					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
+				
+				year = Integer.valueOf(tokens[YEAR_COL]);
+				demand = Double.valueOf(tokens[DEMAND_COL]);
+				
+				data.put(year,  demand);
+			} 
+			reader.close(); 
+		
+		} catch (IOException e) {
+			LogWriter.printlnError("Failed in reading timber demand data");
+			LogWriter.print(e);
+		}
+		globalTimberDemand = data;
+		LogWriter.println("Processed " + filename);
+	}
+	
+	public double getGlobalTimberDemand(int year) {
+		int downYear = (year/5) * 5;
+		int upYear = downYear + 5;
+		Double lowerD = globalTimberDemand.get(downYear);
+		Double upperD = globalTimberDemand.get(upYear);
+		double factor = ((double)(year - downYear)) / (upYear - downYear);
+		Double d = Interpolator.interpolate(lowerD, upperD, factor);
+		return d;
+	}
+}