diff --git a/.gitignore b/.gitignore
index 41bf33c0caa240badd5fb8b349a0db294277f8d3..9b98b0db951b82ff8942f9619ce890547879343d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,5 @@
 /bin/
 plumv2/output/
+/.classpath
+/build.xml
+/cluster.asc
diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 8a750b8c6da1c604feb509ca200c3714cd98260b..6bdee3ba14e6ad83c055a94cd11d5e2f2d9daff9 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -1,5 +1,5 @@
 
- SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, wood, carbonCredits/;
+ SET all_types / monogastrics, ruminants, cereals, oilcropspulses, wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside, roundwood, fuelwood, carbonCredits/;
 
  SET crop(all_types)                                             / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops, pasture, setaside /;
  SET crop_less_pasture(crop)                                     / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops,          setaside/;
@@ -8,7 +8,7 @@
  SET oilpulse_crop(crop)                                                             / oilcrops, pulses /;
  SET feed_crop(crop)                                             / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg,                     pasture/;
  SET feed_crop_less_pasture(feed_crop)                           / wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg /;
- SET import_types(all_types)   / monogastrics, ruminants,          wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops,                    wood, carbonCredits/;
+ SET import_types(all_types)   / monogastrics, ruminants,          wheat, maize, rice, oilcrops, pulses, starchyRoots, fruitveg, sugar, energycrops,                    roundwood, fuelwood, carbonCredits/;
  SET animal(all_types)     / monogastrics, ruminants /;
 
  SET land_cover / cropland, pasture, timberForest, carbonForest, natural /;
@@ -17,6 +17,8 @@
  SET wood_producing(land_cover) / timberForest, natural /;
  ALIAS (land_cover, land_cover_before);
  ALIAS (land_cover, land_cover_after);
+ 
+ SET wood_type(all_types) / roundwood, fuelwood /;
 
  SET location;
  SET year;
@@ -52,10 +54,10 @@
  
  PARAMETER previousLandCoverArea(land_cover, location)                      land cover area in Mha;
  PARAMETER carbonCreditRate(land_cover_before, land_cover_after, location) potential carbon credits from LULUCF - tC-eq per ha;
- 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);
+ PARAMETER maxRotationIntensity(wood_type);
 
  SCALAR meatEfficency                        efficiency of converting feed and pasture into animal products;
  SCALAR fertiliserUnitCost                   fert cost at max fert rate;
@@ -80,7 +82,7 @@ $load previousCropArea, previousFertIntensity, previousIrrigIntensity, previousO
 $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, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam
+$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMax, woodYieldParam, maxRotationIntensity
 $load forestManagementCost, vegClearingCostRate, carbonForestMaxProportion, maxFertChange, maxIrrigChange
 $gdxin
 
@@ -146,15 +148,16 @@ $gdxin
        landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
        totalConversionCost(location)      land cover conversion cost - $1000 per ha
 
-       woodYieldRota(location)
-       rotationPeriod(location)
+       timberForestArea(wood_type, location)
+       woodYieldRota(wood_type, location)
+       rotationIntensity(wood_type, location)
        forestryCost(location)
        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, animalProd, woodYieldRota, rotationPeriod;
+                   landCoverArea, landCoverChange, totalConversionCost, animalProd, woodYieldRota, rotationIntensity, timberForestArea;
 
 
 * POSITIVE VARIABLE A;
@@ -193,10 +196,11 @@ $gdxin
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
        CONVERSION_COST_CALC(location)                          cost of land cover conversion
 
-       WOOD_YIELD_CALC(location)
-       WOOD_DEMAND_CONSTRAINT
-       ROTATION_PERIOD_MIN_CONSTRAINT(location)
-       ROTATION_PERIOD_MAX_CONSTRAINT(location)
+       WOOD_YIELD_CALC(wood_type, location)
+       WOOD_DEMAND_CONSTRAINT(wood_type)
+       ROTA_INTENSITY_MAX_CONSTRAINT(wood_type, location)
+       ROTA_INTENSITY_MIN_CONSTRAINT(wood_type, location)
+       TIMBER_FOREST_AREA_CONSTRAINT(location)
        FORESTRY_COST_CALC(location)
  
        CARBON_CREDIT_CALC(location)
@@ -281,15 +285,18 @@ $gdxin
  
 ************* Forestry ***********************************
  
- WOOD_YIELD_CALC(location) .. woodYieldRota(location) =E= woodYieldMax(location) * (1 - exp(-woodYieldParam(location) * rotationPeriod(location))) / rotationPeriod(location);
+ WOOD_YIELD_CALC(wood_type, location) .. woodYieldRota(wood_type, location) =E= woodYieldMax(location) * (1 - exp(-woodYieldParam(location) / rotationIntensity(wood_type, location))) * rotationIntensity(wood_type, location);
+ 
+ ROTA_INTENSITY_MAX_CONSTRAINT(wood_type, location) .. rotationIntensity(wood_type, location) =L= maxRotationIntensity(wood_type);
+ 
+ ROTA_INTENSITY_MIN_CONSTRAINT(wood_type, location) .. rotationIntensity(wood_type, location) =G= 1 / 160;
  
- ROTATION_PERIOD_MIN_CONSTRAINT(location) .. rotationPeriod(location) =G= 5;
- ROTATION_PERIOD_MAX_CONSTRAINT(location) .. rotationPeriod(location) =L= 160;
+ TIMBER_FOREST_AREA_CONSTRAINT(location) .. sum(wood_type, timberForestArea(wood_type, location)) =E= landCoverArea('timberForest', location);
  
- WOOD_DEMAND_CONSTRAINT .. demand('wood') + exportAmount('wood') - importAmount('wood') =L= sum(location, woodYieldRota(location) * landCoverArea('timberForest', location));
+ WOOD_DEMAND_CONSTRAINT(wood_type) .. demand(wood_type) + exportAmount(wood_type) - importAmount(wood_type) =L= sum(location, woodYieldRota(wood_type, location) * timberForestArea(wood_type, location));
  
  FORESTRY_COST_CALC(location) .. forestryCost(location) =E= sum(managed_forest, landCoverArea(managed_forest, location) * 0.05) +
-                                                            forestManagementCost / rotationPeriod(location);
+                                                            sum(wood_type, forestManagementCost * rotationIntensity(wood_type, location));
  
 *********** Carbon ***********************************
                                        
@@ -308,8 +315,8 @@ $gdxin
 
 ************ Total cost ******************************
 
- COST_CALC .. total_cost =E= sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) * domesticPriceMarkup +
-                             sum(animal, animalProd(animal) * animalProdCost(animal) * (1 - subsidyRate(animal))) +
+ COST_CALC .. total_cost =E= sum((crop, location), cropArea(crop, location) * unitCost(crop, location) - subsidyRate(crop)) * domesticPriceMarkup +
+                             sum(animal, animalProd(animal) * animalProdCost(animal) - subsidyRate(animal)) +
                              sum(location, totalConversionCost(location)) +
                              sum(location, forestryCost(location)) +
                              sum(import_types, importPrices(import_types) * importAmount(import_types)) -
@@ -326,7 +333,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;
+ rotationIntensity.L(wood_type, location) = 0.02;
 
 * LAND_USE.OptFile = 1;
 
@@ -351,7 +358,7 @@ $gdxin
  parameter productionShock(all_types);
  parameter carbonFlux(location);
  scalar netCarbonCredits;
- scalar woodSupplyRota;
+ parameter woodSupplyRota(wood_type);
 
 * Production quantities based on smaller area (before unhandledCropArea adjustment applied)
  totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
@@ -374,7 +381,7 @@ $gdxin
  
  netCarbonCredits = SUM(location, carbonCredits.L(location));
  
- woodSupplyRota = sum(location, woodYieldRota.L(location) * landCoverArea.L('timberForest', location));
+ woodSupplyRota(wood_type) = sum(location, woodYieldRota.L(wood_type, location) * timberForestArea.L(wood_type, location));
 
  Scalar totalCostsLessLU;
 
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 790fec5c2967dccb33fa36e97003b9e7780dd7e5..cb8be1b44ad423d1cdd29e530e6c7091abc62a2e 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -30,7 +30,7 @@ public class InternationalMarket {
 
 	private Map<CropType, GlobalPrice> worldPrices;
 	private GlobalPrice carbonPrice;
-	private GlobalPrice woodPrice;
+	private Map<WoodType, GlobalPrice> woodPrices;
 	private PriceShockManager priceShockManager;
 	private PriceCapManager priceCapManager;
 	
@@ -48,14 +48,22 @@ public class InternationalMarket {
 			}
 			
 			carbonPrice = GlobalPrice.createInitial(ModelConfig.INIT_CARBON_PRICE, ModelConfig.INIT_CARBON_STOCK);
-			woodPrice = GlobalPrice.createInitial(ModelConfig.INIT_WOOD_PRICE, ModelConfig.INIT_WOOD_STOCK);		
+
+			woodPrices = new HashMap<>();
+			for (WoodType woodType : WoodType.values()) {
+				if (ModelConfig.IS_FORESTRY_ON) {
+					woodPrices.put(woodType, GlobalPrice.createInitial(woodType.getInitialPrice(), ModelConfig.INIT_WOOD_STOCK));
+				} else {
+					woodPrices.put(woodType, GlobalPrice.createInitial(0, 0));
+				}
+			}
 		}
 		else {
 			List<Object> deserializedPrices = deserializeGlobalPrice();
 			worldPrices = (Map<CropType, GlobalPrice>) deserializedPrices.get(0);
-			//carbonPrice = (GlobalPrice) deserializedPrices.get(1);
+			//carbonPrice = (GlobalPrice) deserializedPrices.get(1); TODO
 			carbonPrice = GlobalPrice.createInitial(ModelConfig.INIT_CARBON_PRICE, ModelConfig.INIT_CARBON_STOCK);
-			woodPrice = (GlobalPrice) deserializedPrices.get(2);
+			woodPrices = (Map<WoodType, GlobalPrice>) deserializedPrices.get(2);
 						
 			if (ModelConfig.RESET_ENERGYCROP_PRICE) {
 				double previousExportPrice = worldPrices.get(CropType.ENERGY_CROPS).getExportPrice();
@@ -80,12 +88,8 @@ public class InternationalMarket {
 		}
 	}
 	
-	public GlobalPrice getWoodPrice() {
-		if (ModelConfig.IS_FORESTRY_ON) {
-			return woodPrice;
-		} else {
-			return woodPrice.createPriceAdjustedByFactor(0); 
-		}
+	public Map<WoodType, GlobalPrice> getWoodPrices() {
+		return woodPrices;
 	}
 
 	private Map<CropType, Double> getInitialStockLevels() {
@@ -100,7 +104,11 @@ public class InternationalMarket {
 			if (initialStock != null)
 				entry.getValue().resetStock(initialStock);
 		}
-		woodPrice.resetStock(ModelConfig.INIT_WOOD_STOCK);
+
+		for (Map.Entry<WoodType, GlobalPrice> entry : woodPrices.entrySet()) {
+			entry.getValue().resetStock(ModelConfig.INIT_WOOD_STOCK);
+		}
+
 		carbonPrice.resetStock(0.0);
 	}
 
@@ -166,31 +174,32 @@ public class InternationalMarket {
 		carbonPrice = adjustedCPrice;
 		
 		// Update timber price
-		double totalWoodImport = 0;
-		double totalWoodExport = 0;
-		double totalWoodProduction = 0;
-		
-		for (AbstractCountryAgent ca : countryAgents) {
-			Map<WoodType, WoodUsageData> woodUsageMap = ca.getWoodUsageData();
-			for (WoodUsageData woodUsage : woodUsageMap.values()) {
+		for (WoodType woodType : WoodType.values()) {
+			double totalWoodImport = 0;
+			double totalWoodExport = 0;
+			double totalWoodProduction = 0;
+
+			for (AbstractCountryAgent ca : countryAgents) {
+				WoodUsageData woodUsage = ca.getWoodUsageData().get(woodType);
 				totalWoodProduction += woodUsage.getHarvest();
 				double netImport = woodUsage.getNetImport();
-				double lucHarvest = woodUsage.getLucHarvest(); // assume wood from LUC exported
 				if (netImport >= 0) {
 					totalWoodImport += netImport;
 				} else {
-					totalWoodExport += -netImport + lucHarvest;
+					totalWoodExport += -netImport;
 				}
 			}
+			totalWoodProduction = Math.max(totalWoodProduction, 0.0000001); // avoid division by 0
+			GlobalPrice prevTPrice = woodPrices.get(woodType);
+			LogWriter.println(timestep.getYear() + " Updating " + woodType.getName() + " price", 2);
+			GlobalPrice adjustedTPrice = prevTPrice.createWithUpdatedMarketPrices(totalWoodImport, totalWoodExport,
+					timestep, totalWoodProduction, true);
+			LogWriter.println( String.format("Price for wood updated from %s to %s \n", prevTPrice, adjustedTPrice), 2);
+			if (adjustedTPrice.getStockLevel() < 0)
+				LogWriter.println("Global stocks are below zero wood, " + timestep.getYear(), 2);
+			woodPrices.put(woodType, adjustedTPrice);
 		}
-		totalWoodProduction = Math.max(totalWoodProduction, 0.0000001); // avoid division by 0
-		GlobalPrice prevTPrice = woodPrice;
-		LogWriter.println(timestep.getYear() + " Updating wood price", 2);
-		GlobalPrice adjustedTPrice = prevTPrice.createWithUpdatedMarketPrices(totalWoodImport, totalWoodExport, timestep, totalWoodProduction, true);
-		LogWriter.println( String.format("Price for wood updated from %s to %s \n", prevTPrice, adjustedTPrice), 2);
-		if (adjustedTPrice.getStockLevel() < 0)
-			LogWriter.println("Global stocks are below zero wood, " + timestep.getYear(), 2);
-		woodPrice = adjustedTPrice;
+
 		
 		// Cap prices
 		capPrices();
@@ -219,11 +228,12 @@ public class InternationalMarket {
 			outputFile.newLine();
 		}
 		// Timber price
-		{
+		for (WoodType woodType : WoodType.values()){
 			StringBuffer sbData = new StringBuffer();
-			sbData.append(String.format("%d,%s", timestep.getYear(), "timber"));
-			sbData.append(String.format(",%.1f,%.1f", woodPrice.getImportAmount(), woodPrice.getExportsAfterTransportLosses()));
-			sbData.append(String.format(",%.8f,%.3f", woodPrice.getExportPrice(), woodPrice.getStockLevel()));
+			GlobalPrice price = woodPrices.get(woodType);
+			sbData.append(String.format("%d,%s", timestep.getYear(), woodType.getName()));
+			sbData.append(String.format(",%.1f,%.1f", price.getImportAmount(), price.getExportsAfterTransportLosses()));
+			sbData.append(String.format(",%.8f,%.3f", price.getExportPrice(), price.getStockLevel()));
 	
 			outputFile.write(sbData.toString());
 			outputFile.newLine();
@@ -234,7 +244,7 @@ public class InternationalMarket {
 		List<Object> pricesToSerialize = new ArrayList<Object>();
 		pricesToSerialize.add(worldPrices);
 		pricesToSerialize.add(carbonPrice);
-		pricesToSerialize.add(woodPrice);
+		pricesToSerialize.add(woodPrices);
 
 		try {
 			String fileStr = ModelConfig.IS_CALIBRATION_RUN ? ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE : ModelConfig.CHECKPOINT_INTERNATIONAL_MARKET_FILE;
@@ -300,8 +310,8 @@ public class InternationalMarket {
 		Map<CropType, GlobalPrice> cappedCropPrices = priceCapManager.capCropPrices(worldPrices);
 		worldPrices.putAll(cappedCropPrices);
 		
-		GlobalPrice updatedWoodPrice = priceCapManager.capWoodPrices(woodPrice);
-		woodPrice = updatedWoodPrice;
+		Map<WoodType, GlobalPrice> updatedWoodPrice = priceCapManager.capWoodPrices(woodPrices);
+		woodPrices = updatedWoodPrice;
 		
 		GlobalPrice updatedCarbonPrice = priceCapManager.capCarbonPrices(carbonPrice);
 		carbonPrice = updatedCarbonPrice;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index dec4e05755a17423bac861da93cf5d26ad68bca9..836228bf6da2da40cf29a73aa719754ad0aad16b 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -147,7 +147,7 @@ public class ModelConfig {
 	public static final String BASE_DIR = getProperty("BASE_DIR");  // this must to be set in config file
     public static final String OUTPUT_DIR = getProperty("OUTPUT_DIR", ".");
 	public static final String TEMP_DIR = getProperty("TEMP_DIR", OUTPUT_DIR + File.separator + "GamsTmp");
-	public static final String DATA_DIR = getProperty("DATA_DIR", BASE_DIR + File.separator + "data");
+	public static final String DATA_DIR = getProperty("DATA_DIR", BASE_DIR + File.separator + "data_update"); // TODO: change
 	public static final String GAMS_DIR = getProperty("GAMS_DIR", BASE_DIR + File.separator + "GAMS");
 	public static final boolean CLEANUP_GAMS_DIR = getBooleanProperty("CLEANUP_GAMS_DIR", true);
 
@@ -172,7 +172,7 @@ public class ModelConfig {
 	public static final String SSP_FILE = getProperty("SSP_FILE", DATA_DIR + File.separator + SSP_FILENAME);	
 	public static final String BASELINE_CONSUMP_FILE = DATA_DIR + File.separator + "base_consump.csv";
 	public static final String CALORIE_PER_T_FILE = DATA_DIR + File.separator + "calories_per_t.csv";
-	public static final String COUNTRY_CODES_FILE = DATA_DIR + File.separator + "country_codes4.csv";
+	public static final String COUNTRIES_FILE = DATA_DIR + File.separator + "countries.csv";
 	public static final String COUNTRY_DATA_FILE = DATA_DIR + File.separator + "country_data.csv";
 	public static final String NET_IMPORTS_FILE = DATA_DIR + File.separator + "net_imports.csv";
 	public static final String BIOENERGY_1GEN_BASE_DEMAND_FILE = DATA_DIR + File.separator + "bio_demand.csv";
@@ -193,7 +193,7 @@ public class ModelConfig {
 	public static final String SUBSIDY_RATE_FILE = getProperty("SUBSIDY_RATE_FILE", DATA_DIR + File.separator + SUBSIDY_RATE_FILENAME);
 	public static final String ANIMAL_RATES_FILE = DATA_DIR + File.separator + "animal_numbers.csv";;
 	public static final String INITIAL_CONSUMER_PRICE_FILE = DATA_DIR + File.separator + "consumerprices.csv";;
-	public static final String GDP_FRACTIONS_FILE = DATA_DIR + File.separator + "agriculturalGdpFraction.csv";
+	public static final String GDP_FRACTIONS_FILE = DATA_DIR + File.separator + "agri_gdp_fraction.csv";
 
 	// yield data
 	public static final String YIELD_DIR_BASE = getProperty("YIELD_DIR_BASE");
@@ -299,7 +299,7 @@ public class ModelConfig {
 	public static final String LAND_COVER_AGE_DIST_FILENAME = SPATIAL_DATA_DIR + File.separator + "land_cover_age_dist.txt";
 	public static final int LAND_COVER_INIT_AGE_GROUP_SIZE = getIntProperty("LAND_COVER_INIT_AGE_GROUP_SIZE", 10); // years
 	public static final int CARBON_DATA_TIMESTEP_SIZE = getIntProperty("WOOD_AND_CARBON_TIMESTEP_SIZE", 20); // years
-	public static final double WOOD_YIELD_CALIB_FACTOR = getDoubleProperty("WOOD_YIELD_CALIB_FACTOR", 1.08); 
+	public static final double WOOD_YIELD_CALIB_FACTOR = getDoubleProperty("WOOD_YIELD_CALIB_FACTOR", 1.0); 
 	public static final int CARBON_DATA_MIN_YEAR = getIntProperty("CARBON_DATA_MIN_YEAR", 1850); // years
 
 	// Output
@@ -378,7 +378,7 @@ public class ModelConfig {
 	public static final double DIETARY_CLOSURE = getDoubleProperty("DIETARY_CLOSURE", 0.25); // Amount diet converges in DIETARY_CLOSURE_GDP_CHANGE rate of GDP shift
 	public static final double DIETARY_CLOSURE_GDP_CHANGE = getDoubleProperty("DIETARY_CLOSURE_GDP_CHANGE", 2.0); // 2 is double of GDP
 	public static final double DIETARY_CLOSURE_PARAM = getDoubleProperty("DIETARY_CLOSURE_PARAM", Math.log((1-DIETARY_CLOSURE))/DIETARY_CLOSURE_GDP_CHANGE); // Zero is no dietary closure, number specifies closure rate for changes in GDP per capita, e.g. Math.log(0.5)/2
-	public static final String SSP_SCENARIO = getProperty("SSP_SCENARIO", "SSP1_v9_130325");
+	public static final String SSP_SCENARIO = getProperty("SSP_SCENARIO", "SSP2");
 	public static final ModelFitType DEMAND_ANIMAL_PROD_FIT = ModelFitType.findByName(getProperty("DEMAND_ANIMAL_PROD_FIT", "loglinear"));
 	public static final ModelFitType DEMAND_NON_ANIMAL_PROD_FIT = ModelFitType.findByName(getProperty("DEMAND_NON_ANIMAL_PROD_FIT", "loglinear"));
 	public static final boolean LIMIT_DEMAND_FRACTION = getBooleanProperty("LIMIT_DEMAND_FRACTION", true);
@@ -473,7 +473,7 @@ public class ModelConfig {
 	public static final int LPJG_TIMESTEP_SIZE = getIntProperty("LPJG_TIMESTEP_SIZE", 5);
 	public static final int LPJ_YEAR_OFFSET = getIntProperty("LPJ_YEAR_OFFSET", 0);;
 	
-	public static final int NUM_YIELD_CLUSTERS = getIntProperty("NUM_YIELD_CLUSTERS", 8000);
+	public static final int NUM_YIELD_CLUSTERS = getIntProperty("NUM_YIELD_CLUSTERS", 1000);
 	public static final long RANDOM_SEED = getIntProperty("RANDOM_SEED", 1974329);  // any number will do
 	
 	public static final String CHECKPOINT_YEARS = getProperty("CHECKPOINT_YEARS", null); // 2020,2050
@@ -503,17 +503,18 @@ 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_adj.csv");
+	public static final String WOOD_DEMAND_FILENAME = getProperty("WOOD_DEMAND_FILENAME", "wood_demand.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
 	public static final double FUELWOOD_DEMAND_ELASTICITY = getDoubleProperty("FUELDWOOD_DEMAND_ELASTICITY", -0.3598551);
-	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 WOOD_BIOMASS_CONVERSION_FACTOR = getDoubleProperty("WOOD_BIOMASS_CONVERSION_FACTOR", 0.3); // m3 to tC-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", 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.1) : 0.0; // $1000/tC-eq
+	public static final double INIT_ROUNDWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_ROUNDWOOD_PRICE", 0.36) : 0.0; // $1000/tC-eq
+	public static final double INIT_FUELWOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_FUELWOOD_PRICE", 0.24) : 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
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index c529dc40c8171648b35050b8d05401f2db785c35..94d9f10837c99597fbdbb5ce847a5d59592cf571 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -1,61 +1,16 @@
 package ac.ed.lurg;
 
-import java.io.BufferedWriter;
-import java.io.File;
-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.HashMap;
-import java.util.Map;
-import java.util.Map.Entry;
-
 import ac.ed.lurg.carbon.CarbonDataOutputer;
+import ac.ed.lurg.carbon.CarbonFluxItem;
 import ac.ed.lurg.carbon.CarbonFluxRasterSet;
 import ac.ed.lurg.carbon.CarbonFluxReader;
-import ac.ed.lurg.country.AbstractCountryAgent;
-import ac.ed.lurg.country.AnimalRateManager;
-import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
-import ac.ed.lurg.country.CountryAgentManager;
-import ac.ed.lurg.country.CountryBoundaryRaster;
-import ac.ed.lurg.country.CountryBoundaryReader;
-import ac.ed.lurg.country.CountryPrice;
-import ac.ed.lurg.country.GlobalPrice;
-import ac.ed.lurg.demand.AbstractDemandManager;
-import ac.ed.lurg.demand.BaseConsumpManager;
-import ac.ed.lurg.demand.CalorieManager;
-import ac.ed.lurg.demand.DemandManagerFromFile;
-import ac.ed.lurg.demand.DemandManagerSSP;
-import ac.ed.lurg.demand.ElasticDemandManager;
+import ac.ed.lurg.country.*;
+import ac.ed.lurg.demand.*;
 import ac.ed.lurg.forestry.ForestryDataOutputer;
+import ac.ed.lurg.forestry.WoodYieldItem;
 import ac.ed.lurg.forestry.WoodYieldRasterSet;
-import ac.ed.lurg.landuse.CarbonUsageData;
-import ac.ed.lurg.landuse.ConversionCostReader;
-import ac.ed.lurg.landuse.CropUsageData;
-import ac.ed.lurg.landuse.CropUsageReader;
-import ac.ed.lurg.landuse.FPUManager;
-import ac.ed.lurg.landuse.InitProtectedAreasReader;
-import ac.ed.lurg.landuse.IrrigationConstraintReader;
-import ac.ed.lurg.landuse.IrrigationItem;
-import ac.ed.lurg.landuse.IrrigationMaxAmountReader;
-import ac.ed.lurg.landuse.IrrigationRasterSet;
-import ac.ed.lurg.landuse.IrrigiationCostReader;
-import ac.ed.lurg.landuse.LccKey;
-import ac.ed.lurg.landuse.LandCoverItem;
-import ac.ed.lurg.landuse.LandCoverReader;
-import ac.ed.lurg.landuse.LandCoverTile;
-import ac.ed.lurg.landuse.LandUseItem;
-import ac.ed.lurg.landuse.LandTileReader;
-import ac.ed.lurg.landuse.LandUseBinarySerializer;
-import ac.ed.lurg.landuse.MaxCropAreaReader;
-import ac.ed.lurg.landuse.ProtectedAreasReader;
-import ac.ed.lurg.landuse.RunOffReader;
-import ac.ed.lurg.landuse.WoodUsageData;
-import ac.ed.lurg.landuse.WoodUsageReader;
 import ac.ed.lurg.forestry.WoodYieldReader;
+import ac.ed.lurg.landuse.*;
 import ac.ed.lurg.output.LandUseOutputer;
 import ac.ed.lurg.output.LpjgOutputer;
 import ac.ed.lurg.types.CommodityType;
@@ -63,15 +18,16 @@ import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.WoodType;
 import ac.ed.lurg.utils.FileWriterHelper;
+import ac.ed.lurg.utils.KDTree;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.LPJYieldResponseMapReader;
 import ac.ed.lurg.yield.YieldRaster;
-import ac.sac.raster.IntegerRasterItem;
-import ac.sac.raster.IntegerRasterReader;
-import ac.sac.raster.RasterHeaderDetails;
-import ac.sac.raster.RasterKey;
-import ac.sac.raster.RasterOutputer;
-import ac.sac.raster.RasterSet;
+import ac.ed.lurg.yield.YieldResponsesItem;
+import ac.sac.raster.*;
+
+import java.io.*;
+import java.util.*;
+import java.util.Map.Entry;
 
 public class ModelMain {
 
@@ -79,10 +35,9 @@ public class ModelMain {
 	private CountryBoundaryRaster countryBoundaryRaster;
 	private AbstractDemandManager demandManager;
 	private AnimalRateManager animalRateManager;
-	private CompositeCountryManager compositeCountryManager;
-	LPJYieldResponseMapReader lpjYieldReader;
 	private RasterHeaderDetails desiredProjection;
-
+	private LPJYieldResponseMapReader lpjYieldReader;
+	private YieldRaster yieldSurfaces;
 	private InternationalMarket internationalMarket;
 	private IrrigationRasterSet currentIrrigationData;
 	private RasterSet<LandUseItem> globalLandUseRaster;
@@ -103,19 +58,19 @@ public class ModelMain {
 	/* setup models, reading inputs, etc. */
 	private void setup() {
 		desiredProjection = RasterHeaderDetails.getGlobalHeaderFromCellSize(ModelConfig.CELL_SIZE_X, ModelConfig.CELL_SIZE_Y, "999");
-
+		
 		BaseConsumpManager baseConsumpManager = new BaseConsumpManager();
 		CalorieManager calorieManager = new CalorieManager();
-		compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
+		
 		lpjYieldReader = new LPJYieldResponseMapReader(desiredProjection);
-		animalRateManager = new AnimalRateManager(compositeCountryManager);
+		animalRateManager = new AnimalRateManager();
 		
 		if (ModelConfig.DEMAND_FROM_FILE)
-			demandManager = new DemandManagerFromFile(compositeCountryManager,calorieManager);
+			demandManager = new DemandManagerFromFile(calorieManager);
 		else if (ModelConfig.PRICE_ELASTIC_DEMAND)
-			demandManager = new ElasticDemandManager(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager, compositeCountryManager);
+			demandManager = new ElasticDemandManager(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager);
 		else
-			demandManager = new DemandManagerSSP(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager, compositeCountryManager);
+			demandManager = new DemandManagerSSP(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager);
 			
 		currentIrrigationData = getFixedIrrigationData();
 		countryBoundaryRaster = getCountryBoundaryRaster();
@@ -127,7 +82,7 @@ public class ModelMain {
 		woodYieldReader = new WoodYieldReader(desiredProjection);
 		carbonFluxReader = new CarbonFluxReader(desiredProjection);
 		
-		createCountryAgents(compositeCountryManager.getAll());
+		createCountryAgents(CountryManager.getInstance().getAllCompositeCountries());
 	}
 
 	/* run the model */
@@ -145,9 +100,8 @@ public class ModelMain {
 	
 	private void doTimestep(Timestep timestep) {
 		LogWriter.println("Timestep: " + timestep.toString());
-		LogWriter.println("Memory usage 1: " + (Runtime.getRuntime().totalMemory() - Runtime.getRuntime().freeMemory()) / (1024.0*1024.0*1024.0));
 
-		YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
+		yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
 		getUpdateIrrigationData(timestep, yieldSurfaces); // updating currentIrrigationData
 
 		// When running half earth we can to alter protected areas data at a point in time
@@ -163,6 +117,8 @@ public class ModelMain {
 
 		ConversionCostReader conCostReader = new ConversionCostReader(timestep);
 		Map<LccKey, Double> conversionCosts = conCostReader.getConversionCosts();
+
+		handleMissingData();
 		
 		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, carbonFluxData, woodYieldData, conversionCosts);
 		
@@ -415,7 +371,7 @@ public class ModelMain {
 					double prod = woodUsage.getHarvest();
 					double netImports = woodUsage.getNetImport() - woodUsage.getLucHarvest(); // does not include transport losses
 
-					CountryPrice px = country.getCurrentCountryWoodPrice();
+					CountryPrice px = country.getCurrentCountryWoodPrices().get(woodType);
 					double importPrice = px.getImportPrice(); 
 					double exportPrice = px.getExportPrice();
 
@@ -630,22 +586,23 @@ public class ModelMain {
 	}
 
 	public CountryBoundaryRaster getCountryBoundaryRaster() {
-		CountryBoundaryRaster countryBoundaries = new CountryBoundaryRaster(desiredProjection, compositeCountryManager);
+		CountryBoundaryRaster countryBoundaries = new CountryBoundaryRaster(desiredProjection);
 		CountryBoundaryReader countryReader = new CountryBoundaryReader(countryBoundaries);
 		countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);
 		return countryBoundaries;
 	}
 
 	public void createCountryAgents(Collection<CompositeCountry> countryGrouping) {
-		countryAgents = new CountryAgentManager(compositeCountryManager, demandManager, countryBoundaryRaster, internationalMarket, clusterIdRaster, globalLandUseRaster);
+		countryAgents = new CountryAgentManager(demandManager, countryBoundaryRaster, internationalMarket, clusterIdRaster, globalLandUseRaster);
 		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = getInitialCropUsageData();
-		Map<CompositeCountry, Map<WoodType, WoodUsageData>> woodUsageDataMap = getInitialWoodUsageData();
+		Map<CompositeCountry, Map<WoodType, WoodUsageData>> woodUsageDataMap = new HashMap<CompositeCountry, Map<WoodType, WoodUsageData>>();
+		if (ModelConfig.IS_FORESTRY_ON) woodUsageDataMap = getInitialWoodUsageData();
 		Map<CompositeCountry, CarbonUsageData> carbonUsageDataMap = getInitialCarbonUsageData();
 		RasterSet<LandUseItem> initLU = getInitialLandUse();
+		globalLandUseRaster.putAll(initLU);
 		
 		for (CompositeCountry cc : countryGrouping) {
 			countryAgents.addForCountry(cc, cropUsageDataMap, initLU, woodUsageDataMap, carbonUsageDataMap);
-			globalLandUseRaster.putAll(initLU);
 		}
 	}
 
@@ -662,7 +619,7 @@ public class ModelMain {
 	private Map<CompositeCountry, Map<CropType, CropUsageData>> getInitialCropUsageData() {
 		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap;
 		if (ModelConfig.IS_CALIBRATION_RUN || ModelConfig.USE_INITIAL_CROP_USAGE_DATA)
-			cropUsageDataMap = new CropUsageReader(compositeCountryManager).getCommodityData();
+			cropUsageDataMap = new CropUsageReader().getCommodityData();
 		else
 			cropUsageDataMap = deserializeCropUsage();
 
@@ -672,7 +629,7 @@ public class ModelMain {
 	private Map<CompositeCountry, Map<WoodType, WoodUsageData>> getInitialWoodUsageData() {
 		Map<CompositeCountry, Map<WoodType, WoodUsageData>> woodUsageDataMap;
 		if (ModelConfig.IS_CALIBRATION_RUN) {
-			woodUsageDataMap = new WoodUsageReader(compositeCountryManager).getWoodUsageData();
+			woodUsageDataMap = new WoodUsageReader().getWoodUsageData();
 
 		} else {
 			woodUsageDataMap = deserializeWoodUsage();
@@ -684,7 +641,7 @@ public class ModelMain {
 		Map<CompositeCountry, CarbonUsageData> carbonUsageDataMap;
 		if (ModelConfig.IS_CALIBRATION_RUN) {
 			carbonUsageDataMap = new HashMap<CompositeCountry, CarbonUsageData>();
-			for (CompositeCountry cc : compositeCountryManager.getAll()) {
+			for (CompositeCountry cc : CountryManager.getInstance().getAllCompositeCountries()) {
 				CarbonUsageData cuData = new CarbonUsageData(0, 0, 0);
 				carbonUsageDataMap.put(cc, cuData);
 			}
@@ -742,7 +699,7 @@ public class ModelMain {
 		RasterSet<LandUseItem> rasterToSerialise = new RasterSet<LandUseItem>(desiredProjection);
 		for (Map.Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
 			LandUseItem newLuItem = new LandUseItem(entry.getValue());
-			newLuItem.overwriteLandCoverAreas(new HashMap<LandCoverType, LandCoverTile>()); // Clear data
+			newLuItem.overwriteLandCoverAreas(new LandCoverTile()); // Clear data
 			rasterToSerialise.put(entry.getKey(), newLuItem);
 		}
 		try {
@@ -859,13 +816,14 @@ public class ModelMain {
 			carbonFluxData = carbonFluxReader.getCarbonFluxes(globalLandUseRaster, timestep);
 		} else {
 			carbonFluxData = new CarbonFluxRasterSet(desiredProjection);
+			carbonFluxData.fillWithDefaults(globalLandUseRaster.keySet());
 		}
 	}
 	
 	/** Get wood yield data */
 	private void getWoodYieldData(Timestep timestep) {
 		if (ModelConfig.IS_FORESTRY_ON || ModelConfig.ENABLE_VEGETATION_CLEARANCE_COST) {
-			woodYieldData = woodYieldReader.getWoodYields(globalLandUseRaster, timestep, internationalMarket.getWoodPrice().getExportPrice());	
+			woodYieldData = woodYieldReader.getWoodYields(globalLandUseRaster, timestep, internationalMarket.getWoodPrices().get(WoodType.IND_ROUNDWOOD).getExportPrice());
 		} else {
 			woodYieldData = new WoodYieldRasterSet(desiredProjection);
 		}		
@@ -890,4 +848,37 @@ public class ModelMain {
 		countryAgents.serializeCarbonUsageForAll();
 		internationalMarket.serializeGlobalPrices();
 	}
+
+	// Checks for data coverage. Fills missing data with defaults. Currently this is zeroes but plan to use
+	// nearest neighbour interpolation eventually.
+	private void handleMissingData() {
+		int totalCount = 0;
+		int coverage = 0;
+		Set<RasterKey> completeKeys = new HashSet<>();
+		for (RasterKey key : globalLandUseRaster.keySet()) {
+			totalCount ++;
+			boolean hasCropYields = yieldSurfaces.containsKey(key);
+			boolean hasIrrig = currentIrrigationData.containsKey(key);
+			boolean hasWoodYields = !ModelConfig.IS_FORESTRY_ON || woodYieldData.containsKey(key);
+			boolean hasCarbonFluxes = !ModelConfig.IS_CARBON_ON || carbonFluxData.containsKey(key);
+			boolean isComplete = (hasCropYields && hasIrrig && hasWoodYields && hasCarbonFluxes);
+			if (isComplete) {
+				coverage++;
+				completeKeys.add(key);
+			}
+			if (!hasCropYields)
+				yieldSurfaces.put(key, YieldResponsesItem.getDefault());
+
+			if (!hasIrrig)
+				currentIrrigationData.put(key, IrrigationItem.getDefault());
+
+			if (!hasWoodYields)
+				woodYieldData.put(key, WoodYieldItem.getDefault());
+
+			if (!hasCarbonFluxes)
+				carbonFluxData.put(key, CarbonFluxItem.getDefault());
+		}
+		double coveragePc = (double) coverage / totalCount * 100;
+		LogWriter.println(String.format("Data coverage: %.2f%%", coveragePc), 1);
+	}
 }
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxItem.java b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
index da8979e8327780268d5a5c273ad1dee09b85acef..d020dde9df65cb9f5ee3664c0295610383f45f2f 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxItem.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxItem.java
@@ -19,19 +19,17 @@ public class CarbonFluxItem implements RasterItem, Serializable {
 	Map<LccKey, Double> carbonFluxRate = new HashMap<LccKey, Double>(); 
 	Map<LccKey, Double> carbonCreditRate = new HashMap<LccKey, Double>();
 	
-
-	public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
-		LandCoverTile tiles = landUseTiles.get(key.getFromLc());
-		double totalArea = tiles.getTotalArea(LandProtectionType.CONVERTIBLE);
+	public void calcConversionCarbonFlux(LccKey key, Double[] cFluxes, LandCoverTile landTile, int maxAge) {
+		double totalArea = landTile.getTotalArea(key.getFromLc(), LandProtectionType.CONVERTIBLE);
 		
 		if (totalArea == 0) {
 			carbonFluxRate.put(key, 0.0);
 			carbonCreditRate.put(key, 0.0);
 		} else {
 			double totalFlux = 0;
-			for (int age : tiles.getAgeKeys()) {
+			for (int age : landTile.getAgeKeys()) {
 				int ageCapped = Math.min(age, maxAge - 1);
-				totalFlux += cFluxes[ageCapped] * tiles.getArea(LandProtectionType.CONVERTIBLE, age);
+				totalFlux += cFluxes[ageCapped] * landTile.getArea(key.getFromLc(), LandProtectionType.CONVERTIBLE, age);
 			}
 			double meanFlux = totalFlux / totalArea;
 			carbonFluxRate.put(key, meanFlux);
@@ -39,50 +37,48 @@ public class CarbonFluxItem implements RasterItem, Serializable {
 		}
 	}
 	
-	public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
-		LandCoverTile tiles = landUseTiles.get(lcType);
+	public void calcEcosystemCarbonFlux(LandCoverType lcType, Double[] cFluxes, LandCoverTile landTile, int maxAge) {
 		LccKey key = new LccKey(lcType, lcType);
 		
-		double totalArea = tiles.getTotalArea(LandProtectionType.CONVERTIBLE);
+		double totalArea = landTile.getTotalArea(lcType, LandProtectionType.CONVERTIBLE);
 			
 		if (totalArea == 0) {
 			carbonFluxRate.put(key, 0.0);
 		} else {
 			// Flux rate
 			double totalFlux = 0;
-			for (int age : tiles.getAgeKeys()) {
+			for (int age : landTile.getAgeKeys()) {
 				int ageCapped = Math.min(age, maxAge - 1);
-				totalFlux += cFluxes[ageCapped] * tiles.getArea(LandProtectionType.CONVERTIBLE, age);
+				totalFlux += cFluxes[ageCapped] * landTile.getArea(lcType, LandProtectionType.CONVERTIBLE, age);
 			}
 			double meanFlux = totalFlux / totalArea;
 			carbonFluxRate.put(key, meanFlux);
 			
 			// Carbon credit rate for unchanged land cover
-			for (int age : tiles.getAgeKeys()) {
-				for (int i = age; i < ModelConfig.CARBON_HORIZON + age; i++) {
-					int ageCapped = Math.min(age, maxAge - 1);
-					totalFlux += cFluxes[ageCapped] * tiles.getArea(LandProtectionType.CONVERTIBLE, age);
-				}
-			}
-			double annualisedFlux = (totalFlux / totalArea) / ModelConfig.CARBON_HORIZON; 
-			carbonCreditRate.put(key, -annualisedFlux); // switching sign as we want +ve to mean C offset
+
+
 		}
 		
-		// Credit rate for destination land cover
-		{
-			double totalFlux = 0;
-			for (int i = 0; i < ModelConfig.CARBON_HORIZON; i++) {
-				totalFlux += cFluxes[i];
+		double fromLcCumulFlux = 0;
+		for (int age : landTile.getAgeKeys()) {
+			for (int i = age; i <= age + ModelConfig.CARBON_HORIZON; i++) {
+				fromLcCumulFlux += cFluxes[i];
 			}
-			double annualisedFlux = -totalFlux / ModelConfig.CARBON_HORIZON; // switching sign as we want +ve to mean C offset
-			for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
-				if (!fromLc.equals(lcType)) {
-					LccKey lccKey = new LccKey(fromLc, lcType);
-					double fluxSoFar = carbonCreditRate.get(lccKey);
-					double updatedFlux = fluxSoFar + annualisedFlux;
-					carbonCreditRate.put(lccKey, updatedFlux);
-				}
-			} 
+		}
+		double toLcCumulFlux = 0;
+		for (int age = 0; age <= ModelConfig.CARBON_HORIZON; age++) {
+			toLcCumulFlux += cFluxes[age];
+		}
+		
+		// flip signs as we want positive to mean net offset
+		fromLcCumulFlux *= -1;
+		toLcCumulFlux *= -1;
+		
+		for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
+			LccKey fromKey = new LccKey(lcType, lc);
+			LccKey toKey = new LccKey(lc, lcType);
+			carbonCreditRate.merge(fromKey, fromLcCumulFlux, (a, b) -> a - b);
+			carbonCreditRate.merge(toKey, toLcCumulFlux, (a, b) -> a + b);
 		}
 	}
 
@@ -123,5 +119,17 @@ public class CarbonFluxItem implements RasterItem, Serializable {
 		return carbonCreditRate;
 	}
 
+	public static CarbonFluxItem getDefault() {
+		CarbonFluxItem item = new CarbonFluxItem();
+		for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
+			for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
+				LccKey key = new LccKey(fromLc, toLc);
+				item.setCarbonFluxRate(key, 0.0);
+				item.setCarbonCreditRate(key, 0.0);
+			}
+		}
+		return item;
+	}
+
 
 }
diff --git a/src/ac/ed/lurg/carbon/CarbonFluxRasterSet.java b/src/ac/ed/lurg/carbon/CarbonFluxRasterSet.java
index 3289872a0277b89877eaf398886432c08e5ac976..e3fa2d7efcd20d8f761ca62b1a7071dc9ab1e25a 100644
--- a/src/ac/ed/lurg/carbon/CarbonFluxRasterSet.java
+++ b/src/ac/ed/lurg/carbon/CarbonFluxRasterSet.java
@@ -1,6 +1,7 @@
 package ac.ed.lurg.carbon;
 
 import java.util.Collection;
+import java.util.Set;
 
 import ac.sac.raster.RasterHeaderDetails;
 import ac.sac.raster.RasterKey;
@@ -24,4 +25,10 @@ public class CarbonFluxRasterSet extends RasterSet<CarbonFluxItem> {
 		popSubsetForKeys(subsetCarbonFluxRaster, keys);		
 		return subsetCarbonFluxRaster;
 	}
+
+	public void fillWithDefaults(Set<RasterKey> keys) {
+		for (RasterKey key : keys) {
+			put(key, CarbonFluxItem.getDefault());
+		}
+	}
 }
diff --git a/src/ac/ed/lurg/country/AbstractCountryAgent.java b/src/ac/ed/lurg/country/AbstractCountryAgent.java
index d308ca2279c5510dc9cd55f530a545560dd612ad..c76461b9b07b782f168fcdda96d6a1f915fd424f 100644
--- a/src/ac/ed/lurg/country/AbstractCountryAgent.java
+++ b/src/ac/ed/lurg/country/AbstractCountryAgent.java
@@ -24,14 +24,14 @@ public abstract class AbstractCountryAgent {
 	protected Map<CropType, CountryPrice> currentCountryPrices;
 	private Map<CropType, CountryPrice> previousCountryPrices;
 	private Map<CommodityType, Double> currentConsumerCommodityPrices;
-	protected GlobalPrice globalWoodPrice;
+	protected Map<WoodType, GlobalPrice> globalWoodPrices;
 	protected GlobalPrice globalCarbonPrice;
 	protected CountryPrice currentCarbonPrice;
-	protected CountryPrice currentTimberPrice;
+	protected Map<WoodType, CountryPrice> currentWoodPrices;
 	protected Timestep currentTimestep;
 	protected Map<CommodityType, Map<CropType, Double>> currentDemandFract;
 	protected double currentGen2EcDemand;
-	protected Map<Integer, Map<WoodType, Double>> currentWoodDemand;
+	protected Map<WoodType, Double> currentWoodDemand;
 	protected double currentCarbonDemand;
 	protected double exportTaxRate;
 	
@@ -61,7 +61,13 @@ public abstract class AbstractCountryAgent {
 				
 		currentCountryPrices = countryPrices;
 		currentCarbonPrice = new CountryPrice(globalCarbonPrice); // no trade barriers
-		currentTimberPrice = new CountryPrice(globalWoodPrice, ModelConfig.WOOD_TRADE_BARRIER);
+		Map<WoodType, CountryPrice> countryWoodPrices = new HashMap<>();
+		for (WoodType woodType : WoodType.values()) {
+			GlobalPrice globPrice = globalWoodPrices.get(woodType);
+			CountryPrice price = new CountryPrice(globPrice, ModelConfig.WOOD_TRADE_BARRIER);
+			countryWoodPrices.put(woodType, price);
+		}
+		currentWoodPrices = countryWoodPrices;
 	}
 	
 	protected void savePreviousProducerCropPrices() {
@@ -91,10 +97,10 @@ public abstract class AbstractCountryAgent {
 	
 	abstract protected CountryPrice createCountryPrices(CropType crop, GlobalPrice worldPrice);
 	
-	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice, boolean outputGamsDemand) {
+	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, Map<WoodType, GlobalPrice> woodPrices, boolean outputGamsDemand) {
 		LogWriter.println("\ncalculateCountryPricesAndDemand for " + country + " "+ currentTimestep.getTimestep(), 3);
 		currentWorldPrices = worldPrices;
-		globalWoodPrice = timberPrice;
+		globalWoodPrices = woodPrices;
 		globalCarbonPrice = carbonPrice;
 		currentDemandFract = getDemandFraction();
 		calculateCountryPrices();
@@ -244,8 +250,8 @@ public abstract class AbstractCountryAgent {
 	
 	abstract public Map<WoodType, WoodUsageData> getWoodUsageData();
 	
-	public CountryPrice getCurrentCountryWoodPrice() {
-		return currentTimberPrice;
+	public Map<WoodType, CountryPrice> getCurrentCountryWoodPrices() {
+		return currentWoodPrices;
 	}
 	
 	public CountryPrice getCurrentCountryCarbonPrice() {
diff --git a/src/ac/ed/lurg/country/AgriculturalGDPManager.java b/src/ac/ed/lurg/country/AgriculturalGDPManager.java
index d7865c5012ac15627a335fcaabac89a4da75a9ac..1ce115ca9ed0cf8e0a9772933b68d4512088ff87 100644
--- a/src/ac/ed/lurg/country/AgriculturalGDPManager.java
+++ b/src/ac/ed/lurg/country/AgriculturalGDPManager.java
@@ -32,7 +32,7 @@ public class AgriculturalGDPManager {
 		String filename = ModelConfig.GDP_FRACTIONS_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName;
+			String line, countryCode;
 			Double proportion;
 			fitReader.readLine(); // read header
 
@@ -42,10 +42,10 @@ public class AgriculturalGDPManager {
 				if (tokens.length < 2)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[COUNTRY_COL];
+				countryCode = tokens[COUNTRY_COL];
 				proportion = Double.parseDouble(tokens[FRACTION_COL]);
 				
-				SingleCountry country = CountryManager.getForName(countryName);
+				SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 				
 				fractionsMap.put(country,proportion);
 				
diff --git a/src/ac/ed/lurg/country/AnimalRateManager.java b/src/ac/ed/lurg/country/AnimalRateManager.java
index 5b6997915d4bfebb42d8c4e299492cf501d92ed0..06a2875a46d548778b31f3fa2866e0ff5871045f 100644
--- a/src/ac/ed/lurg/country/AnimalRateManager.java
+++ b/src/ac/ed/lurg/country/AnimalRateManager.java
@@ -11,12 +11,10 @@ import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.utils.StringTabularReader;
 
 public class AnimalRateManager {
-	protected CompositeCountryManager compositeCountryManager;
 	private StringTabularReader tabularReader;
 
-	public AnimalRateManager(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
-		tabularReader = new StringTabularReader(",", new String[]{"Country","plumCropItem","LivestockItem","animalRate", "prod"});
+	public AnimalRateManager() {
+		tabularReader = new StringTabularReader(",", new String[]{"Iso3","plumDemandItem","LivestockItem","animalRate", "prod"});
 		tabularReader.read(ModelConfig.ANIMAL_RATES_FILE);
 	}
 
@@ -34,7 +32,7 @@ public class AnimalRateManager {
 		Map<String, Double> animalRates = new HashMap<String, Double>();
 		Map<String, Double> prods = new HashMap<String, Double>();
 
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			Map<String, AnimalRateData> singleCMap = getAnimalRates(c, animalType);
 
 			for (Entry<String, AnimalRateData> entry : singleCMap.entrySet()) {
@@ -61,7 +59,7 @@ public class AnimalRateManager {
 		Map<String, AnimalRateData> rates = new HashMap<String, AnimalRateData>();
 
 		Map<String, String> queryMap = new HashMap<String, String>();
-		queryMap.put("Country", c.getCountryName());
+		queryMap.put("Iso3", c.getCountryCode());
 		queryMap.put("plumCropItem", animalType.getFaoName());
 
 		try {
diff --git a/src/ac/ed/lurg/country/CompositeCountryManager.java b/src/ac/ed/lurg/country/CompositeCountryManager.java
deleted file mode 100644
index 7a0930f219d59fa63e7f5111f06b88e506376bc1..0000000000000000000000000000000000000000
--- a/src/ac/ed/lurg/country/CompositeCountryManager.java
+++ /dev/null
@@ -1,104 +0,0 @@
-package ac.ed.lurg.country;
-
-import java.util.Collection;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.List;
-import java.util.Map;
-import java.util.Set;
-import java.util.stream.Collectors;
-
-import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.demand.BaseConsumpManager;
-import ac.ed.lurg.utils.CollectionHelper;
-import ac.ed.lurg.utils.StringTabularReader;
-
-public class CompositeCountryManager {
-
-	private Map<String, CompositeCountry> mapFromName = new HashMap<String, CompositeCountry>();
-	private Map<SingleCountry, CompositeCountry> mapFromSingleCountry;
-	private Map<CompositeCountry, List<SingleCountry>> mapFromCompositeCountry;
-	private StringTabularReader countryGroups;
-	private Set<String> countryExclusionList = new HashSet<String>();
-	
-	public CompositeCountryManager(BaseConsumpManager baseConsumpManager) {
-		super();
-		if (ModelConfig.PREDEFINED_COUNTRY_GROUPING) {
-			countryGroups = new StringTabularReader(",", new String[]{"Country", "PlumGroup"});
-			countryGroups.read(ModelConfig.COUNTRY_GROUPING_FILE);
-		}
-		if (ModelConfig.EXCLUDE_COUNTRIES_IN_LIST) {
-			countryExclusionList = new StringTabularReader(",", new String[] {"Country"}).read(ModelConfig.EXCLUDED_COUNTRIES_FILE)
-					.stream().map(m -> m.get("Country")).collect(Collectors.toSet());
-		}
-		populate(baseConsumpManager);
-	}
-	
-	private String getCountryGroup(SingleCountry aCountry) {		
-		Map<String, String> queryMap = new HashMap<String, String>();
-		queryMap.put("Country", aCountry.getCountryName());
-		Map<String, String> row = countryGroups.querySingle(queryMap);
-		return row.get("PlumGroup");
-	}
-	
-	
-	private void populate(BaseConsumpManager baseConsumpManager) {
-		mapFromSingleCountry = new HashMap<SingleCountry, CompositeCountry>();
-		
-		for (SingleCountry c : baseConsumpManager.getAllCountries()) {
-			CompositeCountry cc;
-			
-			if (countryExclusionList.contains(c.getCountryName()))
-				continue;			
-			else if (ModelConfig.PREDEFINED_COUNTRY_GROUPING)
-				cc = getForName(getCountryGroup(c),c.getRegion());
-			else if (baseConsumpManager.getPopForGrouping(c) < ModelConfig.POPULATION_AGGREG_LIMIT)
-				cc = getForName(c.getRegion() + "_other", c.getRegion());
-			else 
-				cc = getForName(c.getCountryName(),c.getRegion());
-			
-			mapFromSingleCountry.put(c, cc);
-		}
-		
-		mapFromCompositeCountry = CollectionHelper.invertMap(mapFromSingleCountry);
-	}
-
-	private CompositeCountry getForName(String name, String region) {
-		CompositeCountry cc = mapFromName.get(name);
-		
-		if (cc == null) {
-			cc = new CompositeCountry(name, region);
-			mapFromName.put(name, cc);
-		}
-		return cc;
-	}
-	
-	public CompositeCountry getForSingleCountry(SingleCountry c) {
-		return mapFromSingleCountry.get(c);
-	}
-	
-	public Collection<CompositeCountry> getAll() {
-		return mapFromName.values();
-	}
-	
-	public Collection<SingleCountry> getAllSingleCountries() {
-		Collection<SingleCountry> worldCountries = getAll()
-				.stream()
-				.flatMap(cc -> getAllForCompositeCountry(cc)
-				.stream())
-				.collect(Collectors.toList());
-		
-		return worldCountries;
-	}
-	
-	public Collection<SingleCountry> getAllForCompositeCountry(CompositeCountry cc) {
-		return mapFromCompositeCountry.get(cc);
-	}
-	
-	public String getIncomeGroup(SingleCountry aCountry) {		
-		Map<String, String> queryMap = new HashMap<String, String>();
-		queryMap.put("Country", aCountry.getCountryName());
-		Map<String, String> row = countryGroups.querySingle(queryMap);
-		return row.get("IncomeGroup");
-	}
-}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 59f94d17498dc38188eadf4cb8fca5a6f4cb1d9f..81d7aca6c20e032a11815174419d3cbd39c39364 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -13,7 +13,6 @@ import ac.ed.lurg.country.gams.GamsRasterOptimiser;
 import ac.ed.lurg.country.gams.GamsRasterOutput;
 import ac.ed.lurg.demand.AbstractDemandManager;
 import ac.ed.lurg.forestry.WoodYieldItem;
-import ac.ed.lurg.forestry.WoodYieldRasterSet;
 import ac.ed.lurg.landuse.CarbonUsageData;
 import ac.ed.lurg.landuse.CropUsageData;
 import ac.ed.lurg.landuse.IrrigationItem;
@@ -77,8 +76,8 @@ public class CountryAgent extends AbstractCountryAgent {
 				IrrigationItem irrigItem = irrigData.get(key);
 				
 				LandUseItem luItem = previousGamsRasterOutput.getLandUses().get(key);
-				double totalLand = luItem.getTotalLandCoverArea()-luItem.getLandCoverArea(LandCoverType.URBAN);
-				double protectedAreaFrac = (totalLand <= 0) ? 0.0 : luItem.getTotalLandCoverArea(LandProtectionType.PROTECTED) / totalLand;
+				double totalLand = luItem.getTotalLandArea()-luItem.getLandCoverArea(LandCoverType.URBAN);
+				double protectedAreaFrac = (totalLand <= 0) ? 0.0 : luItem.getTotalLandArea(LandProtectionType.PROTECTED) / totalLand;
 				
 				WoodYieldItem wyItem = woodYieldData.get(key);
 				
@@ -112,10 +111,10 @@ public class CountryAgent extends AbstractCountryAgent {
 	public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData, 
 			Map<CropType, GlobalPrice> worldPrices, RasterSet<CarbonFluxItem> carbonFluxData, 
 			RasterSet<WoodYieldItem> woodYieldData, Map<LccKey, Double> conversionCosts,
-			GlobalPrice carbonPrice, GlobalPrice woodPrice) {
+			GlobalPrice carbonPrice, Map<WoodType, GlobalPrice> globalWoodPrices) {
 
 		// project demand
-		calculateCountryPricesAndDemand(worldPrices, carbonPrice, woodPrice, false);
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, globalWoodPrices, false);
 		resetLccAreas();
 		
 		if (ModelConfig.APPLY_EXPORT_TAXES && !ModelConfig.IS_CALIBRATION_RUN) {
@@ -154,9 +153,10 @@ public class CountryAgent extends AbstractCountryAgent {
 		throw new RuntimeException("Skipping optimisation of country " + country);
 	}
 	
-	public void recalculateDemandAndNetImports(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice) {
+	public void recalculateDemandAndNetImports(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice,
+											   Map<WoodType, GlobalPrice> woodPrices) {
 		shockGDP();
-		calculateCountryPricesAndDemand(worldPrices, carbonPrice, timberPrice, true);
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, globalWoodPrices, true);
 		updateNetImportsFromProdAndDemand(currentProjectedDemand, currentDemandFract, previousGamsRasterOutput.getCropUsageData());
 	}
 	
@@ -227,11 +227,11 @@ public class CountryAgent extends AbstractCountryAgent {
 			carbonTradeConstraint = new TradeConstraint(0, 0);
 		}
 
-		// Timber import/export constraints
 		Map<WoodType, TradeConstraint> woodTradeConstraints = new HashMap<WoodType, TradeConstraint>();
+
+		// Timber import/export constraints
 		Map<WoodType, WoodUsageData> woodUsageData = previousGamsRasterOutput.getWoodUsageData();
-		double totalWoodProd = woodUsageData.values().stream().mapToDouble(o -> o.getHarvest()).reduce(0.0, Double::sum);
-		
+
 		for (Map.Entry<WoodType, WoodUsageData> entry : woodUsageData.entrySet()) {
 			WoodType woodType = entry.getKey();
 			if (!ModelConfig.IS_FORESTRY_ON) {
@@ -239,33 +239,17 @@ public class CountryAgent extends AbstractCountryAgent {
 				continue;
 			}
 			WoodUsageData woodUsage = entry.getValue();
-			
 			double baseTrade = woodUsage.getNetImport();
-			
-			// Make sure country can import sufficient wood
-			// Some countries produce wood from irrigated forests e.g. Egypt which are not modelled here.
-			// Therefore we need to adjust the initial FAO imports/exports.
-			
-			if ((ModelConfig.IS_CALIBRATION_RUN || ModelConfig.USE_INITIAL_CROP_USAGE_DATA) && currentTimestep.isInitialTimestep()) {
-				double prodFract = woodUsage.getHarvest() / totalWoodProd; 
-				double potentialMaxYield = WoodYieldRasterSet.getPotentialTimberProduction(woodYieldData, previousGamsRasterOutput.getLandUses()) * 0.9;
-				double potentialAdjYield = potentialMaxYield * prodFract;
-				double currentDemand = currentWoodDemand.get(0).get(woodType);
-				double importsNeeded = currentDemand - potentialAdjYield;
-				if (baseTrade > 0 & importsNeeded > baseTrade)
-					baseTrade = importsNeeded;
-			}
-			
 			double changeUp;
 			double changeDown;
-			
+
 			double maxOfProdOrSupply = woodUsage.getHarvest() + Math.max(baseTrade, 0);
-			
+
 			if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
 				changeUp = changeDown = 0.0;
 			} else {
-				changeUp = maxOfProdOrSupply * globalWoodPrice.getMaxImportChange();
-				changeDown = maxOfProdOrSupply * globalWoodPrice.getMaxExportChange();
+				changeUp = maxOfProdOrSupply * globalWoodPrices.get(woodType).getMaxImportChange();
+				changeDown = maxOfProdOrSupply * globalWoodPrices.get(woodType).getMaxExportChange();
 			}
 
 			woodTradeConstraints.put(woodType, new TradeConstraint(baseTrade - changeDown, baseTrade + changeUp));
@@ -273,7 +257,7 @@ public class CountryAgent extends AbstractCountryAgent {
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
 				previousGamsRasterOutput.getCropUsageData(), currentDemandFract, subsidyRates, currentGen2EcDemand, currentCarbonDemand, currentCarbonPrice, carbonTradeConstraint, 
-				currentWoodDemand, currentTimberPrice, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData());	
+				currentWoodDemand, currentWoodPrices, woodTradeConstraints, previousGamsRasterOutput.getWoodUsageData());
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, 
 				woodYieldData, carbonFluxData, conversionCosts, countryLevelInputs);
 
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index c989277d128c5ea47a141380b20a5153ca2635da..771ec23089087bea10d7b4a4063abbb34780f6e8 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -56,7 +56,7 @@ public class CountryAgentManager {
 	private Collection<CraftyCountryAgent> craftyCountryAgents = new ArrayList<CraftyCountryAgent>();
 	private CraftyProdManager craftyManager;
 	
-	public CountryAgentManager (CompositeCountryManager compositeCountryManager, AbstractDemandManager demandManager, 
+	public CountryAgentManager (AbstractDemandManager demandManager, 
 			CountryBoundaryRaster countryBoundaryRaster, InternationalMarket internationalMarket,
 			RasterSet<IntegerRasterItem> clusterIdRaster, RasterSet<LandUseItem> globalLandUseRaster) {
 		this.demandManager = demandManager;
@@ -64,8 +64,8 @@ public class CountryAgentManager {
 		this.clusterIdRaster = clusterIdRaster;
 		this.internationalMarket = internationalMarket;
 		this.globalLandUseRaster = globalLandUseRaster;
-		tradeBarrierManager = new TradeManager(compositeCountryManager);
-		subsidyRateManager = new SubsidyRateManager(compositeCountryManager);
+		tradeBarrierManager = new TradeManager();
+		subsidyRateManager = new SubsidyRateManager();
 		exportRestrictionManager = new ExportRestrictionManager();
 		minMaxNetImportManager = new MinMaxNetImportManager();
 		craftyManager = new CraftyProdManager();
@@ -147,7 +147,7 @@ public class CountryAgentManager {
 
 					try {
 						ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), carbonFluxData, woodYieldData, 
-								conversionCosts, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice());
+								conversionCosts, internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrices());
 
 						// update global rasters
 						globalLandUseRaster.putAll(ca.getLandUses());
@@ -173,13 +173,14 @@ public class CountryAgentManager {
 
 		if (craftyCountryAgents.size() > 0) {
 			craftyManager.updateWithCraftyData(craftyCountryAgents, timestep, internationalMarket.getWorldPrices(), 
-					internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice());  // this will wait for the marker file from CRAFTY
+					internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrices());  // this will wait for the marker file from CRAFTY
 		}
 	}
 	
 	public void recalculateDemandAndNetImportsForAll() {
 		for (CountryAgent ca : gamsCountryAgents)
-			ca.recalculateDemandAndNetImports(internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice(), internationalMarket.getWoodPrice());
+			ca.recalculateDemandAndNetImports(internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice(),
+					internationalMarket.getWoodPrices());
 	}	
 	
 	public void updateProtectedAreasForAll(RasterSet<LandUseItem> newProtectedAreas) {
diff --git a/src/ac/ed/lurg/country/CountryBoundaryRaster.java b/src/ac/ed/lurg/country/CountryBoundaryRaster.java
index 223a8d91e80482c5773e3b0073d4aee48970cbd3..952ecdc11cd70015cae820654c52ec5aee301210 100644
--- a/src/ac/ed/lurg/country/CountryBoundaryRaster.java
+++ b/src/ac/ed/lurg/country/CountryBoundaryRaster.java
@@ -12,11 +12,9 @@ import ac.sac.raster.RasterSet;
 
 public class CountryBoundaryRaster extends RasterSet<CountryBoundaryItem> {
 	Collection<CompositeCountry> countryGrouping;
-	CompositeCountryManager compositeCountryManager;
 	
-	public CountryBoundaryRaster(RasterHeaderDetails header, CompositeCountryManager compositeCountryManager) {
+	public CountryBoundaryRaster(RasterHeaderDetails header) {
 		super(header);
-		this.compositeCountryManager = compositeCountryManager;
 	}
 	
 	private static final long serialVersionUID = -8449000692429399253L;
@@ -36,7 +34,7 @@ public class CountryBoundaryRaster extends RasterSet<CountryBoundaryItem> {
 				if (c == null)
 					continue;
 	
-				CompositeCountry cc = compositeCountryManager.getForSingleCountry(c);
+				CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(c);
 				List<RasterKey> keys = countryMap.get(cc);
 	
 				if (keys == null) {
diff --git a/src/ac/ed/lurg/country/CountryBoundaryReader.java b/src/ac/ed/lurg/country/CountryBoundaryReader.java
index 22c8d9f647d8b105ee5279bf51905041e23d36aa..955b47bccd9dcd0fe5c4d615402a4b99ed325867 100644
--- a/src/ac/ed/lurg/country/CountryBoundaryReader.java
+++ b/src/ac/ed/lurg/country/CountryBoundaryReader.java
@@ -16,8 +16,8 @@ public class CountryBoundaryReader extends AbstractTabularRasterReader<CountryBo
 
 	@Override
 	public void setData(RasterKey key, CountryBoundaryItem item, Map<String, Double> rowValues) {
-		int isoCode = (int) getValueForCol(rowValues, "FaoCode");
-		SingleCountry c = CountryManager.getForCode(isoCode);
+		int isoCode = (int) getValueForCol(rowValues, "M49");
+		SingleCountry c = CountryManager.getInstance().getForM49Code(isoCode);
 		item.setCountry(c);
 	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryManager.java b/src/ac/ed/lurg/country/CountryManager.java
index 0707d113cabf6b5e0da5f0aed5f59c3f68d87ddb..6de952025b6593bdb2d14c96ebddc73d68af1d9b 100644
--- a/src/ac/ed/lurg/country/CountryManager.java
+++ b/src/ac/ed/lurg/country/CountryManager.java
@@ -1,74 +1,85 @@
 package ac.ed.lurg.country;
 
-import java.io.BufferedReader;
-import java.io.FileReader;
+import java.util.ArrayList;
+import java.util.Collection;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
-
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.utils.StringTabularReader;
 
 
 public class CountryManager {
-	private static final int NAME_COL = 0; 
-	private static final int CHAR_CODE_COL = 1; 
-	private static final int REGION_COL = 2; 
-	private static final int INCOME_COL = 3;
-	private static final int NUM_CODE_COL = 4; 
 
-	private final Map<String, SingleCountry> nameMap = new HashMap<String, SingleCountry>();
-	private final Map<Integer, SingleCountry> numCodeMap = new HashMap<Integer, SingleCountry>();
+	private Map<String, SingleCountry> codeMap = new HashMap<String, SingleCountry>();
+	private Map<Integer, SingleCountry> m49CodeMap = new HashMap<Integer, SingleCountry>();
+	private Map<SingleCountry, CompositeCountry> mapFromSingleCountry = new HashMap<SingleCountry, CompositeCountry>();
+	private Map<CompositeCountry, List<SingleCountry>> mapFromCompositeCountry = new HashMap<CompositeCountry, List<SingleCountry>>();
+	private List<SingleCountry> countryExclusionList = new ArrayList<SingleCountry>();
 	
-	private static CountryManager instance = null;
-
-	private static CountryManager getInstance() {
-		if(instance == null) {
-			instance = new CountryManager();
-		}
+	private static final CountryManager instance = new CountryManager();
+	
+	public static CountryManager getInstance() {
 		return instance;
 	}
 	
-	protected CountryManager() {
-		String filename = ModelConfig.COUNTRY_CODES_FILE;
-		try {
-			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, name, charCode, region, incomeGroup;
-			int numCode;
-			fitReader.readLine(); // read header
-
-			while ((line=fitReader.readLine()) != null) {
-				String[] tokens = line.split(",");
-				
-				if (tokens.length < 4)
-					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
-				
-				name = tokens[NAME_COL];
-				charCode = tokens[CHAR_CODE_COL];
-				region = tokens[REGION_COL];
-				numCode = Integer.parseInt(tokens[NUM_CODE_COL]);
-				incomeGroup = tokens[INCOME_COL];
-
-				SingleCountry c = new SingleCountry(name, charCode, numCode, region, incomeGroup);
-				nameMap.put(name, c);
-				numCodeMap.put(numCode, c);
-			} 
-			fitReader.close(); 
-		
-		} catch (Exception e) {
-			LogWriter.printlnError("Failed in reading commodity demand fits");
-			LogWriter.print(e);
+	private CountryManager() {
+		readCountries();
+	}
+	
+	private void readCountries() {
+		String filename = ModelConfig.COUNTRIES_FILE;
+		StringTabularReader reader = new StringTabularReader(",", new String[] {"Iso3", "Area", "M49", "Population", "FaoArea", "PlumGroup"});
+		reader.read(filename);
+		List<Map<String, String>> rowList = reader.getRowList();
+		for (Map<String, String> row : rowList) {
+			String countryName = row.get("Area");
+			String countryCode = row.get("Iso3");
+			Integer m45Code = Integer.parseInt(row.get("M49"));
+			String region = "";
+			String incomeGroup = "";
+			String plumGroup = row.get("PlumGroup");
+			
+			SingleCountry sc = new SingleCountry(countryName, countryCode, region, incomeGroup);
+			CompositeCountry cc;
+			if (ModelConfig.PREDEFINED_COUNTRY_GROUPING) {
+				cc = new CompositeCountry(plumGroup, region);
+			} else {
+				cc = new CompositeCountry(region, region);
+			}
+			codeMap.put(countryCode, sc);
+			m49CodeMap.put(m45Code, sc);
+			mapFromSingleCountry.put(sc, cc);
+			List<SingleCountry> scList = mapFromCompositeCountry.computeIfAbsent(cc, k -> new ArrayList<SingleCountry>());
+			scList.add(sc);
 		}
-		LogWriter.println("Processed " + filename + ", create " + nameMap.size() + " country codes", 2);
-
 	}
 	
-	public static SingleCountry getForName(String name) {
-		SingleCountry c = getInstance().nameMap.get(name);
+
+	public SingleCountry getForCode(String isoCode) {
+		SingleCountry c = codeMap.get(isoCode);
 		return c;
 	}
-
-	public static SingleCountry getForCode(int isoCode) {
-		SingleCountry c = getInstance().numCodeMap.get(isoCode);
+	
+	public SingleCountry getForM49Code(int m49) {
+		SingleCountry c = m49CodeMap.get(m49);
 		return c;
 	}
+	
+	public CompositeCountry getForSingleCountry(SingleCountry c) {
+		return mapFromSingleCountry.get(c);
+	}
+	
+	public Collection<SingleCountry> getAllForCompositeCountry(CompositeCountry cc) {
+		return mapFromCompositeCountry.get(cc);
+	}
+	
+	public Collection<SingleCountry> getAllSingleCountries() {
+		return codeMap.values();
+	}
+	
+	public Collection<CompositeCountry> getAllCompositeCountries() {
+		return mapFromCompositeCountry.keySet();
+	}
+
 }
diff --git a/src/ac/ed/lurg/country/LandCoverChangeItem.java b/src/ac/ed/lurg/country/LandCoverChangeItem.java
index 348b93e29e3fb72f8274bbd3c4196d185ea8b478..e464e43946f215be0f2a56b76923bac7e11af1a3 100644
--- a/src/ac/ed/lurg/country/LandCoverChangeItem.java
+++ b/src/ac/ed/lurg/country/LandCoverChangeItem.java
@@ -5,6 +5,7 @@ import ac.ed.lurg.types.LandCoverType;
 public class LandCoverChangeItem {
 	private LandCoverType fromLandCover;
 	private LandCoverType toLandCover;
+	private int age;
 	private double area;
 	
 	public LandCoverChangeItem() {}
@@ -15,6 +16,13 @@ public class LandCoverChangeItem {
 		this.area = area;
 	}
 
+	public LandCoverChangeItem(LandCoverType fromLandCover, LandCoverType toLandCover, int age, double area) {
+		this.fromLandCover = fromLandCover;
+		this.toLandCover = toLandCover;
+		this.age = age;
+		this.area = area;
+	}
+
 	public LandCoverType getFromLandCover() {
 		return fromLandCover;
 	}
@@ -23,6 +31,10 @@ public class LandCoverChangeItem {
 		return toLandCover;
 	}
 
+	public int getAge() {
+		return age;
+	}
+
 	public double getArea() {
 		return area;
 	}
diff --git a/src/ac/ed/lurg/country/PriceCapManager.java b/src/ac/ed/lurg/country/PriceCapManager.java
index a204a599348a5729a9271ed8a84b412fc67ec1dd..fef59c47d63243df5c2812ae00288d4875c29439 100644
--- a/src/ac/ed/lurg/country/PriceCapManager.java
+++ b/src/ac/ed/lurg/country/PriceCapManager.java
@@ -7,6 +7,7 @@ import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.types.WoodType;
 import ac.ed.lurg.utils.LogWriter;
 
 public class PriceCapManager {
@@ -95,10 +96,16 @@ public class PriceCapManager {
 		return updatedPrices;
 	}
 	
-	public GlobalPrice capWoodPrices(GlobalPrice price) {
-		double factor = calcCappedPriceFactor(price.getExportPrice(), woodPriceCaps.get(PriceCapType.MIN_PRICE), woodPriceCaps.get(PriceCapType.MAX_PRICE));
-		GlobalPrice updatedPrice = price.createPriceAdjustedByFactor(factor);
-		return updatedPrice;
+	public Map<WoodType, GlobalPrice> capWoodPrices(Map<WoodType, GlobalPrice> woodPrices) {
+		Map<WoodType, GlobalPrice> updatedPrices = new HashMap<>();
+		for (WoodType woodType : WoodType.values()) {
+			GlobalPrice price = woodPrices.get(woodType);
+			double factor = calcCappedPriceFactor(price.getExportPrice(), woodPriceCaps.get(PriceCapType.MIN_PRICE),
+					woodPriceCaps.get(PriceCapType.MAX_PRICE));
+			updatedPrices.put(woodType, price.createPriceAdjustedByFactor(factor));
+		}
+
+		return updatedPrices;
 	}
 	
 	public GlobalPrice capCarbonPrices(GlobalPrice price) {
diff --git a/src/ac/ed/lurg/country/SingleCountry.java b/src/ac/ed/lurg/country/SingleCountry.java
index 759dd7c8b135aff8aed14b98ec31a57777f408de..5ecdd97ef73ec4a4b7cbb02bb2083313df6637ce 100644
--- a/src/ac/ed/lurg/country/SingleCountry.java
+++ b/src/ac/ed/lurg/country/SingleCountry.java
@@ -2,27 +2,21 @@ package ac.ed.lurg.country;
 
 public class SingleCountry {
 
+	private String iso3;
 	private String countryName;
-	private String iso3CharCode;
-	private int iso3NumCode;
 	private String region;
 	private String incomeGroup;
 	
-	public SingleCountry(String countryName, String iso3CharCode, int iso3NumCode, String region, String incomeGroup) {
+	public SingleCountry(String countryName, String iso3CharCode, String region, String incomeGroup) {
 		super();
 		this.countryName = countryName;
-		this.iso3CharCode = iso3CharCode;
-		this.iso3NumCode = iso3NumCode;
+		this.iso3 = iso3CharCode;
 		this.region = region;
 		this.incomeGroup = incomeGroup;
 	}
 
-	public String getIso3CharCode() {
-		return iso3CharCode;
-	}
-
-	public int getIso3NumCode() {
-		return iso3NumCode;
+	public String getCountryCode() {
+		return iso3;
 	}
 
 	public String getCountryName() {
diff --git a/src/ac/ed/lurg/country/SubsidyRateManager.java b/src/ac/ed/lurg/country/SubsidyRateManager.java
index f896bf165bf618e00f791d3091a93c205cff5a28..5c4757693c556a9ac8f282794041111e999f1f6e 100644
--- a/src/ac/ed/lurg/country/SubsidyRateManager.java
+++ b/src/ac/ed/lurg/country/SubsidyRateManager.java
@@ -10,26 +10,24 @@ import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.utils.StringTabularReader;
 
 public class SubsidyRateManager {
-	protected CompositeCountryManager compositeCountryManager;
 	private StringTabularReader tabularReader;
 
-	public SubsidyRateManager(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
-		tabularReader = new StringTabularReader(",", new String[]{"country", "crop", "rate"});
+	public SubsidyRateManager() {
+		tabularReader = new StringTabularReader(",", new String[]{"Iso3", "PlumCropItem", "SubsPerT"});
 		tabularReader.read(ModelConfig.SUBSIDY_RATE_FILE);
 	}
 
 	public Map<CropType, Double> getSubsidyRates(CompositeCountry cc) {
 
 		Map<CropType, Double> subsidyRates = new HashMap<CropType, Double>();
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			Map<CropType, Double> singleCMap = getSubsidyRates(c);
 
 			for (CropType crop : CropType.getAllItems()) {
 				Double newRateD = singleCMap.get(crop);
 				Double prevRateD = subsidyRates.get(crop);
 
-				// just takes max which isn't right. Should instead weight by production or something else
+				// TODO just takes max which isn't right. Should instead weight by production or something else
 				double newRate = Math.max(getDoubleOrZero(newRateD), getDoubleOrZero(prevRateD));
 				subsidyRates.put(crop, newRate);
 			}
@@ -46,19 +44,19 @@ public class SubsidyRateManager {
 
 		for (CropType crop : CropType.getAllItems()) {
 			Map<String, String> queryMap = new HashMap<String, String>();
-			queryMap.put("country", c.getCountryName());
-			queryMap.put("crop", crop.getFaoName());
+			queryMap.put("Iso3", c.getCountryCode());
+			queryMap.put("PlumCropItem", crop.getFaoName());
 			try {
 				double subsidy = 0 ;
 				List<Map<String, String>> rows = tabularReader.query(queryMap);
 				if (rows.size() == 1) {
-					String subsidyS = rows.get(0).get("rate");
+					String subsidyS = rows.get(0).get("SubsPerT");
 					subsidy = Double.valueOf(subsidyS);
 				}
 				rates.put(crop, subsidy);
 			}
 			catch (Exception e) {
-					LogWriter.println("Problem finding cereal demand: " + crop.getFaoName() + ", " + c.getCountryName());
+					LogWriter.println("Problem finding subsidy rate: " + crop.getFaoName() + ", " + c.getCountryName());
 			}
 		}
 		return rates;
diff --git a/src/ac/ed/lurg/country/TradeManager.java b/src/ac/ed/lurg/country/TradeManager.java
index b32dbfcf89eb0f6519100f592756f7375aefdd7a..0cc3242fc8a0cffedd6d63dfb54ec9720f9c597e 100644
--- a/src/ac/ed/lurg/country/TradeManager.java
+++ b/src/ac/ed/lurg/country/TradeManager.java
@@ -19,15 +19,15 @@ public class TradeManager {
 	private static final int ITEM_COL = 1; 
 	private static final int TRADE_BARRIER_COL = 2; 
 	
-	private CompositeCountryManager compositeCountryManager;
 	private Map<CompositeCountry, CropToDoubleMap> tradeMap = new HashMap<CompositeCountry, CropToDoubleMap>();
 
-	public TradeManager(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
+	public TradeManager() {
 		if(ModelConfig.ACTIVE_TRADE_BARRIERS)
 			read();
 	}
 	
+	// TODO weird summing
+	
 	public CropToDoubleMap getTradeBarriers(CompositeCountry cc) {
 
 		if (tradeMap.get(cc) == null) { 
@@ -53,7 +53,7 @@ public class TradeManager {
 	
 		try {
 			BufferedReader reader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, itemName;
+			String line, countryCode, itemName;
 			double barrierValue;
 			reader.readLine(); // read header
 
@@ -65,15 +65,15 @@ public class TradeManager {
 				
 				try {
 
-					countryName = tokens[COUNTRY_COL];
+					countryCode = tokens[COUNTRY_COL];
 					itemName = tokens[ITEM_COL];
 					barrierValue = Double.parseDouble(tokens[TRADE_BARRIER_COL]);
 					
-					SingleCountry country = CountryManager.getForName(countryName);
+					SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 					
 					CropType crop = CropType.getForFaoName(itemName);
 					
-					CompositeCountry cc = compositeCountryManager.getForSingleCountry(country);
+					CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(country);
 					if (cc == null) {
 						continue;
 					}
@@ -100,15 +100,6 @@ public class TradeManager {
 			LogWriter.print(e);
 		} 
 		LogWriter.println("Processed " + filename + ", create " + tradeMap.size() + " barriers entries");
-		
-		for (Entry<CompositeCountry, CropToDoubleMap> entry : tradeMap.entrySet()) {
-			
-			entry.getValue().divideBy(compositeCountryManager.getAllForCompositeCountry(entry.getKey()).size());
-		
-//			LogWriter.println("number of countries " + compositeCountryManager.getAllForCompositeCountry(entry.getKey()).size());
-//			LogWriter.println("crops " + entry.getValue());
-		
-		}
 
 	}
 	
diff --git a/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java b/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
index b4fcc09e0d80498d4aa7953a20b6a07155e54224..1b9ac19e644d51baf05066efb21dc8682b408a37 100644
--- a/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
+++ b/src/ac/ed/lurg/country/crafty/CraftyCountryAgent.java
@@ -29,9 +29,10 @@ public class CraftyCountryAgent extends AbstractCountryAgent {
 		return cropUsageData;
 	}
 
-	public void updateProduction(Map<CropType, CropUsageData> cropUsageMap, Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, GlobalPrice timberPrice) {
+	public void updateProduction(Map<CropType, CropUsageData> cropUsageMap, Map<CropType, GlobalPrice> worldPrices,
+								 GlobalPrice carbonPrice, Map<WoodType, GlobalPrice> woodPrices) {
 		this.cropUsageData = cropUsageMap;
-		calculateCountryPricesAndDemand(worldPrices, carbonPrice, timberPrice, false);
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, woodPrices, false);
 		updateNetImportsFromProdAndDemand(currentProjectedDemand, currentDemandFract, cropUsageMap);
 	}
 	
diff --git a/src/ac/ed/lurg/country/crafty/CraftyProdManager.java b/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
index 12265272db3d0c290dadab4c802b211ae8b6184b..92797a8ee16c43153b3d1ae48ef7ac1d72fe1309 100644
--- a/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
+++ b/src/ac/ed/lurg/country/crafty/CraftyProdManager.java
@@ -13,6 +13,7 @@ import ac.ed.lurg.country.CompositeCountry;
 import ac.ed.lurg.country.GlobalPrice;
 import ac.ed.lurg.landuse.CropUsageData;
 import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.types.WoodType;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.utils.StringTabularReader;
 import ac.ed.lurg.utils.WatchForFile;
@@ -33,8 +34,8 @@ public class CraftyProdManager {
 		return craftyCountries;
 	}
 
-	public void updateWithCraftyData(Collection<CraftyCountryAgent> craftyCountryAgents, Timestep timestep, Map<CropType, GlobalPrice> worldPrices, 
-			GlobalPrice carbonPrice, GlobalPrice timberPrice) {
+	public void updateWithCraftyData(Collection<CraftyCountryAgent> craftyCountryAgents, Timestep timestep, Map<CropType, GlobalPrice> worldPrices,
+									 GlobalPrice carbonPrice, Map<WoodType, GlobalPrice> woodPrices) {
 		String rootDir = ModelConfig.CRAFTY_PRODUCTION_DIR + File.separator + timestep.getYear();
 		long startTime = System.currentTimeMillis();
 		
@@ -76,7 +77,7 @@ public class CraftyProdManager {
 				if (cropUsageMap.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 " + cropUsageMap.size());
 				}
-				cca.updateProduction(cropUsageMap, worldPrices, carbonPrice, timberPrice);
+				cca.updateProduction(cropUsageMap, worldPrices, carbonPrice, woodPrices);
 			}
 			catch (Exception e) {
 				LogWriter.printlnWarning("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 f621e4fabfd08139e109c5dff045c8b5f3a73f02..bfd2bd02db59f2c04f5084807d3817e40699fb0d 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -26,16 +26,16 @@ public class GamsCountryInput {
 	private double carbonDemand;
 	private CountryPrice carbonPrice;
 	private TradeConstraint carbonTradeConstraint;
-	private Map<Integer, Map<WoodType, Double>> woodDemand;
-	private CountryPrice woodPrice;
+	private Map<WoodType, Double> woodDemand;
+	private Map<WoodType, CountryPrice> woodPrices;
 	private Map<WoodType, TradeConstraint> woodTradeConstraints;
 	private Map<WoodType, WoodUsageData> previousWoodUsageData;
 
 	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<Integer, Map<WoodType, Double>> woodDemand, 
-			CountryPrice woodPrice, Map<WoodType, TradeConstraint> woodTradeConstraints, Map<WoodType, WoodUsageData> previousWoodUsageData) {
+			double carbonDemand, CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint, Map<WoodType, Double> woodDemand,
+			Map<WoodType, CountryPrice> woodPrices, Map<WoodType, TradeConstraint> woodTradeConstraints, Map<WoodType, WoodUsageData> previousWoodUsageData) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
@@ -49,7 +49,7 @@ public class GamsCountryInput {
 		this.carbonPrice = carbonPrice;
 		this.carbonTradeConstraint = carbonTradeConstraint;
 		this.woodDemand = woodDemand;
-		this.woodPrice = woodPrice;
+		this.woodPrices = woodPrices;
 		this.woodTradeConstraints = woodTradeConstraints;
 		this.previousWoodUsageData = previousWoodUsageData;
 	}
@@ -120,12 +120,12 @@ public class GamsCountryInput {
 		return carbonTradeConstraint;
 	}
 
-	public Map<Integer, Map<WoodType, Double>> getWoodDemand() {
+	public Map<WoodType, Double> getWoodDemand() {
 		return woodDemand;
 	}
 
-	public CountryPrice getWoodPrice() {
-		return woodPrice;
+	public Map<WoodType, CountryPrice> getWoodPrices() {
+		return woodPrices;
 	}
 
 	public Map<WoodType, TradeConstraint> getWoodTradeConstraints() {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index c86ce1119d60db62d4fb46f43e56913f2955538d..12846d67f6d944dce83c232a9f029a222ed66a03 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -7,12 +7,8 @@ import java.nio.file.Files;
 import java.nio.file.Path;
 import java.nio.file.Paths;
 import java.nio.file.StandardCopyOption;
-import java.util.ArrayList;
-import java.util.HashMap;
-import java.util.Map;
+import java.util.*;
 import java.util.Map.Entry;
-import java.util.Set;
-import java.util.Vector;
 
 import com.gams.api.GAMSDatabase;
 import com.gams.api.GAMSException;
@@ -362,28 +358,32 @@ public class GamsLocationOptimiser {
 		addScalar(inDB, "maxIrrigChange", ModelConfig.MAX_IRRIG_CHANGE, 5);
 		
 		// Forestry
-		
-		// Not simulating production of different wood products so need to aggregate
-		double woodMinNetImport = 0;
-		double woodMaxNetImport = 0;
+		GAMSParameter maxRotaIntensityP = inDB.addParameter("maxRotationIntensity", 1);
+
 		Map<WoodType, TradeConstraint> woodTradeConstraints = countryInput.getWoodTradeConstraints();
-		for (TradeConstraint tc: woodTradeConstraints.values()) {
-			woodMinNetImport += tc.getMinConstraint();
-			woodMaxNetImport += tc.getMaxConstraint();
+		Map<WoodType, CountryPrice> woodPrices = countryInput.getWoodPrices();
+		Map<WoodType, Double> woodDemandMap = countryInput.getWoodDemand();
+
+		for (WoodType woodType : WoodType.values()) {
+			// Demand
+			setGamsParamValue(demandP.addRecord(woodType.getName()), woodDemandMap.get(woodType), -1);
+			// Trade limits
+			TradeConstraint tc = woodTradeConstraints.get(woodType);
+			setGamsParamValue(minTradeP.addRecord(woodType.getName()), tc.getMinConstraint(), 5);
+			setGamsParamValue(maxTradeP.addRecord(woodType.getName()), tc.getMaxConstraint(), 5);
+			// Prices
+			setGamsParamValue(importPriceP.addRecord(woodType.getName()), woodPrices.get(woodType).getImportPrice(), 5);
+			setGamsParamValue(exportPriceP.addRecord(woodType.getName()), woodPrices.get(woodType).getExportPrice(), 5);
+			// Rotation limits
+			setGamsParamValue(maxRotaIntensityP.addRecord(woodType.getName()), 1.0 / woodType.getMinRotation(), 6);
 		}
-		setGamsParamValue(minTradeP.addRecord("wood"), woodMinNetImport, 5);
-		setGamsParamValue(maxTradeP.addRecord("wood"), woodMaxNetImport, 5);
-		
-		setGamsParamValue(importPriceP.addRecord("wood"), countryInput.getWoodPrice().getImportPrice(), 5);
-		setGamsParamValue(exportPriceP.addRecord("wood"), countryInput.getWoodPrice().getExportPrice(), 5);
 		
 		addScalar(inDB, "forestManagementCost", ModelConfig.FOREST_MANAGEMENT_COST, 5);
 		addScalar(inDB, "vegClearingCostRate", ModelConfig.VEGETATION_CLEARING_COST, 5);
 		
 		// Yield from timber forest rotation
 		Map<Integer, ? extends WoodYieldData> woodYieldData = inputData.getWoodYields();
-		
-		GAMSParameter woodYieldLucP = inDB.addParameter("woodYieldLUC", 3);
+
 		GAMSParameter woodRespMaxYieldP = inDB.addParameter("woodYieldMax", 1);
 		GAMSParameter woodRespSlopeP = inDB.addParameter("woodYieldParam", 1);
 		
@@ -394,19 +394,6 @@ public class GamsLocationOptimiser {
 			String locString = Integer.toString(locationId);
 			WoodYieldData wYield = entry.getValue();
 
-			for (LandCoverType fromLc : LandCoverType.getForestedTypes()) {
-				for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {	
-				    Vector<String> w = new Vector<String>();
-					w.add(fromLc.getName());
-					w.add(toLc.getName());
-					w.add(locString);
-					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);
@@ -415,20 +402,6 @@ public class GamsLocationOptimiser {
 			}			
 			
 		}
-		
-		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);
@@ -475,7 +448,8 @@ public class GamsLocationOptimiser {
 				String fromName = key.getFromLc().getName();
 				String toName = key.getToLc().getName();
 				double baseCost = entry.getValue();
-				Double vegClearCost = woodYieldData.get(location).getYieldLuc(key.getFromLc(), key.getToLc()) * ModelConfig.VEGETATION_CLEARING_COST;
+				double totalVegYield = woodYieldData.get(location).getYieldLuc(key);
+				Double vegClearCost = totalVegYield * ModelConfig.VEGETATION_CLEARING_COST;
 				double totalCost = (vegClearCost.isNaN()) ? baseCost : baseCost + vegClearCost;
 				Vector<String> v = new Vector<String>();
 				v.add(fromName);
@@ -606,7 +580,7 @@ public class GamsLocationOptimiser {
 			
 		// Land cover change
 		GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
-		Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges = new HashMap<Integer, ArrayList<LandCoverChangeItem>>();
+		Map<Integer, List<LandCoverChangeItem>> landCoverChanges = new HashMap<>();
 		
 		for (GAMSVariableRecord rec : varLandCoverChange) {
 			String fromLcStr = rec.getKey(0);
@@ -622,37 +596,37 @@ public class GamsLocationOptimiser {
 			LandCoverType fromLc = LandCoverType.getForName(fromLcStr);
 			LandCoverType toLc = LandCoverType.getForName(toLcStr);	
 			
-			ArrayList<LandCoverChangeItem> changesList = landCoverChanges.computeIfAbsent(locId, k -> new ArrayList<LandCoverChangeItem>());
+			List<LandCoverChangeItem> changesList = landCoverChanges.computeIfAbsent(locId, k -> new ArrayList<>());
 			changesList.add(new LandCoverChangeItem(fromLc, toLc, change));
 		}
 		
-		// Timber harvest
-		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.
-		Map<WoodType, WoodUsageData> newWoodUsageMap = new HashMap<WoodType, WoodUsageData>();
-		Map<WoodType, WoodUsageData> previousWoodUsageMap = inputData.getCountryInput().getPreviousWoodUsageData();
-		double previousNetImportSum = previousWoodUsageMap.values().stream().mapToDouble(o -> o.getNetImport()).sum();
-		double netWoodImportChange = netWoodImport - previousNetImportSum;
-		double previousWoodHarvestSum = previousWoodUsageMap.values().stream().mapToDouble(o -> o.getHarvest()).sum();
-		int numWoodTypes = WoodType.values().length;
-		for (WoodType wt : WoodType.values()) {
-			double previousNetImport = previousWoodUsageMap.get(wt).getNetImport();
-			double newNetImport = previousNetImport + netWoodImportChange / numWoodTypes; // assume equal split for trade change
-			double previousHarvest = previousWoodUsageMap.get(wt).getHarvest();
-			
-			// assuming equal split if no previous harvest
-			double newHarvest = (previousWoodHarvestSum > 0) ? woodSupplyRota * (previousHarvest / previousWoodHarvestSum) : woodSupplyRota / numWoodTypes;
-			newWoodUsageMap.put(wt, new WoodUsageData(newHarvest, newNetImport, 0));
+		// Wood production
+		Map<WoodType, WoodUsageData> newWoodUsageMap = new HashMap<>();
+		GAMSParameter parWoodSupplyRota = outDB.getParameter("woodSupplyRota");
+		for (WoodType w : WoodType.values()) {
+			double supply = getParmValue(parWoodSupplyRota, w.getName());
+			double netImports = getParmValue(parmNetImports, w.getName());
+			newWoodUsageMap.put(w, new WoodUsageData(supply, netImports));
 		}
-		
-		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);
+
+		Map<Integer, Map<WoodType, Double>> rotationPeriods = new HashMap<>();
+		GAMSVariable varRota = outDB.getVariable("rotationIntensity");
+		for (GAMSVariableRecord rec : varRota) {
+			WoodType woodType = WoodType.getForName(rec.getKey(0));
+			int loc = Integer.parseInt(rec.getKey(1));
+			double intensity = rec.getLevel();
+			Map<WoodType, Double> map = rotationPeriods.computeIfAbsent(loc, k -> new HashMap<>());
+			map.put(woodType, 1 / intensity);
+		}
+
+		Map<Integer, Map<WoodType, Double>> forestAreaMap = new HashMap<>();
+		GAMSVariable varForestArea = outDB.getVariable("timberForestArea");
+		for (GAMSVariableRecord rec : varForestArea) {
+			WoodType woodType = WoodType.getForName(rec.getKey(0));
+			int loc = Integer.parseInt(rec.getKey(1));
+			double forestArea = rec.getLevel();
+			Map<WoodType, Double> map = forestAreaMap.computeIfAbsent(loc, k -> new HashMap<>());
+			map.put(woodType, forestArea);
 		}
 		
 		// Carbon
@@ -660,9 +634,8 @@ public class GamsLocationOptimiser {
 		double netCarbonImport = getParmValue(parmNetImports, "carbonCredits");
 		CarbonUsageData carbonUsageData = new CarbonUsageData(carbonSequestered, netCarbonImport, 0.0);
 
-		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, carbonUsageData, 
-				newWoodUsageMap, rotationPeriods);
-		return results;
+		return new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, carbonUsageData,
+				newWoodUsageMap, rotationPeriods, forestAreaMap);
 	}
 
 	private double getParmValue(GAMSParameter aParm, String itemName) {
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 309b8eeee60156f97fe6017aea8a318e1063b3b4..4b3bce17f452d0a455ecfc710d597a9234d15bad 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -1,6 +1,7 @@
 package ac.ed.lurg.country.gams;
 
 import java.util.ArrayList;
+import java.util.List;
 import java.util.Map;
 
 import com.gams.api.GAMSGlobals.ModelStat;
@@ -18,17 +19,19 @@ public class GamsLocationOutput {
 	
 	Map<Integer, LandUseItem> landUses;  // data mapped from id (not raster)
 	private Map<CropType, CropUsageData> cropUsageData;
-	Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges;
+	Map<Integer, List<LandCoverChangeItem>> landCoverChanges;
 	private CarbonUsageData carbonUsageData;
 	private Map<WoodType, WoodUsageData> woodUsageData;
-	private Map<Integer, Double> rotationPeriods;
+	private Map<Integer, Map<WoodType, Double>> rotationPeriods;
+	private Map<Integer, Map<WoodType, Double>> forestAreas;
 	
-	public GamsLocationOutput(ModelStat status, 
-			Map<Integer, LandUseItem> landUses, 
-			Map<CropType, CropUsageData> cropUsageData,
-			Map<Integer, ArrayList<LandCoverChangeItem>> landCoverChanges,
-			CarbonUsageData carbonUsageData, Map<WoodType, WoodUsageData> woodUsageData,
-			Map<Integer, Double> rotationPeriods) {
+	public GamsLocationOutput(ModelStat status,
+                              Map<Integer, LandUseItem> landUses,
+                              Map<CropType, CropUsageData> cropUsageData,
+                              Map<Integer, List<LandCoverChangeItem>> landCoverChanges,
+                              CarbonUsageData carbonUsageData, Map<WoodType, WoodUsageData> woodUsageData,
+                              Map<Integer, Map<WoodType, Double>> rotationPeriods,
+							  Map<Integer, Map<WoodType, Double>> forestAreas) {
 		super();
 		this.status = status;
 		this.landUses = landUses;
@@ -37,6 +40,7 @@ public class GamsLocationOutput {
 		this.carbonUsageData = carbonUsageData;
 		this.woodUsageData = woodUsageData;
 		this.rotationPeriods = rotationPeriods;
+		this.forestAreas = forestAreas;
 	}
 	
 	public ModelStat getStatus() {
@@ -50,7 +54,7 @@ public class GamsLocationOutput {
 		return cropUsageData;
 	}
 	
-	public Map<Integer, ArrayList<LandCoverChangeItem>> getLandCoverChanges() {
+	public Map<Integer, List<LandCoverChangeItem>> getLandCoverChanges() {
 		return landCoverChanges;
 	}
 
@@ -62,8 +66,11 @@ public class GamsLocationOutput {
 		return woodUsageData;
 	}
 
-	public Map<Integer, Double> getRotationPeriods() {
+	public Map<Integer, Map<WoodType, Double>> getRotationPeriods() {
 		return rotationPeriods;
 	}
-	
+
+	public Map<Integer, Map<WoodType, Double>> getForestAreas() {
+		return forestAreas;
+	}
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 492cd248af63543fdcc581a6e814e173726c3ac9..5e355f306d442fd18f0e131cd02dd7206d9d210a 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -1,29 +1,15 @@
 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.*;
 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.ForestManager;
 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;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.types.LandProtectionType;
-import ac.ed.lurg.types.YieldType;
+import ac.ed.lurg.landuse.*;
+import ac.ed.lurg.types.*;
 import ac.ed.lurg.utils.LazyTreeMap;
 import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.yield.YieldRaster;
@@ -63,26 +49,6 @@ 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());
@@ -119,7 +85,7 @@ public class GamsRasterOptimiser {
 		double suitableArea=0, protectedArea=0;
 		for (LandUseItem a : areaRaster.values()) {
 			if (a != null) {
-				protectedArea += a.getTotalLandCoverArea(LandProtectionType.PROTECTED);
+				protectedArea += a.getTotalLandArea(LandProtectionType.PROTECTED);
 				suitableArea += a.getSuitableArea();
 			}
 		}
@@ -130,6 +96,32 @@ public class GamsRasterOptimiser {
 	private RasterSet<LandUseItem> allocAreas(Map<Integer, ? extends LandUseItem> prevAreasAgg, GamsLocationOutput gamsOutput, int year) {
 		RasterSet<LandUseItem> newLandUseRaster = createWithSameLandCovers(rasterInputData.getPreviousLandUses());
 
+		// Order is important here
+
+		// Disaggregate cluster land cover changes to raster grid
+		Map<RasterKey, List<LandCoverChangeItem>> lcChanges = disaggregateLandCoverChanges(gamsOutput.getLandCoverChanges(),
+				newLandUseRaster);
+
+		// Calculate amount of wood produced from land cover change
+		Map<WoodType, Double> lccWoodAmount = ForestManager.handleLccWood(rasterInputData.getWoodYields(),
+				lcChanges, rasterInputData.getPreviousLandUses());
+
+		// Allocate land cover changes
+		allocLandCoverChanges(newLandUseRaster, lcChanges);
+
+		// Calculate land cover changes for timber forest rotation and clear-cut
+		Map<RasterKey, List<LandCoverChangeItem>> lcChangesRota = ForestManager.handleRotaWood(
+				newLandUseRaster, rasterInputData.getCountryInput().getWoodDemand(),
+				rasterInputData.getWoodYields(), gamsOutput, mapping);
+
+		// Add wood produced from LCC to output
+		Map<WoodType, WoodUsageData> woodUsageMap = gamsOutput.getWoodUsageData();
+		for (WoodType wType : lccWoodAmount.keySet()) {
+			double amount = lccWoodAmount.get(wType);
+			WoodUsageData woodUsage = woodUsageMap.get(wType);
+			woodUsage.setLucHarvest(amount);
+		}
+
 		for (Map.Entry<Integer, LandUseItem> entry : gamsOutput.getLandUses().entrySet()) {
 			Integer locId = entry.getKey();
 			LandUseItem newLandUseAggItem = entry.getValue();
@@ -143,19 +135,6 @@ public class GamsRasterOptimiser {
 			}
 			
 			RasterSet<LandUseItem> landUseItemsForLocation = newLandUseRaster.createSubsetForKeys(keys);
-			
-			ArrayList<LandCoverChangeItem> landCoverChanges = gamsOutput.getLandCoverChanges().get(locId);
-			
-			if (landCoverChanges != null) {	
-				// Allocate land cover change to grid cells
-				for (LandCoverChangeItem lccItem : landCoverChanges) {
-					LandCoverType fromLC = lccItem.getFromLandCover();
-					LandCoverType toLC = lccItem.getToLandCover();
-					double change = lccItem.getArea();
-					allocAllLandCropsTop(newLandUseRaster, keys, toLC, fromLC, change, locId, false);
-					
-				}
-			}
 
 			for (RasterKey key : keys) {
 				LandUseItem newLandUseItem = newLandUseRaster.get(key);
@@ -172,128 +151,64 @@ public class GamsRasterOptimiser {
 
 		return newLandUseRaster;
 	}
-	
-	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;
-
-		// reverse direction if negative
-		if (change < 0) {
-			change = -change;
-			toType = fromLC;
-			fromType = toLC;
-		}
 
-		double shortfall = allocAllLandCrops(newLandUseRaster, keys, toType, fromType, change, natToTimb);
+	private Map<RasterKey, List<LandCoverChangeItem>> disaggregateLandCoverChanges(
+			Map<Integer, List<LandCoverChangeItem>> aggLandCoverChanges, RasterSet<LandUseItem> luRaster) {
 
-		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, boolean natToTimb) {
-		if (DEBUG) LogWriter.println("allocAllLandCrops: from " + fromType + " to " + toType + " change " + change, 3);
-		if (change < 0.00001)
-			return 0;
-
-		double totalShortfall = 0;
-		Set<RasterKey> keysWithSpace = new HashSet<RasterKey>();
-		double totalUnprotectedLC = 0;
-		
-		for (RasterKey key : keys) {
-			LandUseItem newLandUseItem = newLandUseRaster.get(key);
-			totalUnprotectedLC += newLandUseItem.getLandCoverArea(fromType, LandProtectionType.CONVERTIBLE);
-		}
+		Map<RasterKey, List<LandCoverChangeItem>> rasterChanges = new HashMap<>();
 
-		for (RasterKey key : keys) {
-			//if (DEBUG) LogWriter.println("Processing raster key " + key);
-			LandUseItem newLandUseItem = newLandUseRaster.get(key);
-
-			if (newLandUseItem!=null) {
-				double cellChange = (totalUnprotectedLC > 0) ? change * newLandUseItem.getLandCoverArea(fromType, LandProtectionType.CONVERTIBLE)/totalUnprotectedLC : change / keys.size();
-				double shortfall = 0;
-				if (natToTimb) {
-					shortfall = newLandUseItem.convertNaturalToTimber(cellChange);
-				} else {
-					shortfall = newLandUseItem.moveAreas(toType, fromType, cellChange);
+		for (Integer locId : aggLandCoverChanges.keySet()) {
+			List<LandCoverChangeItem> changes = aggLandCoverChanges.get(locId);
+			Set<RasterKey> keys = new HashSet<>();
+			for (RasterKey key : mapping.keySet()) {
+				if (mapping.get(key).getInt() == locId) {
+					keys.add(key);
 				}
-				if (shortfall == 0)
-					keysWithSpace.add(key);
-				else
-					totalShortfall += shortfall;
-				
-				//LogWriter.println(key.toString() + "|" + cellChange);
 			}
-		}
-
-		if (DEBUG) LogWriter.println("allocAllLandCrops: from " + fromType + " to " + toType + " change " + change + " totalShortfall " + totalShortfall, 3);
-		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);
+			for (LandCoverChangeItem item : changes) {
+				LandCoverType fromLc = item.getFromLandCover();
+				LandCoverType toLc = item.getToLandCover();
+				double totalChangeArea = item.getArea();
+				double totalArea = 0;
+				for (RasterKey key : keys) {
+					totalArea += luRaster.get(key).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
+				}
+				for (RasterKey key : keys) {
+					double area = luRaster.get(key).getLandCoverArea(fromLc, LandProtectionType.CONVERTIBLE);
+					double fraction = area / totalArea;
+					// Splitting up the area changes in clusters proportionally by areas in raster cells
+					double changedArea = fraction * totalChangeArea;
+					if (changedArea < 2 * Double.MIN_VALUE) {
+						continue;
 					}
+					LandCoverChangeItem lccItem = new LandCoverChangeItem(fromLc, toLc, changedArea);
+					List<LandCoverChangeItem> lccList = rasterChanges.computeIfAbsent(key, k -> new ArrayList<>());
+					lccList.add(lccItem);
 				}
 			}
 		}
-
-		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 rasterChanges;
+	}
+	
+	private void allocLandCoverChanges(RasterSet<LandUseItem> newLandUseRaster,
+									  Map<RasterKey, List<LandCoverChangeItem>> landCoverChanges) {
+		for (RasterKey key : landCoverChanges.keySet()) {
+			List<LandCoverChangeItem> changes = landCoverChanges.get(key);
+			for (LandCoverChangeItem lccItem : changes) {
+				LandCoverType fromLc = lccItem.getFromLandCover();
+				LandCoverType toLc = lccItem.getToLandCover();
+				double areaChanged = lccItem.getArea();
+				double shortfall = newLandUseRaster.get(key).moveAreas(toLc, fromLc, areaChanged);
+				if (shortfall > 0.00001) {
+					LogWriter.printlnError("This should never happen, due to GAMS constraint. Not able to incorporate all changes: from " +
+							fromLc + " to " + toLc + " " + key + ": " + shortfall);
 				}
 			}
 		}
-		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.
 		YieldRaster yieldRaster = rasterInputData.getYields();
 		RasterSet<LandUseItem> landUseRaster = rasterInputData.getPreviousLandUses();
 
@@ -305,18 +220,24 @@ public class GamsRasterOptimiser {
 
 		// look for inconsistencies
 		for (Entry<RasterKey, YieldResponsesItem> entry : yieldRaster.entrySet()) {
-
 			YieldResponsesItem yresp = entry.getValue();
+			if (yresp == null) {
+				int foo = 1;
+			}
 			RasterKey key = entry.getKey();
 			CropType crop = CropType.WHEAT;
+			double fertMaxIrrigMax = yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop);
+			double fertMaxNoIrrig = yresp.getYield(YieldType.FERT_MAX_NO_IRRIG, crop);
+			double fertNoneIrrigMax = yresp.getYield(YieldType.NO_FERT_IRRIG_MAX, crop);
 
-			if (yresp != null && (yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) < yresp.getYield(YieldType.FERT_MAX_NO_IRRIG, crop) || yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) < yresp.getYield(YieldType.NO_FERT_IRRIG_MAX, crop))) {
-				logWarningWithCoord("Inconsistency F only:" + yresp.getYield(YieldType.FERT_MAX_NO_IRRIG, crop) + ", I only" + yresp.getYield(YieldType.NO_FERT_IRRIG_MAX, crop) + ", max " + yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) + " at ", key, yieldRaster);
+			if (fertMaxIrrigMax < fertMaxNoIrrig || fertNoneIrrigMax < fertMaxIrrigMax) {
+				logWarningWithCoord("Inconsistency F only:" + fertMaxNoIrrig + ", I only" + fertNoneIrrigMax
+						+ ", max " + fertMaxIrrigMax + " at ", key, yieldRaster);
 			}
 		}
 
 		// Aggregation
-		LazyTreeMap<Integer, YieldResponsesItem> aggregatedYields = new LazyTreeMap<Integer, YieldResponsesItem>() { 
+		LazyTreeMap<Integer, YieldResponsesItem> aggregatedYields = new LazyTreeMap<Integer, YieldResponsesItem>() {
 			protected YieldResponsesItem createValue() { return new YieldResponsesItem(); }
 		};
 		LazyTreeMap<Integer, LandUseItem> aggregatedAreas = new LazyTreeMap<Integer, LandUseItem>() { 
@@ -332,42 +253,18 @@ public class GamsRasterOptimiser {
 			protected WoodYieldData createValue() { return new WoodYieldData(); }
 		};
 
-		int yRespFound = 0, yRespMissing = 0;
-		int mappingsFound = 0, mappingsMissing = 0;
-
 		for (RasterKey key : landUseRaster.keySet()) {
-			if (mapping.get(key) == null) {
-				mappingsMissing++;
-				continue;
-			} else {
-				mappingsFound++;
-			}
-			
 			int clusterId = mapping.get(key).getInt();
 
-			YieldResponsesItem yresp = yieldRaster.get(key);
-			if (yresp == null) {
-				yRespMissing++;
-				//logErrorWithCoord("No YieldResponsesItem for ", key, yieldRaster);
-			}
-			else {
-				yRespFound++;
-				//logErrorWithCoord("YieldResponsesItem found for ", key, yieldRaster);
-			}
-			
 			LandUseItem landUseItem  = landUseRaster.get(key);
-			if (landUseItem == null) {
-				LogWriter.printlnError("GamsRasterOptimiser: Can't find landUseItem for " + key);
-				continue;
-			}
-
+			YieldResponsesItem yresp = yieldRaster.get(key);
 			IrrigationItem irrigItem = irrigRaster.get(key);
 			WoodYieldItem woodYieldItem = woodYieldRaster.get(key);
 			CarbonFluxItem carbonFluxItem = carbonFluxRaster.get(key);
 
 			YieldResponsesItem aggYResp = aggregatedYields.lazyGet(clusterId);
 			LandUseItem aggLandUse = aggregatedAreas.lazyGet(clusterId);
-			IrrigationItem aggIrig = aggregatedIrrigCosts.lazyGet(clusterId);
+			IrrigationItem aggIrrig = aggregatedIrrigCosts.lazyGet(clusterId);
 			WoodYieldData aggWYield = aggregatedWoodYields.lazyGet(clusterId);
 			CarbonFluxItem aggCFlux = aggregatedCarbonFluxes.lazyGet(clusterId);
 
@@ -375,103 +272,86 @@ public class GamsRasterOptimiser {
 			double suitableAreaSoFar = aggLandUse.getSuitableArea();
 
 			// Irrigation cost
-			if (irrigItem!= null) {
-				aggIrig.setIrrigCost( aggregateMean(aggIrig.getIrrigCost(), suitableAreaSoFar, irrigItem.getIrrigCost(), suitableAreaThisTime));
-				aggIrig.setIrrigConstraint(aggregateMean(aggIrig.getIrrigConstraint(), suitableAreaSoFar, irrigItem.getIrrigConstraint(), suitableAreaThisTime));
-				//LogWriter.println(id + ", " + key + ", " + avgIrigCost.getIrrigCost() + ", from " + areaSoFar + ", " + suitableAreaThisTime + ", " + irrigCost.getIrrigCost());
-			}
+			aggIrrig.setIrrigCost( aggregateMean(aggIrrig.getIrrigCost(), suitableAreaSoFar, irrigItem.getIrrigCost(), suitableAreaThisTime));
+			aggIrrig.setIrrigConstraint(aggregateMean(aggIrrig.getIrrigConstraint(), suitableAreaSoFar, irrigItem.getIrrigConstraint(), suitableAreaThisTime));
 
 			// Aggregate carbon fluxes and wood yields
-			if (carbonFluxItem != null) {
-				for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
-					for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
-						double areaThisTime = landUseItem.getLandCoverArea(fromLc);
-						double areaSoFar = aggLandUse.getLandCoverArea(toLc);
-						
-						double cFluxThisTime = carbonFluxItem.getCarbonFluxRate(fromLc, toLc);
-						double cFluxSoFar = aggCFlux.getCarbonFluxRate(fromLc, toLc);
-						double cFluxAgg = aggregateMean(cFluxSoFar, areaSoFar, cFluxThisTime, areaThisTime);
-						aggCFlux.setCarbonFluxRate(fromLc, toLc, cFluxAgg);
-						
-						double cCreditThisTime = carbonFluxItem.getCarbonCreditRate(fromLc, toLc);
-						double cCreditSoFar = aggCFlux.getCarbonCreditRate(fromLc, toLc);
-						double cCreditAgg = aggregateMean(cCreditSoFar, areaSoFar, cCreditThisTime, areaThisTime);
-						aggCFlux.setCarbonCreditRate(fromLc, toLc, cCreditAgg);
-					}
-				}	
-			}			
-			
-			if (woodYieldItem != null) {
-				// wood from LUC
-				for (Map.Entry<LccKey, Map<Integer, Double>> wyEntry : woodYieldItem.getYieldLucMap().entrySet()) {
-					LccKey lccKey = wyEntry.getKey();
-					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 wYieldAgg = wYieldSoFar + wYieldThisTime;
-					aggWYield.setYieldLuc(lccKey, 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);
+			for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
+				for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
+					double areaThisTime = landUseItem.getLandCoverArea(fromLc);
+					double areaSoFar = aggLandUse.getLandCoverArea(toLc);
+
+					double cFluxThisTime = carbonFluxItem.getCarbonFluxRate(fromLc, toLc);
+					double cFluxSoFar = aggCFlux.getCarbonFluxRate(fromLc, toLc);
+					double cFluxAgg = aggregateMean(cFluxSoFar, areaSoFar, cFluxThisTime, areaThisTime);
+					aggCFlux.setCarbonFluxRate(fromLc, toLc, cFluxAgg);
+
+					double cCreditThisTime = carbonFluxItem.getCarbonCreditRate(fromLc, toLc);
+					double cCreditSoFar = aggCFlux.getCarbonCreditRate(fromLc, toLc);
+					double cCreditAgg = aggregateMean(cCreditSoFar, areaSoFar, cCreditThisTime, areaThisTime);
+					aggCFlux.setCarbonCreditRate(fromLc, toLc, cCreditAgg);
 				}
+			}
 
+			// wood from LUC
+			for (LccKey lccKey : woodYieldItem.getYieldLucMap().keySet()) {
+				double areaThisTime = landUseItem.getLandCoverArea(lccKey.getFromLc(), LandProtectionType.CONVERTIBLE);
+				double areaSoFar = aggLandUse.getLandCoverArea(lccKey.getFromLc(), LandProtectionType.CONVERTIBLE);
+				double wYieldThisTime = woodYieldItem.getAggregatedYieldLuc(lccKey, landUseItem);
+				double wYieldSoFar = aggWYield.getYieldLuc(lccKey);
+				double wYieldAgg = aggregateMean(wYieldSoFar, areaSoFar, wYieldThisTime, areaThisTime);
+				aggWYield.setYieldLuc(lccKey, wYieldAgg);
 			}
 
-			// Crops yields and area fractions
-			if (yresp != null) {
-				for (CropType crop : CropType.getNonMeatTypes()) {
-					if (!crop.equals(CropType.SETASIDE)) {
-						if (irrigItem!= null) {							
-							double irrigMax = irrigItem.getMaxIrrigAmount(crop);
-
-							if (Double.isNaN(irrigMax))
-								logWarningWithCoord("Can't find irrig max amount for ", key, yieldRaster, crop);
-							else 
-								aggIrig.setMaxIrrigAmount(crop, aggregateMean(aggIrig.getMaxIrrigAmount(crop), suitableAreaSoFar, irrigMax, suitableAreaThisTime));
-						}
+			for (Map.Entry<Integer, Double> entry : woodYieldItem.getYieldResponseMap().entrySet()) {
+				double areaThisTime = landUseItem.getLandCoverArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE);
+				double areaSoFar = aggLandUse.getLandCoverArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE);
+				int age = entry.getKey();
+				double yieldThisTime = entry.getValue();
+				double yieldSoFar = aggWYield.getYieldRota(age);
+				double yieldAgg = aggregateMean(yieldSoFar, areaSoFar, yieldThisTime, areaThisTime);
+				aggWYield.setYieldRota(age, yieldAgg);
+			}
 
-						for (YieldType yieldType : YieldType.values()) {
-							double y = yresp.getYield(yieldType, crop);
-							if (Double.isNaN(y))
-								logWarningWithCoord("Problem getting yield for type:" + yieldType + ", ", key, yieldRaster, crop);
-							else {
-								double yieldSoFar = aggYResp.getYield(yieldType, crop);
-								double updatedYield = aggregateMean(yieldSoFar, suitableAreaSoFar, y, suitableAreaThisTime);
-								aggYResp.setYield(yieldType, crop, updatedYield);
-							}
+			// Crops yields and area fractions
+			for (CropType crop : CropType.getNonMeatTypes()) {
+				if (!crop.equals(CropType.SETASIDE)) {
+
+					double irrigMax = irrigItem.getMaxIrrigAmount(crop);
+
+					if (Double.isNaN(irrigMax))
+						logWarningWithCoord("Can't find irrig max amount for ", key, yieldRaster, crop);
+					else
+						aggIrrig.setMaxIrrigAmount(crop, aggregateMean(aggIrrig.getMaxIrrigAmount(crop), suitableAreaSoFar, irrigMax, suitableAreaThisTime));
+
+					for (YieldType yieldType : YieldType.values()) {
+						double y = yresp.getYield(yieldType, crop);
+						if (Double.isNaN(y))
+							logWarningWithCoord("Problem getting yield for type:" + yieldType + ", ", key, yieldRaster, crop);
+						else {
+							double yieldSoFar = aggYResp.getYield(yieldType, crop);
+							double updatedYield = aggregateMean(yieldSoFar, suitableAreaSoFar, y, suitableAreaThisTime);
+							aggYResp.setYield(yieldType, crop, updatedYield);
 						}
+					}
 
-						double s = yresp.getShockRate(crop);
-						if (!Double.isNaN(s)) {
-							double shockSoFar = aggYResp.getShockRate(crop);
-							double updatedShock = aggregateMean(shockSoFar, suitableAreaSoFar, s, suitableAreaThisTime);
-							aggYResp.setShockRate(crop, updatedShock);
-						}
+					double s = yresp.getShockRate(crop);
+					if (!Double.isNaN(s)) {
+						double shockSoFar = aggYResp.getShockRate(crop);
+						double updatedShock = aggregateMean(shockSoFar, suitableAreaSoFar, s, suitableAreaThisTime);
+						aggYResp.setShockRate(crop, updatedShock);
 					}
+				}
 
-					double areaSoFar = aggLandUse.getCropArea(crop);
-					double areaThisTime = landUseItem.getCropArea(crop);
-					double updateCropland = (landUseItem.getLandCoverArea(LandCoverType.CROPLAND) + aggLandUse.getLandCoverArea(LandCoverType.CROPLAND));
+				double areaSoFar = aggLandUse.getCropArea(crop);
+				double areaThisTime = landUseItem.getCropArea(crop);
+				double updateCropland = (landUseItem.getLandCoverArea(LandCoverType.CROPLAND) + aggLandUse.getLandCoverArea(LandCoverType.CROPLAND));
 
-					if (updateCropland > 0)
-						aggLandUse.setCropFraction(crop, (areaSoFar + areaThisTime) / updateCropland);
+				if (updateCropland > 0)
+					aggLandUse.setCropFraction(crop, (areaSoFar + areaThisTime) / updateCropland);
 
-					Intensity intensityThisTime = landUseItem.getIntensity(crop);
-					aggLandUse.setIntensity(crop, intensityThisTime);  // intensity currently is always the same within a location, so can just that any.  If this changes need to average here.
-				}
+				Intensity intensityThisTime = landUseItem.getIntensity(crop);
+				aggLandUse.setIntensity(crop, intensityThisTime);  // intensity currently is always the same within a location, so can just that any.  If this changes need to average here.
 			}
 
 			// Land covers areas
@@ -486,9 +366,6 @@ public class GamsRasterOptimiser {
 			}			
 		}	
 
-		LogWriter.println("Mapping: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + mappingsFound + ", countMissing=" + mappingsMissing, 3);
-		LogWriter.println("YieldResponsesItem: " + rasterInputData.getCountryInput().getCountry() + ", countFound=" + yRespFound + ", countMissing=" + yRespMissing, 3);
-
 		checkedTotalAreas(landUseRaster, "before");
 		checkedTotalAreas(aggregatedAreas, "after");
 
diff --git a/src/ac/ed/lurg/demand/AbstractDemandManager.java b/src/ac/ed/lurg/demand/AbstractDemandManager.java
index d72c4b857be33e6079f008de9adbb38d6c332869..8abebcc9ace3341789abcbc1d11931c5018fd4b9 100644
--- a/src/ac/ed/lurg/demand/AbstractDemandManager.java
+++ b/src/ac/ed/lurg/demand/AbstractDemandManager.java
@@ -6,7 +6,7 @@ import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
+import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
@@ -16,24 +16,22 @@ import ac.ed.lurg.utils.LogWriter;
 public abstract class AbstractDemandManager implements Serializable {
 
 	private static final long serialVersionUID = 8622023691732185744L;
-	protected transient CompositeCountryManager compositeCountryManager;
 	protected transient CalorieManager calorieManager;
 	protected transient BioenergyDemandManager bioenergyDemandManager;
 	protected transient CerealFractionsManager cerealFractionsManager;
 	protected CarbonDemandManager carbonDemandManager;
 	protected WoodDemandManager woodDemandManager;
 
-	public AbstractDemandManager(CompositeCountryManager compositeCountryManager, CalorieManager calorieManager) {
-		setup(compositeCountryManager, calorieManager);
+	public AbstractDemandManager(CalorieManager calorieManager) {
+		setup(calorieManager);
 	}
 	
-	protected void setup(CompositeCountryManager compositeCountryManager, CalorieManager calorieManager) {
-		this.compositeCountryManager = compositeCountryManager;
+	protected void setup(CalorieManager calorieManager) {
 		this.calorieManager = calorieManager; 
 		bioenergyDemandManager = new BioenergyDemandManager();
-		cerealFractionsManager = new CerealFractionsManager(compositeCountryManager);
-		carbonDemandManager = new CarbonDemandManager();
-		woodDemandManager = new WoodDemandManager();
+		cerealFractionsManager = new CerealFractionsManager();
+		if (ModelConfig.IS_CARBON_ON) carbonDemandManager = new CarbonDemandManager();
+		if (ModelConfig.IS_FORESTRY_ON) woodDemandManager = new WoodDemandManager();
 	}
 
 	public Map<CommodityType, Double> getDemand(CompositeCountry cc, int year, Map<CommodityType, Double> prices, boolean outputGamsDemand) {
@@ -44,7 +42,7 @@ public abstract class AbstractDemandManager implements Serializable {
 		Map<CommodityType, Double> demandMap = new HashMap<CommodityType, Double>();
 		Map<CommodityType, Double> foodDemandMap;
 
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 
 			try {
 			foodDemandMap = getFoodDemand(c, year, prices, outputGamsDemand);
@@ -85,7 +83,7 @@ public abstract class AbstractDemandManager implements Serializable {
 			year = ModelConfig.BASE_YEAR;
 
 		double bioenergyDemand = 0.0;
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			
 				double newDemand = bioenergyDemandManager.getFirstGenBioenergyDemand(c, year, commodity);
 				bioenergyDemand += newDemand; // aggregate if composite country
@@ -102,7 +100,7 @@ public abstract class AbstractDemandManager implements Serializable {
 			year = ModelConfig.BASE_YEAR;
 		
 		double bioenergyDemand = 0.0;
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			
 				double newDemand = bioenergyDemandManager.getSecondGenBioenergyDemand(c, year);
 				bioenergyDemand += newDemand; // aggregate if composite country
@@ -114,16 +112,15 @@ public abstract class AbstractDemandManager implements Serializable {
 
 	protected abstract Map<WoodType, Double> getWoodDemand(SingleCountry country, int year);
 	
-	public Map<Integer, Map<WoodType, Double>> getWoodDemandComposite(CompositeCountry cc, int year) {
+	public Map<WoodType, Double> getWoodDemandComposite(CompositeCountry cc, int year) {
 		if (!ModelConfig.CHANGE_DEMAND_YEAR)
 			year = ModelConfig.BASE_YEAR;
-		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);
+			for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
+				singleDemandMap = getWoodDemand(c, year);
 				if (singleDemandMap.isEmpty()) {
 					LogWriter.printlnWarning("Wood demand not found for: " + c.getCountryName());
 					continue;
@@ -135,10 +132,7 @@ public abstract class AbstractDemandManager implements Serializable {
 					compositeDemandMap.put(w, totalDemand);				
 				}
 			}
-			annualDemandMap.put(i, compositeDemandMap);
-		}
-
-		return annualDemandMap;
+		return compositeDemandMap;
 	}
 	
 	public double getWorldCarbonDemand(int year) {
@@ -153,7 +147,7 @@ public abstract class AbstractDemandManager implements Serializable {
 
 		double totalDemand = 0;
 		
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			totalDemand += getCarbonDemand(c, year);
 		}
 		
diff --git a/src/ac/ed/lurg/demand/AbstractSSPDemandManager.java b/src/ac/ed/lurg/demand/AbstractSSPDemandManager.java
index 153d194fa13db7d8eb7d064c4218186362c77268..0c7e3b193382ccaf2359873934fd254ea367bf5e 100644
--- a/src/ac/ed/lurg/demand/AbstractSSPDemandManager.java
+++ b/src/ac/ed/lurg/demand/AbstractSSPDemandManager.java
@@ -6,7 +6,6 @@ import java.util.Map;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.AgriculturalGDPManager;
 import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
@@ -22,8 +21,8 @@ public abstract class AbstractSSPDemandManager extends AbstractDemandManager {
 	private transient SingleCountry usa;
 	private transient AgriculturalGDPManager agriculturalGDPManager;
 
-	public AbstractSSPDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager,CalorieManager calorieManager, CompositeCountryManager compositeCountryManager) {
-		super(compositeCountryManager, calorieManager);
+	public AbstractSSPDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager,CalorieManager calorieManager) {
+		super(calorieManager);
 		setup(ssp_scenario, baseConsumpManager);
 	}
 	
@@ -31,19 +30,19 @@ public abstract class AbstractSSPDemandManager extends AbstractDemandManager {
 		this.ssp_scenario = ssp_scenario;
 		this.baseConsumpManager = baseConsumpManager;
 		sspManager = new SspManager();
-		usa = CountryManager.getForName("United States of America");
+		usa = CountryManager.getInstance().getForCode("USA");
 		agriculturalGDPManager = new AgriculturalGDPManager();
 	}
 
-	protected void setup(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager, CompositeCountryManager compositeCountryManager) {
-		super.setup(compositeCountryManager, calorieManager);
+	protected void setup(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager) {
+		super.setup(calorieManager);
 		setup(ssp_scenario, baseConsumpManager);
 	}
 	
 	@Override
 	public void updateGdpLossesFromShock(CompositeCountry cc, int year, double shockMagnitude) {
 
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			SspData sd = sspManager.get(ssp_scenario, year, c);
 			LogWriter.println("sd data " + sd, 3);
 			if(sd != null) {
@@ -122,7 +121,7 @@ public abstract class AbstractSSPDemandManager extends AbstractDemandManager {
 			return 0.0;
 		}
 		double gdp = sd.getGdp();
-		double worldGDP = sspManager.getWorldGdp(ssp_scenario, year, compositeCountryManager.getAllSingleCountries());
+		double worldGDP = sspManager.getWorldGdp(ssp_scenario, year, CountryManager.getInstance().getAllSingleCountries());
 		double worldCarbonDemand = getWorldCarbonDemand(year);
 		double countryCarbonDemand = worldCarbonDemand * (gdp / worldGDP);
 		return countryCarbonDemand;
diff --git a/src/ac/ed/lurg/demand/BaseConsumpManager.java b/src/ac/ed/lurg/demand/BaseConsumpManager.java
index b768c74114f7fc4e88b3daf2b82b12e8b028aa39..d5ceedc94c758a86b2373f96e924bf4f1ac6ff7c 100644
--- a/src/ac/ed/lurg/demand/BaseConsumpManager.java
+++ b/src/ac/ed/lurg/demand/BaseConsumpManager.java
@@ -14,15 +14,10 @@ import ac.ed.lurg.utils.LogWriter;
 
 public class BaseConsumpManager {
 	private static final int COUNTRY_COL = 0; 
-	private static final int BASE_POP_COL = 1; 
-	private static final int COMMODITY_COL = 2; 
-	private static final int BASE_CPC_COL = 3; 
-	private static final int POP_FOR_GROUPING_COL = 4; 
+	private static final int COMMODITY_COL = 1; 
+	private static final int BASE_CPC_COL = 2;
 
 	private Map<SingleCountry, Map<CommodityType, Double>> baseConsumpMap = new HashMap<SingleCountry, Map<CommodityType, Double>>();
-	// Don't like the duplication of population in this base consumption file, but guess not really a problem
-	private Map<SingleCountry, Double> basePop = new HashMap<SingleCountry, Double>();
-	private Map<SingleCountry, Double> popForGroupingMap = new HashMap<SingleCountry, Double>();
 
 	public BaseConsumpManager() {
 		super();
@@ -33,25 +28,23 @@ public class BaseConsumpManager {
 		String filename = ModelConfig.BASELINE_CONSUMP_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, commodityName;
-			double pop, cpc, popForGrouping;
+			String line, countryCode, commodityName;
+			double cpc;
 			fitReader.readLine(); // read header
 
 			while ((line=fitReader.readLine()) != null) {
 				String[] tokens = line.split(",");
 				
-				if (tokens.length < 4)
+				if (tokens.length < 3)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[COUNTRY_COL];
-				pop = Double.parseDouble(tokens[BASE_POP_COL]);
+				countryCode = tokens[COUNTRY_COL];
 				commodityName = tokens[COMMODITY_COL];
 				cpc = Double.parseDouble(tokens[BASE_CPC_COL]);
-				popForGrouping = Double.parseDouble(tokens[POP_FOR_GROUPING_COL]);
 
-				SingleCountry country = CountryManager.getForName(countryName);
+				SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 				if (country == null)
-					LogWriter.printlnError("Null country for" + countryName);
+					LogWriter.printlnError("Null country for " + countryCode);
 				else {
 					CommodityType commodity = CommodityType.getForFaoName(commodityName);
 					
@@ -62,8 +55,6 @@ public class BaseConsumpManager {
 					}
 					
 					commodityMap.put(commodity, cpc);
-					basePop.put(country, pop);
-					popForGroupingMap.put(country, popForGrouping);
 				}
 			} 
 			fitReader.close(); 
@@ -87,18 +78,5 @@ public class BaseConsumpManager {
 		Double d = commodityMap.get(commodity);
 		return d.doubleValue();
 	}
-	
-	public double getBasePopulation(SingleCountry country) {
-		Double d = basePop.get(country);
-		return (d != null) ? d.doubleValue() : Double.NaN;
-	}
-	
-	public double getPopForGrouping(SingleCountry country) {
-		Double d = popForGroupingMap.get(country);
-		return (d != null) ? d.doubleValue() : Double.NaN;
-	}
-	
-	public Collection<SingleCountry> getAllCountries() {
-		return basePop.keySet();
-	}
+
 }
diff --git a/src/ac/ed/lurg/demand/BioenergyDemandManager.java b/src/ac/ed/lurg/demand/BioenergyDemandManager.java
index 627186c5d741e295893c43ffadbf73606489db8a..0f0a3809fc07ca98250e19f485fabf5aaf61ea0d 100644
--- a/src/ac/ed/lurg/demand/BioenergyDemandManager.java
+++ b/src/ac/ed/lurg/demand/BioenergyDemandManager.java
@@ -16,8 +16,8 @@ import ac.ed.lurg.utils.LogWriter;
 
 public class BioenergyDemandManager {
 	private static final int GEN1_BASE_COUNTRY_COL = 0; 
-	private static final int GEN1_BASE_COMMODITY_COL = 2; 
-	private static final int GEN1_BASE_OTHER_COL = 3; 
+	private static final int GEN1_BASE_COMMODITY_COL = 1; 
+	private static final int GEN1_BASE_OTHER_COL = 2; 
 	private static final int GEN1_FUTURE_SCENARIO_COL = 0; 
 	private static final int GEN1_FUTURE_YEAR_COL = 3; 
 	private static final int GEN1_FUTURE_DEMAND_COL = 4;
@@ -46,21 +46,21 @@ public class BioenergyDemandManager {
 		String filename = ModelConfig.BIOENERGY_1GEN_BASE_DEMAND_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, commodityName;
+			String line, countryCode, commodityName;
 			Double other;
 			fitReader.readLine(); // read header
 
 			while ((line=fitReader.readLine()) != null) {
 				String[] tokens = line.split(",");
 				
-				if (tokens.length < 4)
+				if (tokens.length < 3)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[GEN1_BASE_COUNTRY_COL];
+				countryCode = tokens[GEN1_BASE_COUNTRY_COL];
 				commodityName = tokens[GEN1_BASE_COMMODITY_COL];
 				other = Double.valueOf(tokens[GEN1_BASE_OTHER_COL]);
 
-				SingleCountry country = CountryManager.getForName(countryName);
+				SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 				CommodityType crop = CommodityType.getForFaoName(commodityName);
 				
 				Map<CommodityType, Double> countryData = bioenergyMap.lazyGet(country);
@@ -129,7 +129,7 @@ public class BioenergyDemandManager {
 		String filename = ModelConfig.BIOENERGY_2GEN_DEMAND_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, model, scenario;
+			String line, countryCode, model, scenario;
 			Integer year;
 			Double demand;
 			fitReader.readLine(); // read header
@@ -140,14 +140,14 @@ public class BioenergyDemandManager {
 				if (tokens.length < 7)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[GEN2_COUNTRY_COL];
+				countryCode = tokens[GEN2_COUNTRY_COL];
 				model = tokens[GEN2_MODEL_COL];
 				scenario = tokens[GEN2_SCENARIO_COL];
 				year = Integer.valueOf(tokens[GEN2_YEAR_COL]);
 				demand = Double.valueOf(tokens[GEN2_DEMAND_COL]);
 				
 				if (model.equals(ModelConfig.GEN2_BIOENERGY_MODEL) && scenario.equals(ModelConfig.BIOENERGY_DEMAND_SCENARIO)) {
-					SingleCountry country = CountryManager.getForName(countryName);
+					SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 					Map<Integer, Double> countryData = bioenergyMap.lazyGet(country);
 					countryData.put(year, demand);
 				}
diff --git a/src/ac/ed/lurg/demand/CalorieManager.java b/src/ac/ed/lurg/demand/CalorieManager.java
index da54e6bf9ac04f425c17018ba7534be24e3afc24..66aab003c73ea1d6a007d422413d2e39a307a038 100644
--- a/src/ac/ed/lurg/demand/CalorieManager.java
+++ b/src/ac/ed/lurg/demand/CalorieManager.java
@@ -28,7 +28,7 @@ public class CalorieManager {
 		String filename = ModelConfig.CALORIE_PER_T_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, commodityName;
+			String line, countryCode, commodityName;
 			double kcalPerT;
 			fitReader.readLine(); // read header
 
@@ -38,14 +38,14 @@ public class CalorieManager {
 				if (tokens.length < 3)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[COUNTRY_COL];
+				countryCode= tokens[COUNTRY_COL];
 				commodityName = tokens[COMMODITY_COL];
 				kcalPerT = Double.parseDouble(tokens[CALORIE_PER_G_COL]);
 			
 
-				SingleCountry country = CountryManager.getForName(countryName);
+				SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 				if (country == null)
-					LogWriter.printlnError("Null country for" + countryName);
+					LogWriter.printlnError("Null country for " + countryCode);
 				else {
 					CommodityType commodity = CommodityType.getForFaoName(commodityName);
 					
diff --git a/src/ac/ed/lurg/demand/CerealFractionsManager.java b/src/ac/ed/lurg/demand/CerealFractionsManager.java
index 0f1496cf53070b74659645dacaea36a612fcd610..3f635f6d74a5c8711ced8ad189d2e5f1a8ee3f37 100644
--- a/src/ac/ed/lurg/demand/CerealFractionsManager.java
+++ b/src/ac/ed/lurg/demand/CerealFractionsManager.java
@@ -7,7 +7,7 @@ import java.util.Map.Entry;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
+import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
@@ -15,12 +15,10 @@ import ac.ed.lurg.utils.LogWriter;
 import ac.ed.lurg.utils.StringTabularReader;
 
 public class CerealFractionsManager {
-	protected CompositeCountryManager compositeCountryManager;
 	private StringTabularReader cerealFractionReader;
 
-	public CerealFractionsManager(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
-		cerealFractionReader = new StringTabularReader(",", new String[]{"Country", "plumDemandItem", "Item", "food"});
+	public CerealFractionsManager() {
+		cerealFractionReader = new StringTabularReader(",", new String[]{"Iso3", "plumDemandItem", "plumCropItem", "food"});
 		cerealFractionReader.read(ModelConfig.BASE_DEMAND_FRACT_FILE);
 	}
 
@@ -48,7 +46,7 @@ public class CerealFractionsManager {
 	private Map<CropType, Double> getBaseDemandFracts(CompositeCountry cc, CommodityType commodity) {
 
 		Map<CropType, Double> commodityDemandMap = new HashMap<CropType, Double>();
-		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+		for (SingleCountry c : CountryManager.getInstance().getAllForCompositeCountry(cc)) {
 			Map<CropType, Double> singleCMap = getCommodityFraction(c, commodity);
 
 			for (CropType crop : commodity.getCropTypes()) {
@@ -81,8 +79,8 @@ public class CerealFractionsManager {
 
 		for (CropType crop : commodity.getCropTypes()) {
 			Map<String, String> queryMap = new HashMap<String, String>();
-			queryMap.put("Country", c.getCountryName());
-			queryMap.put("Item", crop.getFaoName());
+			queryMap.put("Iso3", c.getCountryCode());
+			queryMap.put("plumCropItem", crop.getFaoName());
 			queryMap.put("plumDemandItem", commodity.getFaoName());
 
 			try {
diff --git a/src/ac/ed/lurg/demand/DemandManagerFromFile.java b/src/ac/ed/lurg/demand/DemandManagerFromFile.java
index a4c63b6e44708387acc70ff4bc3798d95a64de47..4a97829dba7372de18acc842d38ea2d36d579a99 100644
--- a/src/ac/ed/lurg/demand/DemandManagerFromFile.java
+++ b/src/ac/ed/lurg/demand/DemandManagerFromFile.java
@@ -5,7 +5,6 @@ import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.WoodType;
@@ -18,12 +17,12 @@ public class DemandManagerFromFile extends AbstractDemandManager {
 	private StringTabularReader historicalConsumptionReader;
 	private StringTabularReader histWoodConsumpReader;
 
-	public DemandManagerFromFile(CompositeCountryManager compositeCountryManager, CalorieManager calorieManager) {
-		super(compositeCountryManager, calorieManager);
-		historicalConsumptionReader = new StringTabularReader(",", new String[]{"Year", "Country", "Item", "population", "cpc"});
+	public DemandManagerFromFile(CalorieManager calorieManager) {
+		super(calorieManager);
+		historicalConsumptionReader = new StringTabularReader(",", new String[]{"Year", "Iso3", "Item", "population", "cpc"});
 		historicalConsumptionReader.read(ModelConfig.DEMAND_CONSUMPTION_FILE);
 		
-		histWoodConsumpReader = new StringTabularReader(",", new String[]{"country", "year", "type", "demand_m3", "net_imports", "harvest"});
+		histWoodConsumpReader = new StringTabularReader(",", new String[]{"iso3", "year", "type", "demand_m3", "net_imports", "harvest"});
 		histWoodConsumpReader.read(ModelConfig.WOOD_DEMAND_FILE);
 	}
 
@@ -34,7 +33,7 @@ public class DemandManagerFromFile extends AbstractDemandManager {
 		for (CommodityType commodity : CommodityType.getAllFoodItems()) {
 			Map<String, String> queryMap = new HashMap<String, String>();
 			queryMap.put("Year", Integer.toString(year));
-			queryMap.put("Country", c.getCountryName());
+			queryMap.put("Iso3", c.getCountryCode());
 			queryMap.put("Item", commodity.getFaoName());
 			try {
 				Map<String, String> row = historicalConsumptionReader.querySingle(queryMap);
@@ -72,7 +71,7 @@ public class DemandManagerFromFile extends AbstractDemandManager {
 			
 			Map<String, String> queryMap = new HashMap<String, String>();
 			queryMap.put("year", Integer.toString(year));
-			queryMap.put("country", country.getCountryName());
+			queryMap.put("iso3", country.getCountryName());
 			queryMap.put("type", woodType.getName());
 			try {
 				Map<String, String> row = histWoodConsumpReader.querySingle(queryMap);
diff --git a/src/ac/ed/lurg/demand/DemandManagerSSP.java b/src/ac/ed/lurg/demand/DemandManagerSSP.java
index 67c37e4168521eb491749221d6d4a838a9597b8a..902a2840231eb19afd1d6253a3ca4b1be8dd3e1d 100644
--- a/src/ac/ed/lurg/demand/DemandManagerSSP.java
+++ b/src/ac/ed/lurg/demand/DemandManagerSSP.java
@@ -4,7 +4,6 @@ import java.util.HashMap;
 import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.ModelFitType;
@@ -15,8 +14,8 @@ public class DemandManagerSSP extends AbstractSSPDemandManager {
 	private static final long serialVersionUID = 6500619362109439516L;
 	private DemandCurveManager demandCurveManager;
 
-	public DemandManagerSSP(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager, CompositeCountryManager compositeCountryManager) {
-		super(ssp_scenario, baseConsumpManager, calorieManager, compositeCountryManager);
+	public DemandManagerSSP(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager) {
+		super(ssp_scenario, baseConsumpManager, calorieManager);
 		demandCurveManager = new DemandCurveManager();		
 	}
 
diff --git a/src/ac/ed/lurg/demand/ElasticDemandManager.java b/src/ac/ed/lurg/demand/ElasticDemandManager.java
index 7d53594ab0bf779fc5366ce7207a8571cc194dd9..aa3e634e603f74e4576b8dd53349cd526f28eb5d 100755
--- a/src/ac/ed/lurg/demand/ElasticDemandManager.java
+++ b/src/ac/ed/lurg/demand/ElasticDemandManager.java
@@ -6,7 +6,6 @@ import java.util.HashMap;
 import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.country.gams.GamsCommodityDemand;
 import ac.ed.lurg.country.gams.GamsDemandInput;
@@ -27,8 +26,8 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 	private int baseYearToRebaseConsumption;
 
 	/** Constructor used in calibration only */
-	public ElasticDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager, CompositeCountryManager compositeCountryManager) {
-		super(ssp_scenario, baseConsumpManager, calorieManager, compositeCountryManager);
+	public ElasticDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager) {
+		super(ssp_scenario, baseConsumpManager, calorieManager);
 
 		this.priceMarkups = new PriceMarkUps();
 		baseCpcCache = new HashMap<SingleCountry, Map<CommodityType, Double>>();		
@@ -47,8 +46,8 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 		}
 	}
 
-	public void setup(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager, CompositeCountryManager compositeCountryManager) {
-		super.setup(ssp_scenario, baseConsumpManager, calorieManager, compositeCountryManager);
+	public void setup(String ssp_scenario, BaseConsumpManager baseConsumpManager, CalorieManager calorieManager) {
+		super.setup(ssp_scenario, baseConsumpManager, calorieManager);
 		setup();
 	}
 	
@@ -64,7 +63,9 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 		// convert from producer to consumer prices
 		for (CommodityType commodity : CommodityType.getAllFoodItems()) {
 			double producerPricePlum = producerPrices.get(commodity);
-			double producerPrice = producerPricePlum * (1000/kcalPerT.get(commodity) * 2000 *365);  //price per person per year commodity x per year assuming 2000 kcal a day
+			// price per person per year assuming 2000 kcal a day (in $ per year)
+			// $1000 * $1000/t commodity * (1 / kcalPerT) * 2000 kcal per day * 365 days
+			double producerPrice = 1000 * producerPricePlum * (1 / kcalPerT.get(commodity)) * 2000 * 365; 
 
 			double consumerPrice = priceMarkups.markedupPrice(c, commodity, producerPrice);
 			consumerPrices.put(commodity, consumerPrice);
@@ -164,9 +165,13 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 	public Map<WoodType, Double> getWoodDemandMap(SingleCountry country, double gdpPc, double baseGdpPc, double population, double basePopulation) {
 		Map<WoodType, Double> demandMap = new HashMap<WoodType, Double>();
 		for (WoodType woodType : WoodType.values()) {
-			double baseDemand = woodDemandManager.get(country, woodType);
-			double newDemand = baseDemand * (population / basePopulation) * Math.pow(gdpPc / baseGdpPc, woodType.getIncomeDemandElasticity());
-			demandMap.put(woodType, newDemand);
+			if (!ModelConfig.IS_FORESTRY_ON) {
+				demandMap.put(woodType, 0.0);
+			} else {
+				double baseDemand = woodDemandManager.get(country, woodType);
+				double newDemand = baseDemand * (population / basePopulation) * Math.pow(gdpPc / baseGdpPc, woodType.getIncomeDemandElasticity());
+				demandMap.put(woodType, newDemand);
+			}
 		}
 		return demandMap;
 	}
diff --git a/src/ac/ed/lurg/demand/PriceMarkUps.java b/src/ac/ed/lurg/demand/PriceMarkUps.java
index e8ab6de9e6f9917a923935ac2cf84f3973042696..a011be149c092d31e71548a8ee44c5cef2c21ed0 100644
--- a/src/ac/ed/lurg/demand/PriceMarkUps.java
+++ b/src/ac/ed/lurg/demand/PriceMarkUps.java
@@ -16,7 +16,7 @@ public class PriceMarkUps {
 	public PriceMarkUps() {
 		// could read priceMarkupFactors from file here in some cases
 		
-		intialConsumerPrices = new StringTabularReader(",", new String[]{"Item", "country", "price2010"});
+		intialConsumerPrices = new StringTabularReader(",", new String[]{"Iso3", "Item", "Price"});
 		intialConsumerPrices.read(ModelConfig.INITIAL_CONSUMER_PRICE_FILE);
 	}
 	
@@ -39,12 +39,12 @@ public class PriceMarkUps {
 		
 		double initialPrice;
 		Map<String, String> queryMap = new HashMap<String, String>();
-		queryMap.put("country", country.getCountryName());
+		queryMap.put("Iso3", country.getCountryCode());
 		queryMap.put("Item", commodity.getGamsName());
 
 		try {
 			Map<String, String> row = intialConsumerPrices.querySingle(queryMap);
-			String priceS = row.get("price2010");
+			String priceS = row.get("Price");
 			initialPrice = Double.valueOf(priceS);
 		} catch (Exception e) {
 			LogWriter.println("Problem finding initial consumer price " + commodity.getGamsName() + ", " + commodity.getGamsName() + ", " + country.getCountryName());
diff --git a/src/ac/ed/lurg/demand/SspManager.java b/src/ac/ed/lurg/demand/SspManager.java
index 6f22b54b567c1d5939736c40eddff7af750deb59..b766b468ee84cbe21d19141d317759b60fad0eca 100644
--- a/src/ac/ed/lurg/demand/SspManager.java
+++ b/src/ac/ed/lurg/demand/SspManager.java
@@ -13,11 +13,11 @@ import ac.ed.lurg.utils.LogWriter;
 
 public class SspManager {
 //	private static final int MODEL_COL = 0; 
+	private static final int COUNTRY_COL = 0; 
 	private static final int SCENARIO_COL = 1; 
 	private static final int YEAR_COL = 2; 
 	private static final int POPULATION_COL = 3; 
-	private static final int COUNTRY_COL = 4; 
-	private static final int GDP_COL = 5; 
+	private static final int GDP_COL = 4; 
 
 	private Map<SspKey, SspData> sspMap = new HashMap<SspKey, SspData>();
 
@@ -30,7 +30,7 @@ public class SspManager {
 		String filename = ModelConfig.SSP_FILE;
 		try {
 			BufferedReader reader = new BufferedReader(new FileReader(filename)); 
-			String line, scenario, countryName;
+			String line, scenario, countryCode;
 			double population, gdp;
 			int year;
 			reader.readLine(); // read header
@@ -38,7 +38,7 @@ public class SspManager {
 			while ((line=reader.readLine()) != null) {
 				String[] tokens = line.split(",");
 				
-				if (tokens.length < 6)
+				if (tokens.length < 5)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
 				try {
@@ -46,10 +46,10 @@ public class SspManager {
 					scenario = tokens[SCENARIO_COL];
 					year = Integer.parseInt(tokens[YEAR_COL]);
 					population = Double.parseDouble(tokens[POPULATION_COL]);
-					countryName = tokens[COUNTRY_COL];
+					countryCode = tokens[COUNTRY_COL];
 					gdp = Double.parseDouble(tokens[GDP_COL]);
 	
-					SspKey sspKey = new SspKey(scenario, year, CountryManager.getForName(countryName));
+					SspKey sspKey = new SspKey(scenario, year, CountryManager.getInstance().getForCode(countryCode));
 					SspData sspData = new SspData(population, gdp);
 					sspMap.put(sspKey, sspData);
 				} catch (Exception e) {
diff --git a/src/ac/ed/lurg/demand/WoodDemandManager.java b/src/ac/ed/lurg/demand/WoodDemandManager.java
index 35d237af92501c72931ee0eb1d64e5def3092756..d9fcac3f7325f05f038d7904e4fc666f826f04e9 100644
--- a/src/ac/ed/lurg/demand/WoodDemandManager.java
+++ b/src/ac/ed/lurg/demand/WoodDemandManager.java
@@ -1,9 +1,7 @@
 package ac.ed.lurg.demand;
 
-import java.io.BufferedReader;
-import java.io.FileReader;
-import java.io.IOException;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
@@ -11,68 +9,30 @@ import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.WoodType;
 import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.utils.StringTabularReader;
 
 public class WoodDemandManager {
-	private static final int COUNTY_COL = 0;
-	private static final int YEAR_COL = 1;
-	private static final int WOOD_TYPE_COL = 2; 
-	private static final int DEMAND_COL = 3; 
-	
 	private Map<SingleCountry, Map<WoodType, Double>> baseWoodDemand = new HashMap<SingleCountry, Map<WoodType, Double>>();
 	
 	public WoodDemandManager() {
-		readWoodDemandData();
+		readWoodDemand();
 	}
 	
-	public void readWoodDemandData() {
-		
+	private void readWoodDemand() {
 		String filename = ModelConfig.WOOD_DEMAND_FILE;
-		try {
-			BufferedReader reader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, woodTypeName;
-			Double demand;
-			Integer year;
-			reader.readLine(); // read header
-
-			while ((line=reader.readLine()) != null) {
-				String[] tokens = line.split(",");
-				
-				if (tokens.length < 6)
-					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
-				
-				year = Integer.valueOf(tokens[YEAR_COL]);
-				
-				if (year != ModelConfig.BASE_YEAR)
-					continue;
-				
-				countryName = tokens[COUNTY_COL];
-				woodTypeName = tokens[WOOD_TYPE_COL];
-				demand = Double.valueOf(tokens[DEMAND_COL]);
-				
-
-				SingleCountry country = CountryManager.getForName(countryName);
-				if (country == null)
-					LogWriter.printlnError("Null country for" + countryName);
-				else {
-					WoodType woodType = WoodType.getForName(woodTypeName);
-					
-					Map<WoodType, Double> woodTypeMap = baseWoodDemand.get(country);
-					if (woodTypeMap == null) {
-						woodTypeMap = new HashMap<WoodType, Double>();
-						baseWoodDemand.put(country, woodTypeMap);
-					}
-					
-					double demandBiomassEq = demand * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // converting from m3 to MtC-eq
-					woodTypeMap.put(woodType, demandBiomassEq);		
-				}
-			} 
-			reader.close(); 
-		
-		} catch (IOException e) {
-			LogWriter.printlnError("Failed in reading wood demand data");
-			LogWriter.print(e);
+		StringTabularReader reader = new StringTabularReader(",", new String[] {"Iso3", "Type", "Production", "Demand", "NetImports"});
+		reader.read(filename);
+		List<Map<String, String>> rowList = reader.getRowList();
+		for (Map<String, String> row : rowList) {
+			String countryCode = row.get("Iso3");
+			String woodTypeStr = row.get("Type");
+			double demand = Double.parseDouble(row.get("Demand"));
+			demand *= ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR; // converting from m3 to MtC-eq
+			SingleCountry sc = CountryManager.getInstance().getForCode(countryCode);
+			WoodType woodType = WoodType.getForName(woodTypeStr);
+			Map<WoodType, Double> wMap = baseWoodDemand.computeIfAbsent(sc, k -> new HashMap<WoodType, Double>());
+			wMap.put(woodType, demand);
 		}
-		LogWriter.println("Processed " + filename);
 	}
 	
 	public double get(SingleCountry country, WoodType woodType) {
@@ -83,7 +43,7 @@ public class WoodDemandManager {
 		if (woodTypeMap != null && woodTypeMap.containsKey(woodType)) {
 			return woodTypeMap.get(woodType).doubleValue();
 		} else {
-			LogWriter.printlnError("BaseConsumpManager: can't get value for " + woodType + ", " + country);
+			LogWriter.printlnError("WoodDemandManager: can't get value for " + woodType + ", " + country);
 			return 0.0;
 		}
 	}
diff --git a/src/ac/ed/lurg/forestry/ForestManager.java b/src/ac/ed/lurg/forestry/ForestManager.java
new file mode 100644
index 0000000000000000000000000000000000000000..0f8ccdd760637d06fed46bd2e602174adef7a3b3
--- /dev/null
+++ b/src/ac/ed/lurg/forestry/ForestManager.java
@@ -0,0 +1,153 @@
+package ac.ed.lurg.forestry;
+
+import ac.ed.lurg.country.LandCoverChangeItem;
+import ac.ed.lurg.country.gams.GamsLocationOutput;
+import ac.ed.lurg.landuse.LandCoverTile;
+import ac.ed.lurg.landuse.LandUseItem;
+import ac.ed.lurg.types.LandCoverType;
+import ac.ed.lurg.types.LandProtectionType;
+import ac.ed.lurg.types.WoodType;
+import ac.ed.lurg.utils.LogWriter;
+import ac.sac.raster.IntegerRasterItem;
+import ac.sac.raster.RasterKey;
+import ac.sac.raster.RasterSet;
+
+import java.util.*;
+
+public class ForestManager {
+
+    // Calculates amount of wood harvested from land cover change
+    public static Map<WoodType, Double> handleLccWood(RasterSet<WoodYieldItem> woodYieldRaster,
+                                                       Map<RasterKey, List<LandCoverChangeItem>> landCoverChanges,
+                                                       RasterSet<LandUseItem> inputLandUseRaster) {
+
+        Map<WoodType, Double> woodAmount = new HashMap<>();
+        for (Map.Entry<RasterKey, List<LandCoverChangeItem>> entry : landCoverChanges.entrySet()) {
+            RasterKey key = entry.getKey();
+            List<LandCoverChangeItem> changes = entry.getValue();
+            WoodYieldItem wyItem = woodYieldRaster.get(key);
+            LandCoverTile tiles = inputLandUseRaster.get(key).getLandCoverTiles();
+            for (LandCoverChangeItem lccItem : changes) {
+                LandCoverType fromLc = lccItem.getFromLandCover();
+                LandCoverType toLc = lccItem.getToLandCover();
+                double changedArea = lccItem.getArea();
+                double totalArea = tiles.getTotalArea(fromLc, LandProtectionType.CONVERTIBLE);
+                if (changedArea > totalArea) {
+                    LogWriter.printlnError("Error: ForestManager changed area exceeds total area.");
+                }
+                double fraction = changedArea / totalArea;
+                // Assume stands are converted proportionally by area.
+                // Type of wood harvested depends on stand age.
+                for (int age : tiles.getAgeKeys()) {
+                    double standArea = tiles.getArea(fromLc, LandProtectionType.CONVERTIBLE, age);
+                    double standYield = wyItem.getYieldLuc(fromLc, toLc, age);
+                    double standAmount = standYield * standArea * fraction;
+                    if (age > WoodType.IND_ROUNDWOOD.getMinRotation()) {
+                        woodAmount.merge(WoodType.IND_ROUNDWOOD, standAmount, Double::sum);
+                    } else {
+                        woodAmount.merge(WoodType.FUELWOOD, standAmount, Double::sum);
+                    }
+                }
+            }
+        }
+
+        return woodAmount;
+    }
+
+    // Calculates area of timber forest cleared to meet demand
+    public static Map<RasterKey, List<LandCoverChangeItem>> handleRotaWood(RasterSet<LandUseItem> outputLandUseRaster,
+                                                                            Map<WoodType, Double> woodDemand,
+                                                                           RasterSet<WoodYieldItem> woodYieldRaster,
+                                                                            GamsLocationOutput locationOutput,
+                                                                            RasterSet<IntegerRasterItem> mapping) {
+
+        final List<WoodType> harvestOrder = new ArrayList<>(Arrays.asList(WoodType.IND_ROUNDWOOD, WoodType.FUELWOOD));
+
+        Map<RasterKey, List<LandCoverChangeItem>> allChanges = new HashMap<>();
+
+        for (WoodType wType : harvestOrder) {
+            Map<RasterKey, List<LandCoverChangeItem>> changes = handleRotaWoodForType(outputLandUseRaster, woodDemand,
+                    woodYieldRaster, locationOutput, mapping, wType);
+            for (RasterKey key : changes.keySet()) {
+                List<LandCoverChangeItem> changesList = changes.get(key);
+                List<LandCoverChangeItem> changesListAll = allChanges.computeIfAbsent(key, k -> new ArrayList<>());
+                changesListAll.addAll(changesList);
+            }
+        }
+
+        return  allChanges;
+    }
+
+    private static Map<RasterKey, List<LandCoverChangeItem>> handleRotaWoodForType(RasterSet<LandUseItem> outputLandUseRaster,
+                                                                                   Map<WoodType, Double> woodDemand,
+                                                                                   RasterSet<WoodYieldItem> woodYieldRaster,
+                                                                                   GamsLocationOutput locationOutput,
+                                                                                   RasterSet<IntegerRasterItem> mapping,
+                                                                                   WoodType wType) {
+        Map<Integer, Map<WoodType, Double>> rotationPeriods = locationOutput.getRotationPeriods();
+        Map<RasterKey, List<LandCoverChangeItem>> lcChanges = new HashMap<>();
+
+            double demand = woodDemand.get(wType);
+            if (demand == 0)
+                return lcChanges;
+            Map<RasterKey, Double> stockMap = new HashMap<>();
+            for (RasterKey key : outputLandUseRaster.keySet()) {
+                LandUseItem luItem = outputLandUseRaster.get(key);
+                LandCoverTile tiles = luItem.getLandCoverTiles();
+                WoodYieldItem wyItem = woodYieldRaster.get(key);
+                Integer locId = mapping.get(key).getInt();
+                double stock = 0;
+                double rotationAge = 1 / rotationPeriods.get(locId).get(wType);
+                for (int age : tiles.getAgeKeys()) {
+                    if (age >= rotationAge) {
+                        double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
+                        double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
+                        stock += yield * area;
+                    }
+                }
+                stockMap.put(key, stock);
+            }
+            double totalStock = stockMap.values().stream().reduce(0.0, Double::sum);
+            double harvestFraction = demand / totalStock;
+            for (RasterKey key : outputLandUseRaster.keySet()) {
+                LandUseItem luItem = outputLandUseRaster.get(key);
+                LandCoverTile tiles = luItem.getLandCoverTiles();
+                WoodYieldItem wyItem = woodYieldRaster.get(key);
+                Integer locId = mapping.get(key).getInt();
+                double rotationAge = 1 / rotationPeriods.get(locId).get(wType);
+                double totalAmountToRemove = stockMap.get(key) * harvestFraction;
+                List<Integer> ageList = new ArrayList<>(tiles.getAgeKeys());
+                List<Integer> rotationAgeList = new ArrayList<>();
+                for (int i : ageList) {
+                    if (i >= rotationAge) {
+                        rotationAgeList.add(i);
+                    }
+                }
+                rotationAgeList.sort(Collections.reverseOrder());
+
+                for (int age : rotationAgeList) {
+                    if (totalAmountToRemove < 1e-10)
+                        break;
+                    if (age < rotationAge)
+                        continue;
+                    double area = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
+                    double yield = wyItem.getYieldRota(LandCoverType.TIMBER_FOREST, age);
+                    if (area == 0 || yield == 0)
+                        continue;
+                    double standStock = area * yield;
+                    double amountToRemove = Math.min(totalAmountToRemove, standStock);
+                    double previousArea = tiles.getArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age);
+                    double convertedArea = Math.min(previousArea, amountToRemove / yield);
+                    LandCoverChangeItem lccItem = new LandCoverChangeItem(LandCoverType.TIMBER_FOREST,
+                            LandCoverType.TIMBER_FOREST, age, convertedArea);
+                    List<LandCoverChangeItem> mapForKey = lcChanges.computeIfAbsent(key, k -> new ArrayList<>());
+                    mapForKey.add(lccItem);
+                    luItem.doForestRotation(age, convertedArea);
+                    totalAmountToRemove -= amountToRemove;
+                }
+            }
+
+        return lcChanges;
+    }
+
+}
diff --git a/src/ac/ed/lurg/forestry/ForestryPlanner.java b/src/ac/ed/lurg/forestry/ForestryPlanner.java
deleted file mode 100644
index 7de2c87cbf8cabed5998aa26a18cf423bf7eb914..0000000000000000000000000000000000000000
--- a/src/ac/ed/lurg/forestry/ForestryPlanner.java
+++ /dev/null
@@ -1,7 +0,0 @@
-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 e3a504c10cfa79c90134aee10e3e570ffd54f2f7..baaef2c61361f96be88d3e0bf2033fe8bcc33c6b 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldData.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldData.java
@@ -5,13 +5,14 @@ import java.util.Map;
 
 import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
+import ac.ed.lurg.types.WoodType;
 
 public class WoodYieldData {
-	private Map<LccKey, Double> yieldLuc = new HashMap<LccKey, Double>();
+	private Map<LccKey, Double> yieldLuc = new HashMap<>();
 	private WoodYieldResponse yieldResponse = new WoodYieldResponse();
 	
-	public void setYieldLuc(LccKey key, double yield) {
-		yieldLuc.put(key, yield);
+	public void setYieldLuc(LccKey lccKey, double yield) {
+		yieldLuc.put(lccKey, yield);
 	}
 	
 	public double getYieldLuc(LandCoverType fromLc, LandCoverType toLc) {
@@ -19,8 +20,8 @@ public class WoodYieldData {
 		return getYieldLuc(key);
 	}
 	
-	public double getYieldLuc(LccKey key) {
-		return yieldLuc.getOrDefault(key, 0.0);
+	public double getYieldLuc(LccKey lccKey) {
+		return yieldLuc.getOrDefault(lccKey, Double.NaN);
 	}
 	
 	public Map<LccKey, Double> getYieldLucMap() {
diff --git a/src/ac/ed/lurg/forestry/WoodYieldItem.java b/src/ac/ed/lurg/forestry/WoodYieldItem.java
index aa0aea7c56b16a7c0f57a76f7114750d49e68ce9..55266f5bb2829f5dae14da46a7c29d04c2417a33 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldItem.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldItem.java
@@ -4,8 +4,11 @@ import java.util.HashMap;
 import java.util.Map;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.landuse.LandCoverTile;
+import ac.ed.lurg.landuse.LandUseItem;
 import ac.ed.lurg.landuse.LccKey;
 import ac.ed.lurg.types.LandCoverType;
+import ac.ed.lurg.types.LandProtectionType;
+import ac.ed.lurg.types.WoodType;
 import ac.sac.raster.RasterItem;
 
 public class WoodYieldItem implements RasterItem {	
@@ -13,23 +16,20 @@ public class WoodYieldItem implements RasterItem {
 	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 setYieldLuc(LccKey key, Double[] yields, Map<LandCoverType, LandCoverTile> landUseTiles, int maxAge) {
-		LandCoverTile tiles = landUseTiles.get(key.getFromLc());
+	public void setYieldLuc(LccKey key, Double[] yields, LandCoverTile landUseTiles, int maxAge) {
 		Map<Integer, Double> yieldMap = yieldLuc.computeIfAbsent(key, k -> new HashMap<Integer, Double>());
-		for (Integer age : tiles.getAgeKeys()) {
+		for (Integer age : landUseTiles.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);
+	public void setYieldRota(LandCoverType lcType, Double[] yields, LandCoverTile landUseTiles, int maxAge) {
 		Map<Integer, Double> yieldMap = yieldRota.computeIfAbsent(lcType, k -> new HashMap<Integer, Double>());
-		for (Integer age : tiles.getAgeKeys()) {
+		for (Integer age : landUseTiles.getAgeKeys()) {
 			int ageCapped = Math.min(age, maxAge);
 			yieldMap.put(ageCapped, yields[ageCapped]);
 		}
@@ -37,7 +37,6 @@ public class WoodYieldItem implements RasterItem {
 		// For clustering
 		if (lcType.equals(LandCoverType.TIMBER_FOREST)) {
 			maxYield = yields[ModelConfig.MAX_FOREST_ROTATION_PERIOD];
-			potentialYield = yields[20] / 20;
 		}
 	}
 	
@@ -81,6 +80,60 @@ public class WoodYieldItem implements RasterItem {
 	public Map<LccKey, Map<Integer, Double>> getYieldLucMap() {
 		return yieldLuc;
 	}
+
+	// Average yield assuming proportional conversion of forest stands
+	public double getAggregatedYieldLuc(LccKey key, LandUseItem luItem) {
+		LandCoverType fromLc = key.getFromLc();
+		LandCoverTile tile = luItem.getLandCoverTiles();
+		double area = 0;
+		double stock = 0;
+		for (int age : tile.getAgeKeys()) {
+			double a = tile.getArea(fromLc, LandProtectionType.CONVERTIBLE, age);
+			double yield = getYieldLuc(key, age);
+			stock += a * yield;
+			area += a;
+		}
+		return area > 0 ? stock / area : 0;
+	}
+
+	// Yield for each wood type
+	public Map<LccKey, Map<WoodType, Double>> getYieldLucByType(LandUseItem luItem) {
+		Map<LccKey, Map<WoodType, Double>> yieldByType = new HashMap<>();
+		Map<LccKey, Map<WoodType, Double>> stockByType = new HashMap<>();
+		Map<LccKey, Map<WoodType, Double>> areaByType = new HashMap<>();
+
+		for (LandCoverType fromLc : LandCoverType.values()) {
+			for (LandCoverType toLc : LandCoverType.values()) {
+				LccKey key = new LccKey(fromLc, toLc);
+				LandCoverTile tile = luItem.getLandCoverTiles();
+				for (int age : tile.getAgeKeys()) {
+					double area = tile.getArea(fromLc, LandProtectionType.CONVERTIBLE, age);
+					double yield = getYieldLuc(key, age);
+					Map<WoodType, Double> stockMap = stockByType.computeIfAbsent(key, k -> new HashMap<>());
+					Map<WoodType, Double> areaMap = stockByType.computeIfAbsent(key, k -> new HashMap<>());
+					if (age >= WoodType.IND_ROUNDWOOD.getMinRotation()) {
+						areaMap.merge(WoodType.IND_ROUNDWOOD, area, Double::sum);
+						stockMap.merge(WoodType.IND_ROUNDWOOD, yield * area, Double::sum);
+					} else {
+						areaMap.merge(WoodType.FUELWOOD, area, Double::sum);
+						stockMap.merge(WoodType.FUELWOOD, yield * area, Double::sum);
+					}
+				}
+			}
+		}
+
+		for (Map.Entry<LccKey, Map<WoodType, Double>> entry : stockByType.entrySet()) {
+			LccKey key = entry.getKey();
+			Map<WoodType, Double> stockMap = entry.getValue();
+			for (WoodType w : stockMap.keySet()) {
+				double area = areaByType.get(key).get(w);
+				double yield = area > 0 ? stockMap.get(w) / area : 0;
+				yieldByType.computeIfAbsent(key, k -> new HashMap<>()).put(w, yield);
+			}
+		}
+
+		return yieldByType;
+	}
 	
 	public Map<LandCoverType, Map<Integer, Double>> getYieldRotaMap() {
 		return yieldRota;
@@ -89,9 +142,9 @@ public class WoodYieldItem implements RasterItem {
 	public Map<Integer, Double> getYieldResponseMap() {
 		return yieldResponse.getYieldMap();
 	}
-	
-	public double getPotentialYield() {
-		return potentialYield;
+
+	public static WoodYieldItem getDefault() {
+		return new WoodYieldItem();
 	}
 	
 }
diff --git a/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java b/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
index ae4ab078ad3bfb1941f4b36a9bd83d547d9f9512..bfc8c905b8df9fb3e2413c2e06e064cc6b695136 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldRasterSet.java
@@ -1,9 +1,6 @@
 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,21 +26,4 @@ public class WoodYieldRasterSet extends RasterSet<WoodYieldItem> {
 		popSubsetForKeys(subsetForestGrowthRaster, keys);		
 		return subsetForestGrowthRaster;
 	}
-	
-	public static double getPotentialTimberProduction(RasterSet<WoodYieldItem> woodYieldData, RasterSet<LandUseItem> landUseData) {
-		double total = 0;
-		for (RasterKey key : landUseData.keySet()) {
-			LandUseItem luItem = landUseData.get(key);
-			WoodYieldItem wyItem  = woodYieldData.get(key);
-			if (luItem == null | wyItem == null)
-				continue;
-			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 6a9a40e1981d1186583f197d44db4ed4ad4e9950..53dee70a534ab757819b523831f1fbdd7c4236c2 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldReader.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldReader.java
@@ -18,15 +18,13 @@ 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;
 
 public class WoodYieldReader {
 	private static final double CONVERSION_FACTOR = 10.0 * ModelConfig.WOOD_YIELD_CALIB_FACTOR; // convert kgC/m2 to tC/ha with calib factor
-	private RasterHeaderDetails rasterProj;
+	private final RasterHeaderDetails rasterProj;
 	private int modelYear;
 	private Map<Integer, RasterKey> keyMap;
 	double woodPrice;
@@ -42,51 +40,58 @@ public class WoodYieldReader {
 		this.woodPrice = woodPrice;
 		long startTime = System.currentTimeMillis();
 			// Wood yield from rotation
-			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_FORS_FILENAME,
+					List.of(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))); 
+			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_ROTA_NTRL_FILENAME,
+					List.of(new LccKey(LandCoverType.NATURAL, LandCoverType.NATURAL)));
 			
 			// Wood yield from LCC
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_PAST_FILENAME, 
-					getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.PASTURE)));
+					getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.PASTURE)));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_CROP_FILENAME, 
-					getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.CROPLAND))); 
+					getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.CROPLAND)));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_NTRL_FILENAME, 
-					getLccMappings(LandCoverType.getManagedForestTypes(), Arrays.asList(LandCoverType.NATURAL)));
+					getLccMappings(LandCoverType.getManagedForestTypes(), List.of(LandCoverType.NATURAL)));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_FORS_TO_FORS_FILENAME, 
 					getLccMappings(LandCoverType.getManagedForestTypes(), LandCoverType.getManagedForestTypes()));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_PAST_FILENAME, 
-					getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.PASTURE))); 
+					getLccMappings(List.of(LandCoverType.NATURAL), List.of(LandCoverType.PASTURE)));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_CROP_FILENAME, 
-					getLccMappings(Arrays.asList(LandCoverType.NATURAL), Arrays.asList(LandCoverType.CROPLAND))); 
+					getLccMappings(List.of(LandCoverType.NATURAL), List.of(LandCoverType.CROPLAND)));
 			
 			readLccData(woodYieldData, landUseRaster, ModelConfig.WOOD_YIELD_NTRL_TO_FORS_FILENAME, 
-					getLccMappings(Arrays.asList(LandCoverType.NATURAL), LandCoverType.getManagedForestTypes())); 
+					getLccMappings(List.of(LandCoverType.NATURAL), LandCoverType.getManagedForestTypes()));
 			LogWriter.println("Reading forestry data took " + (System.currentTimeMillis() - startTime) + " ms");
 		return woodYieldData;
 	}
 	
+	/* Data schema
+	 * rows correspond to age class, columns correspond to locations
+	 * serial binary format. data[i + j * nLoc] to get value for location i, age class j, number of locations nLoc
+	 */
+	
 	public void readLccData(WoodYieldRasterSet woodYieldData, RasterSet<LandUseItem> landUseRaster, String fileName, List<LccKey> lccMappings) {
 		long startTime = System.currentTimeMillis();
 		int maxYear = Math.floorDiv(modelYear - 10, ModelConfig.CARBON_DATA_TIMESTEP_SIZE) * ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 10;
 		int maxAge = modelYear - ModelConfig.CARBON_DATA_MIN_YEAR;
-		Map<RasterKey, Double[]> map = new HashMap<RasterKey, Double[]>();
+		Map<RasterKey, Double[]> map = new HashMap<>();
 		for (RasterKey key : keyMap.values()) {
 			map.put(key, new Double[maxAge + 1]);
 		}
+		// We need to extract data from each establishment year group
 		for (int dataYear = maxYear; dataYear >= ModelConfig.CARBON_DATA_MIN_YEAR; dataYear -= ModelConfig.CARBON_DATA_TIMESTEP_SIZE) {
-			int nCol = keyMap.size();
-			int endAge = modelYear - dataYear;
-			int startAge = Math.max(0, endAge - ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 1);
+			int nCol = keyMap.size(); // number of locations
+			int endAge = modelYear - dataYear; // oldest age class we need in this dataset
+			int startAge = Math.max(0, endAge - ModelConfig.CARBON_DATA_TIMESTEP_SIZE + 1); // youngest age class
 			String filePath = getDataDir(fileName, dataYear); 
 			try {
+				// Read binary data.
 				RandomAccessFile reader = new RandomAccessFile(filePath, "r");
 				int startPos = nCol * startAge * Double.BYTES;
 				int nRead = (endAge - startAge + 1) * nCol * Double.BYTES;
@@ -124,19 +129,8 @@ public class WoodYieldReader {
 		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>();
+		Map<Integer, RasterKey> keyMap = new HashMap<>();
 		String filePath = ModelConfig.FORESTRY_DIR + File.separator + "coordinates.dat";
 		byte[] bytes;
 		try {
@@ -170,7 +164,7 @@ public class WoodYieldReader {
 	}
 	
 	private ArrayList<LccKey> getLccMappings(Collection<LandCoverType> fromTypes, Collection<LandCoverType> toTypes) {
-		ArrayList<LccKey> lccMappings = new ArrayList<LccKey>();
+		ArrayList<LccKey> lccMappings = new ArrayList<>();
 		for (LandCoverType fromLc : fromTypes) {
 			for (LandCoverType toLc : toTypes) {
 				if (!fromLc.equals(toLc)) // exclude no change as yield should be 0
diff --git a/src/ac/ed/lurg/forestry/WoodYieldResponse.java b/src/ac/ed/lurg/forestry/WoodYieldResponse.java
index e68902f0c55cd32619c572647f3ce1d32c4d7a5b..a956cc28e472628f5daaaaaab54279c51d8ba3a6 100644
--- a/src/ac/ed/lurg/forestry/WoodYieldResponse.java
+++ b/src/ac/ed/lurg/forestry/WoodYieldResponse.java
@@ -34,7 +34,7 @@ public class WoodYieldResponse {
 			kEst.add(Math.log(1 - y / maxYield) / (-a));
 		}
 		Collections.sort(kEst);
-		slope = kEst.get((int) kEst.size() / 2); 
+		slope = kEst.get(kEst.size() / 2);
 
 		double minA = maxYield * 0.8;
 		double stepA = (maxYield * 1.2 - minA) / 100;
diff --git a/src/ac/ed/lurg/landuse/ConversionCostReader.java b/src/ac/ed/lurg/landuse/ConversionCostReader.java
index 5b18af6bba3d4c13bfd111dbd46d16aa8d3ef323..4017beea3c784dbff7e808090d51b4461e13ad67 100644
--- a/src/ac/ed/lurg/landuse/ConversionCostReader.java
+++ b/src/ac/ed/lurg/landuse/ConversionCostReader.java
@@ -17,8 +17,7 @@ public class ConversionCostReader {
 	private static final int TO_COL = 1;
 	private static final int COST_COL = 2;
 	private Map<LccKey, Double> conversionCosts = new HashMap<LccKey, Double>();
-	
-	
+
 	
 	public ConversionCostReader(Timestep timestep) {
 		if (timestep.getTimestep() <= ModelConfig.END_SECOND_STAGE_CALIBRATION && ModelConfig.IS_CALIBRATION_RUN) {
diff --git a/src/ac/ed/lurg/landuse/CropUsageReader.java b/src/ac/ed/lurg/landuse/CropUsageReader.java
index dd2f2c024cecb1e63b3de37369158bae5ad7a56c..f2b9d85a9fa52d530dd51be4a70f9835cc99928e 100644
--- a/src/ac/ed/lurg/landuse/CropUsageReader.java
+++ b/src/ac/ed/lurg/landuse/CropUsageReader.java
@@ -8,7 +8,6 @@ import java.util.Map;
 
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.CompositeCountry;
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.types.CropType;
@@ -16,15 +15,12 @@ import ac.ed.lurg.utils.LazyHashMap;
 import ac.ed.lurg.utils.LogWriter;
 
 public class CropUsageReader {
-	private static final int COUNTRY_COL = 1; 
-	private static final int COMMODITY_COL = 0; 
-	private static final int NET_IMPORT_COL = 3; 
-	private static final int PRODUCTION_COL = 4; 
-
-	private CompositeCountryManager compositeCountryManager;
+	private static final int COUNTRY_COL = 0; 
+	private static final int COMMODITY_COL = 1; 
+	private static final int NET_IMPORT_COL = 2; 
+	private static final int PRODUCTION_COL = 3; 
 	
-	public CropUsageReader(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
+	public CropUsageReader() {
 	}
 
 	@SuppressWarnings("serial")
@@ -41,7 +37,7 @@ public class CropUsageReader {
 		String filename = ModelConfig.NET_IMPORTS_FILE;
 		try {
 			BufferedReader fitReader = new BufferedReader(new FileReader(filename)); 
-			String line, countryName, commodityName;
+			String line, countryCode, commodityName;
 			double netImports, production;
 			fitReader.readLine(); // read header
 
@@ -51,15 +47,15 @@ public class CropUsageReader {
 				if (tokens.length < 4)
 					LogWriter.printlnError("Too few columns in " + filename + ", " + line);
 				
-				countryName = tokens[COUNTRY_COL];
+				countryCode = tokens[COUNTRY_COL];
 				commodityName = tokens[COMMODITY_COL];
 				netImports = Double.parseDouble(tokens[NET_IMPORT_COL]);
 				production = Double.parseDouble(tokens[PRODUCTION_COL]);
 
-				SingleCountry country = CountryManager.getForName(countryName);
+				SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
 				
 				CropType crop = CropType.getForFaoName(commodityName);
-				CompositeCountry cc = compositeCountryManager.getForSingleCountry(country);
+				CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(country);
 				
 				Map<CropType, CropUsageData> countryData = commodityMap.lazyGet(cc);
 				CropUsageData oldData = countryData.get(crop);
diff --git a/src/ac/ed/lurg/landuse/FPUManager.java b/src/ac/ed/lurg/landuse/FPUManager.java
index 6ca11ee34c548ece8367fcbb49ba9da0d09f9eb3..5baf87c4bd5312f51817caeee4b58bf5285dec99 100644
--- a/src/ac/ed/lurg/landuse/FPUManager.java
+++ b/src/ac/ed/lurg/landuse/FPUManager.java
@@ -64,7 +64,7 @@ public class FPUManager {
 			}			
 		}
 	}
-	
+	// TODO shouldn't be relying on exceptions!
 	private WaterBasin getBasin(Fpu fpu) {
 		Map<String, String> queryMap = new HashMap<String, String>();
 		queryMap.put("fpu", Integer.toString(fpu.getId()));
diff --git a/src/ac/ed/lurg/landuse/IrrigationItem.java b/src/ac/ed/lurg/landuse/IrrigationItem.java
index 49dd5934bf21cb73dba99a4f03d4e414d2e23190..c650e75a0a407c2aba8b3038cf7e466d10ac715a 100644
--- a/src/ac/ed/lurg/landuse/IrrigationItem.java
+++ b/src/ac/ed/lurg/landuse/IrrigationItem.java
@@ -89,4 +89,15 @@ public class IrrigationItem implements RasterItem {
 	public void setConstraintOffsets(double constraintOffset) {
 		this.constraintOffset = constraintOffset;
 	}
+
+	public static IrrigationItem getDefault() {
+		IrrigationItem item = new IrrigationItem();
+		item.setIrrigCost(0.0);
+		item.setBaselineIrrigConstraint(0.0);
+		item.setBaseRunOff(0.0);
+		for (CropType crop: CropType.values()) {
+			item.setMaxIrrigAmount(crop, 0.0);
+		}
+		return item;
+	}
 }
diff --git a/src/ac/ed/lurg/landuse/LandCoverItem.java b/src/ac/ed/lurg/landuse/LandCoverItem.java
index cd292e713a504b7784b44159571a495f96b34921..1e0fce7942302f6aeb5144fc72cb49b3d8d28c1b 100644
--- a/src/ac/ed/lurg/landuse/LandCoverItem.java
+++ b/src/ac/ed/lurg/landuse/LandCoverItem.java
@@ -18,7 +18,7 @@ public class LandCoverItem implements RasterItem {
 	
 	private Map<LandCoverType, Double> landCoverFract = new HashMap<LandCoverType, Double>(); // Land cover fraction
 	private Map<LandCoverType, Map<Integer, Double>> landUseAgeDist = new HashMap<LandCoverType, Map<Integer, Double>>();
-	private Map<LandCoverType, LandCoverTile> landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
+	private LandCoverTile landCoverAreas = new LandCoverTile();
 	private double totalArea;
 	private double unavailableFract; // due to slope
 	private double protectedFraction;
@@ -78,7 +78,7 @@ public class LandCoverItem implements RasterItem {
 		this.protectedFraction = protectedFraction;
 	}
 	
-	public Map<LandCoverType, LandCoverTile> getLandCoverTiles() {
+	public LandCoverTile getLandCoverTiles() {
 		createLandCoverTiles();
 		setUnavailableArea();
 		cleanUp();
@@ -86,11 +86,6 @@ public class LandCoverItem implements RasterItem {
 	}
 	
 	public void createLandCoverTiles() {
-		
-		for (LandCoverType lcType : LandCoverType.values()) {
-			landCoverAreas.put(lcType, new LandCoverTile());
-		}
-		
 		double protectableFraction = 0;
 		for (LandCoverType lcType : LandCoverType.getProtectibleTypes()) {
 			protectableFraction += getLandCoverFract(lcType);
@@ -99,14 +94,13 @@ public class LandCoverItem implements RasterItem {
 		double adjProtectedFraction = Math.min(1.0, protectedFraction / protectableFraction); 
 				
 		for (LandCoverType lcType : LandCoverType.values()) {
-			LandCoverTile tiles = landCoverAreas.get(lcType);
 			double totalArea = getLandCoverArea(lcType);
 			if (totalArea == 0) 
 				continue;
 
 			// static land covers
 			if (!lcType.isConvertible()) { 
-				tiles.addArea(LandProtectionType.UNAVAILABLE, 0, totalArea); // Assumes land cover age is 0
+				landCoverAreas.addArea(lcType, LandProtectionType.UNAVAILABLE, 0, totalArea); // Assumes land cover age is 0
 				continue;
 			}
 			// convertible land covers			
@@ -149,10 +143,10 @@ public class LandCoverItem implements RasterItem {
 				for (int a = minAge; a <= maxAge; a+=2) {
 					double unprotectedA = unprotectedArea * prop * groupProp;
 					if (unprotectedA > 0)
-						landCoverAreas.get(lcType).setArea(LandProtectionType.CONVERTIBLE, a, unprotectedA);
+						landCoverAreas.setArea(lcType ,LandProtectionType.CONVERTIBLE, a, unprotectedA);
 					double protectedA = protectedArea * prop * groupProp;
 					if (protectedA > 0)
-						landCoverAreas.get(lcType).setArea(LandProtectionType.PROTECTED, a, protectedA);		
+						landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, a, protectedA);
 				}	
 			}
 		}
@@ -167,24 +161,27 @@ public class LandCoverItem implements RasterItem {
 	}
 	
 	private void setUnavailableArea() {
-		double alreadyUnavail = landCoverAreas.values().stream().mapToDouble(o -> o.getTotalArea(LandProtectionType.UNAVAILABLE)).sum();
+		double alreadyUnavail = 0;
+		for (LandCoverType lcType : LandCoverType.values()) {
+			alreadyUnavail += landCoverAreas.getTotalArea(lcType, LandProtectionType.UNAVAILABLE);
+		}
 		double unavailableArea = unavailableFract * totalArea;
 		double shortfall = Math.max(0, unavailableArea - alreadyUnavail);
 		if (shortfall == 0)
 			return;
 
 		for (LandCoverType lcType : unavailPriorityList) {
-			double cArea = landCoverAreas.get(lcType).getTotalArea(LandProtectionType.CONVERTIBLE);
-			double pArea = landCoverAreas.get(lcType).getTotalArea(LandProtectionType.PROTECTED);
+			double cArea = landCoverAreas.getTotalArea(lcType, LandProtectionType.CONVERTIBLE);
+			double pArea = landCoverAreas.getTotalArea(lcType, LandProtectionType.PROTECTED);
 			double totalArea = cArea + pArea;
 			double cProp = cArea / totalArea;
 			
 			double additionalUnavail = Math.min(shortfall, totalArea);
 			if (additionalUnavail > 0) {
 				// all additional unavailable area moved to NATURAL land cover
-				landCoverAreas.get(LandCoverType.NATURAL).addArea(LandProtectionType.UNAVAILABLE, 0, additionalUnavail);
-				landCoverAreas.get(lcType).removeArea(LandProtectionType.CONVERTIBLE, additionalUnavail * cProp);
-				landCoverAreas.get(lcType).removeArea(LandProtectionType.PROTECTED, additionalUnavail * (1 - cProp));
+				landCoverAreas.addArea(LandCoverType.NATURAL, LandProtectionType.UNAVAILABLE, 0, additionalUnavail);
+				landCoverAreas.removeArea(lcType, LandProtectionType.CONVERTIBLE, additionalUnavail * cProp);
+				landCoverAreas.removeArea(lcType, LandProtectionType.PROTECTED, additionalUnavail * (1 - cProp));
 				shortfall -= additionalUnavail;
 			}
 			
@@ -197,21 +194,20 @@ public class LandCoverItem implements RasterItem {
 		}
 	}
 	
-	public void runAreaCheck(Map<LandCoverType, Double> lcAreas, Map<LandCoverType, LandCoverTile> landCoverAreas) {
+	public void runAreaCheck(Map<LandCoverType, Double> lcAreas, LandCoverTile landCoverAreas) {
 		final double THRESHOLD = 1e-7;
-		for (LandCoverType lcType : landCoverAreas.keySet()) {
+		for (LandCoverType lcType : LandCoverType.values()) {
 			double initArea = lcAreas.get(lcType);
-			double processedArea = landCoverAreas.get(lcType).getTotalArea();
+			double processedArea = landCoverAreas.getTotalArea(lcType);
 			if (Math.abs(initArea - processedArea) > THRESHOLD) {
-				LogWriter.printlnError("LandUseItem: processed area and initial area different: " + lcType.getName() + " " + processedArea + ", " + initArea);
+				LogWriter.printlnError("LandUseItem: processed area and initial area different: "
+						+ lcType.getName() + " " + processedArea + ", " + initArea);
 			}
 		}
 	}
 	
 	private void cleanUp() {
-		for (LandCoverTile tile : landCoverAreas.values()) {
-			tile.cleanUp();
-		}
+		landCoverAreas.cleanUp();
 	}
 	
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/landuse/LandCoverTile.java b/src/ac/ed/lurg/landuse/LandCoverTile.java
index fef1a772453fe4c84c7dd5ca07b2bdd557f8a89d..246adc279a0c3fc8ebf686fe13d9db41d57f4eb1 100644
--- a/src/ac/ed/lurg/landuse/LandCoverTile.java
+++ b/src/ac/ed/lurg/landuse/LandCoverTile.java
@@ -7,6 +7,7 @@ import java.util.Map;
 import java.util.Set;
 
 import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.utils.LogWriter;
 import ac.sac.raster.InterpolatingRasterItem;
 import ac.ed.lurg.types.LandProtectionType;
@@ -18,8 +19,8 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 	
 	public LandCoverTile() {}
 
-	public double getArea(LandProtectionType lpType, int age) {
-		LandKey key = new LandKey(lpType, age);
+	public double getArea(LandCoverType lcType, LandProtectionType lpType, int age) {
+		LandKey key = new LandKey(lcType, lpType, age);
 		return getArea(key);
 	}
 	
@@ -27,8 +28,8 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 		return areas.getOrDefault(key, 0.0);
 	}
 	
-	public void setArea(LandProtectionType lpType, int age, double a) {
-		LandKey key = new LandKey(lpType, age);
+	public void setArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) {
+		LandKey key = new LandKey(lcType, lpType, age);
 		setArea(key, a);
 	}
 	
@@ -36,26 +37,28 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 		areas.put(key, a);
 	}
 	
-	public void addArea(LandProtectionType lpType, int age, double a) {
-		LandKey key = new LandKey(lpType, age);
+	public void addArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) {
+		LandKey key = new LandKey(lcType, lpType, age);
+		if (Double.isNaN(a)) {
+			int foo = 1;
+		}
 		addArea(key, a);
 	}
 	
 	public void addArea(LandKey key, double a) {
-		double prevArea = getArea(key);
-		areas.put(key, prevArea + a);
+		areas.merge(key, a, Double::sum);
 	}
 	
 	// This is used during calibration to preserve Land cover age distribution
-	public void addArea(LandProtectionType lpType, double a) {
-		double totalArea = getTotalArea(lpType);
+	public void addArea(LandCoverType lcType, LandProtectionType lpType, double a) {
+		double totalArea = getTotalArea(lcType, lpType);
 		if (totalArea == 0) {
-			addArea(lpType, 0, a); // assuming all age=0 if no previous land cover
+			addArea(lcType, lpType, 0, a); // assuming all age=0 if no previous land cover
 			
 		} else {
 			double factor = (a + totalArea) / totalArea; // proportional allocation
 			for (LandKey key : areas.keySet()) {
-				if (key.getLpType().equals(lpType)) {
+				if (key.getLpType().equals(lpType) && key.getLcType().equals(lcType)) {
 					double prevArea = getArea(key);
 					setArea(key, prevArea * factor);
 				}
@@ -63,19 +66,20 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 		}
 	}
 
-	public double getTotalArea(LandProtectionType lpType) {
-		double a = areas.keySet()
-				.stream()
-				.filter(k -> k.getLpType().equals(lpType))
-				.mapToDouble(k -> areas.get(k))
-				.sum();
+	public double getTotalArea(LandCoverType lcType, LandProtectionType lpType) {
+		double a = 0;
+		for (LandKey key : areas.keySet()) {
+			if (key.getLpType().equals(lpType) && key.getLcType().equals(lcType)) {
+				a += areas.get(key);
+			}
+		}
 		return a;
 	}
 	
-	public void removeArea(LandProtectionType lpType, int age, double a) {
+	public void removeArea(LandCoverType lcType, LandProtectionType lpType, int age, double a) {
 		if (a == 0)
 			return;
-		LandKey key = new LandKey(lpType, age);
+		LandKey key = new LandKey(lcType, lpType, age);
 		double prevArea = getArea(key);
 		if (a > prevArea) {
 			LogWriter.printlnError("LandCoverTile: attempting to remove too much area.");
@@ -84,35 +88,35 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 		setArea(key, prevArea - a);
 	}
 
-	public void removeArea(LandProtectionType lpType, double a) { // Subtracts area proportionally
+	public void removeArea(LandCoverType lcType, LandProtectionType lpType, double a) { // Subtracts area proportionally
 		if (a == 0)
 			return;
-		double totalArea = getTotalArea(lpType);
+		double totalArea = getTotalArea(lcType, lpType);
 		double areaToRemove = Math.min(a, totalArea);
 		if (Math.abs(a - areaToRemove) > 1e-10) {
 			LogWriter.printlnError("LandCoverTile: attempting to remove too much area.");
 		}
 		double factor = (totalArea - areaToRemove) / totalArea;
 		for (LandKey key : areas.keySet()) {
-			if (key.getLpType().equals(lpType)) {
+			if (key.getLcType().equals(lcType) && key.getLpType().equals(lpType)) {
 				double prevArea = getArea(key);
 				setArea(key, prevArea * factor);
 			}
 		}
 	}
 	
-	public double getTotalArea() {
+	public double getTotalArea(LandCoverType lcType) {
 		double area = 0;
 		for (LandProtectionType lpType : LandProtectionType.values()) {
-			area += getTotalArea(lpType);
+			area += getTotalArea(lcType, lpType);
 		}
 		return area;
 	}
 
-	public void removeAllArea(LandProtectionType lpType) {
+	public void removeAllArea(LandCoverType lcType, LandProtectionType lpType) {
 		Map<LandKey, Double> newAreas = new HashMap<LandKey, Double>();
 		for (LandKey key : areas.keySet()) {
-			if (!key.getLpType().equals(lpType)) {
+			if (!(key.getLpType().equals(lpType) && key.getLcType().equals(lcType))) {
 				newAreas.put(key, getArea(key));;
 			}
 		}
@@ -123,7 +127,7 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 	public void incrementAge() {
 		Map<LandKey, Double> newAreas = new HashMap<LandKey, Double>();
 		for (LandKey key : areas.keySet()) {
-			LandKey newKey = new LandKey(key.getLpType(), key.getAge() + ModelConfig.TIMESTEP_SIZE);
+			LandKey newKey = new LandKey(key.getLcType(), key.getLpType(), key.getAge() + ModelConfig.TIMESTEP_SIZE);
 			newAreas.put(newKey, getArea(key));
 		}
 		areas.clear();
@@ -131,18 +135,18 @@ public class LandCoverTile implements Serializable, InterpolatingRasterItem<Land
 	}
 
 	// Moves area between land protection types
-	public void moveArea(LandProtectionType fromPt, LandProtectionType toPt, double area) {
-		double fromA = getTotalArea(fromPt);
+	public void moveArea(LandCoverType lcType, LandProtectionType fromPt, LandProtectionType toPt, double area) {
+		double fromA = getTotalArea(lcType, fromPt);
 		if (area > fromA) {
 			LogWriter.printlnWarning("LandCoverTile.moveArea cannot protect all required area");
 			area = Math.min(fromA, area);
 		}
 		double fractToMove = area / fromA;
 		for (LandKey key : areas.keySet()) {
-			if (key.getLpType().equals(fromPt)) {
+			if (key.getLpType().equals(fromPt) & key.getLcType().equals(lcType)) {
 				double areaToMove = getArea(key) * fractToMove;
-				addArea(toPt, key.getAge(), areaToMove);
-				removeArea(fromPt, key.getAge(), areaToMove);
+				addArea(lcType, toPt, key.getAge(), areaToMove);
+				removeArea(lcType, fromPt, key.getAge(), areaToMove);
 			}
 		}
 	}
diff --git a/src/ac/ed/lurg/landuse/LandKey.java b/src/ac/ed/lurg/landuse/LandKey.java
index aae009aea20747fc517b8101cbe78c97f5cafa7d..287b3c54d185dcaabfb676346422b4da2085d3f7 100644
--- a/src/ac/ed/lurg/landuse/LandKey.java
+++ b/src/ac/ed/lurg/landuse/LandKey.java
@@ -3,20 +3,26 @@ package ac.ed.lurg.landuse;
 import java.io.Serializable;
 import java.util.Objects;
 
+import ac.ed.lurg.types.LandCoverType;
 import ac.ed.lurg.types.LandProtectionType;
 
 public class LandKey implements Serializable {
 	private static final long serialVersionUID = -1173322425708587699L;
-	
+
+	private LandCoverType lcType;
 	private LandProtectionType lpType;
 	private int age; // years since conversion
 	
-	public LandKey(LandProtectionType lpType, int age) {
+	public LandKey(LandCoverType lcType, LandProtectionType lpType, int age) {
 		super();
+		this.lcType = lcType;
 		this.lpType = lpType;
 		this.age = age;
 	}
-	
+
+	public LandCoverType getLcType() {
+		return lcType;
+	}
 	public LandProtectionType getLpType() {
 		return lpType;
 	}
@@ -26,20 +32,15 @@ public class LandKey implements Serializable {
 	}
 
 	@Override
-	public int hashCode() {
-		return Objects.hash(age, lpType);
+	public boolean equals(Object o) {
+		if (this == o) return true;
+		if (o == null || getClass() != o.getClass()) return false;
+		LandKey landKey = (LandKey) o;
+		return age == landKey.age && lcType == landKey.lcType && lpType == landKey.lpType;
 	}
 
 	@Override
-	public boolean equals(Object obj) {
-		if (this == obj)
-			return true;
-		if (obj == null)
-			return false;
-		if (getClass() != obj.getClass())
-			return false;
-		LandKey other = (LandKey) obj;
-		return age == other.age && lpType == other.lpType;
+	public int hashCode() {
+		return Objects.hash(lcType, lpType, age);
 	}
-	
 }
diff --git a/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java b/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java
index 320441b88226bf90cb3834985ba396998e5a4072..3a8519ec3b18d22b115b83e338d8bc1739cd0934 100644
--- a/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java
+++ b/src/ac/ed/lurg/landuse/LandUseBinarySerializer.java
@@ -46,28 +46,25 @@ public class LandUseBinarySerializer {
 				writer.writeDouble(key.getCol());
 				writer.writeDouble(key.getRow());
 
-				Map<LandCoverType, LandCoverTile> tiles = luItem.getLandCoverTiles();
+				LandCoverTile tiles = luItem.getLandCoverTiles();
 
 				// Convertible
 				for (LandCoverType lcType : LandCoverType.values()) {
-					LandCoverTile tile = tiles.get(lcType);
 					for (int a = 0; a < NUM_AGE_CLASSES; a++) {
-						writer.writeDouble(tile.getArea(LandProtectionType.CONVERTIBLE, a));
+						writer.writeDouble(tiles.getArea(lcType, LandProtectionType.CONVERTIBLE, a));
 					}
 				}
 
 				// Protected
 				for (LandCoverType lcType : LandCoverType.values()) {
-					LandCoverTile tile = tiles.get(lcType);
 					for (int a = 0; a < NUM_AGE_CLASSES; a++) {
-						writer.writeDouble(tile.getArea(LandProtectionType.PROTECTED, a));
+						writer.writeDouble(tiles.getArea(lcType, LandProtectionType.PROTECTED, a));
 					}
 				}
 
 				// Unavailable
-				for (LandCoverType lcType : LandCoverType.values()) {
-					LandCoverTile tile = tiles.get(lcType);
-					writer.writeDouble(tile.getTotalArea(LandProtectionType.UNAVAILABLE));
+				for (LandCoverType lcType : LandCoverType.values()) {;
+					writer.writeDouble(tiles.getTotalArea(lcType, LandProtectionType.UNAVAILABLE));
 				}
 			}
 			writer.close();
@@ -116,10 +113,7 @@ public class LandUseBinarySerializer {
 				int col = (int) reader.readDouble();
 				int row = (int) reader.readDouble();
 
-				Map<LandCoverType, LandCoverTile> landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
-				for (LandCoverType lcType: LandCoverType.values()) {
-					landCoverAreas.put(lcType, new LandCoverTile());
-				}
+				LandCoverTile landCoverAreas = new LandCoverTile();
 
 				LandUseItem luItem = luRaster.get(new RasterKey(col, row));
 
@@ -132,7 +126,7 @@ public class LandUseBinarySerializer {
 							for (int age = 0; age < NUM_AGE_CLASSES; age++) {
 								double area = reader.readDouble();
 								if (area > 0)
-									landCoverAreas.get(lcType).setArea(LandProtectionType.CONVERTIBLE, age, area);
+									landCoverAreas.setArea(lcType, LandProtectionType.CONVERTIBLE, age, area);
 							}
 						}
 						break;
@@ -141,7 +135,7 @@ public class LandUseBinarySerializer {
 							for (int age = 0; age < NUM_AGE_CLASSES; age++) {
 								double area = reader.readDouble();
 								if (area > 0)
-									landCoverAreas.get(lcType).setArea(LandProtectionType.PROTECTED, age, area);
+									landCoverAreas.setArea(lcType, LandProtectionType.PROTECTED, age, area);
 							}
 						}	
 						break;
@@ -149,7 +143,7 @@ public class LandUseBinarySerializer {
 						for (LandCoverType lcType : LandCoverType.values()) {
 							double area = reader.readDouble();
 							if (area > 0)
-								landCoverAreas.get(lcType).setArea(LandProtectionType.UNAVAILABLE, 0, area); // TODO Do we need to keep track of age?
+								landCoverAreas.setArea(lcType, LandProtectionType.UNAVAILABLE, 0, area); // TODO Do we need to keep track of age?
 						}
 					}
 				}
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 4ba22ac7fa28549f5669cf535054255e9735d572..99661df0003b43b2acc43322af1bc0fe22abea60 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -5,10 +5,7 @@ import java.util.Collection;
 import java.util.HashMap;
 import java.util.Map;
 import ac.ed.lurg.ModelConfig;
-import ac.ed.lurg.types.CropToDouble;
-import ac.ed.lurg.types.CropType;
-import ac.ed.lurg.types.LandCoverType;
-import ac.ed.lurg.types.LandProtectionType;
+import ac.ed.lurg.types.*;
 import ac.ed.lurg.utils.Interpolator;
 import ac.sac.raster.InterpolatingRasterItem;
 
@@ -17,15 +14,12 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	
 	private Map<CropType, Intensity> intensityMap = new HashMap<CropType, Intensity>();
 	private Map<CropType, Double> cropFractions = new HashMap<CropType, Double>();
-	private Map<LandCoverType, LandCoverTile> landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
+	private LandCoverTile landCoverAreas = new LandCoverTile();
+	private Map<WoodType, Double> forestFractions = new HashMap<>();
 	private Map<LccKey, Double> landCoverChange = new HashMap<LccKey, Double>(); // previous timestep LCC
 	private double protectedFraction; // protected fraction of total area
 	
-	public LandUseItem() {
-		for (LandCoverType lcType : LandCoverType.values()) {
-			landCoverAreas.put(lcType, new LandCoverTile());
-		}
-	}
+	public LandUseItem() {}
 
 	public LandUseItem(LandCoverItem landCover) {
 		this();
@@ -33,7 +27,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 			setCropFraction(CropType.WHEAT, 0.5); // random start as don't have better data
 			setCropFraction(CropType.MAIZE, 0.5);
 			setProtectedFraction(landCover.getProtectedFraction());
-			this.landCoverAreas.putAll(landCover.getLandCoverTiles());
+			this.landCoverAreas = (landCover.getLandCoverTiles());
 		}
 	}
 
@@ -41,16 +35,15 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		this();
 		intensityMap.putAll(luItemToCopy.intensityMap);
 		cropFractions.putAll(luItemToCopy.cropFractions);
-		landCoverAreas.putAll(luItemToCopy.landCoverAreas);
+		landCoverAreas = (luItemToCopy.landCoverAreas);
 		protectedFraction = (luItemToCopy.protectedFraction);
 	}
 
-	public void overwriteLandCoverAreas(Map<LandCoverType, LandCoverTile> landCoverAreas) {
-		this.landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
-		this.landCoverAreas.putAll(landCoverAreas);
+	public void overwriteLandCoverAreas(LandCoverTile landCoverAreas) {
+		this.landCoverAreas = (landCoverAreas);
 	}
 	
-	public Map<LandCoverType, LandCoverTile> getLandCoverTiles() {
+	public LandCoverTile getLandCoverTiles() {
 		return landCoverAreas;
 	}
 	
@@ -71,19 +64,19 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	}
 	
 	public double getLandCoverArea(LandCoverType lcType, LandProtectionType lpType) {
-		return landCoverAreas.get(lcType).getTotalArea(lpType);
+		return landCoverAreas.getTotalArea(lcType, lpType);
 	}
 	
 	public double getLandCoverArea(LandCoverType lcType, LandProtectionType lpType, int age) {
-		return landCoverAreas.get(lcType).getArea(lpType, age);
+		return landCoverAreas.getArea(lcType, lpType, age);
 	}
 	
 	public double getLandCoverArea(LandCoverType lcType) {
-		Double d = landCoverAreas.get(lcType).getTotalArea();
+		Double d = landCoverAreas.getTotalArea(lcType);
 		return d == null ? 0.0 : d;
 	}
 	
-	public double getTotalLandCoverArea() {
+	public double getTotalLandArea() {
 		double d = 0;
 		for (LandCoverType l : LandCoverType.values()) 
 			d += getLandCoverArea(l);
@@ -91,7 +84,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		return d;
 	}
 	
-	public double getTotalLandCoverArea(LandProtectionType lpType) {
+	public double getTotalLandArea(LandProtectionType lpType) {
 		double d = 0;
 		for (LandCoverType l : LandCoverType.values()) 
 			d += getLandCoverArea(l, lpType);
@@ -100,7 +93,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	}
 	
 	public double getLandCoverFract(LandCoverType c) {
-		double totalArea = getTotalLandCoverArea();
+		double totalArea = getTotalLandArea();
 		return totalArea==0 ? 0.0 : getLandCoverArea(c) / totalArea;
 	}
 
@@ -127,14 +120,12 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
         }
            
 		// Note if toType == fromType, area will be removed from fromType and added back in with age 0
-		LandCoverTile toTiles = landCoverAreas.get(toType);
-		LandCoverTile fromTiles = landCoverAreas.get(fromType);
 
-		fromTiles.removeArea(LandProtectionType.CONVERTIBLE, areaMoved);	
+		landCoverAreas.removeArea(fromType, LandProtectionType.CONVERTIBLE, areaMoved);
 		if (ModelConfig.IS_CALIBRATION_RUN) {
-			toTiles.addArea(LandProtectionType.CONVERTIBLE, areaMoved); // preserve age distribution
+			landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, areaMoved); // preserve age distribution
 		} else {
-			toTiles.addArea(LandProtectionType.CONVERTIBLE, 0, areaMoved); // new area starts at age=0
+			landCoverAreas.addArea(toType, LandProtectionType.CONVERTIBLE, 0, areaMoved); // new area starts at age=0
 		}
 		
 		LccKey lccKey = new LccKey(fromType, toType);
@@ -142,26 +133,10 @@ 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 doForestRotation(int age, double area) {
+		landCoverAreas.removeArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, age, area);
+		landCoverAreas.addArea(LandCoverType.TIMBER_FOREST, LandProtectionType.CONVERTIBLE, 0, area);
 	}
 	
 	public void resetLccAreas() {
@@ -178,25 +153,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 
 		double landCover = (d < 0.0) ? 0.0 : d;
 		
-		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>();
-		for (LandCoverType lcType : LandCoverType.values()) {
-			areaMap.put(lcType, getLandCoverArea(lcType));
-		}
-		return areaMap;
+		landCoverAreas.removeAllArea(c, p);
+		landCoverAreas.addArea(c, p, 0, landCover);
 	}
 	
 	public void updateProtectedFraction(int year) {
@@ -214,25 +172,24 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	
 	private void updateProtectedArea() {
 		
-		double previousProtectedTotal = landCoverAreas.values().stream().mapToDouble(o -> o.getTotalArea(LandProtectionType.PROTECTED)).sum();
-		double newProtectedTotal = protectedFraction * getTotalLandCoverArea();
+		double previousProtectedTotal = 0;
+		for (LandCoverType lc : LandCoverType.values()) {
+			previousProtectedTotal += landCoverAreas.getTotalArea(lc, LandProtectionType.PROTECTED);
+		}
+		double newProtectedTotal = protectedFraction * getTotalLandArea();
 		double factor = newProtectedTotal / previousProtectedTotal;
 		
 		for (LandCoverType lcType : LandCoverType.values()) {
 			double prevArea = getLandCoverArea(lcType, LandProtectionType.PROTECTED);
 			double newArea = prevArea * factor;
 			if (newArea > prevArea) {
-				landCoverAreas.get(lcType).moveArea(LandProtectionType.CONVERTIBLE, LandProtectionType.PROTECTED, newArea);
+				landCoverAreas.moveArea(lcType, LandProtectionType.CONVERTIBLE, LandProtectionType.PROTECTED, newArea);
 			} else {
-				landCoverAreas.get(lcType).moveArea(LandProtectionType.PROTECTED, LandProtectionType.CONVERTIBLE, newArea);
+				landCoverAreas.moveArea(lcType, LandProtectionType.PROTECTED, LandProtectionType.CONVERTIBLE, newArea);
 			}			
 		}		
 	}
 	
-	public Map<CropType, Intensity> getIntensityMap() {
-		return intensityMap;
-	}
-	
 	public Intensity getIntensity(CropType crop) {
 		return intensityMap.get(crop);
 	}
@@ -253,18 +210,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		double area = getCropArea(c);
 		return rate * area;
 	}
-	
-	public double getFertiliserAverageRate(CropType... crops) {
-		double fertTotal = 0;
-		double areaTotal = 0;
-		
-		for (CropType c : crops) {
-			fertTotal += getFertiliserAmount(c);
-			areaTotal += getCropArea(c);
-		}
-	
-		return areaTotal > 0 ? fertTotal / areaTotal : 0;
-	}
 
 	public static double getFertiliserTotal(Collection<? extends LandUseItem> items, CropType... crops) {			
 		double total = 0;
@@ -397,23 +342,7 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	}
 
 	public void incrementTilesAge() {
-		for (LandCoverTile tile : landCoverAreas.values()) {
-			tile.cleanUp();
-			tile.incrementAge();
-		}
-	}	
-	
-	public void doForestRotation(int optimalRotation) {
-		LandCoverTile forestTile = landCoverAreas.get(LandCoverType.TIMBER_FOREST);
-		for (int age : forestTile.getAgeKeys()) {
-			if (age >= optimalRotation) {
-				LandKey key = new LandKey(LandProtectionType.CONVERTIBLE, age);
-				double area = forestTile.getArea(key);
-				forestTile.setArea(key, 0.0); // remove area from rotated age class
-				forestTile.addArea(LandProtectionType.CONVERTIBLE, 0, area); // add area to age=0
-			} 
-		}
-		forestTile.cleanUp(); // get rid of 0 area parcels
+		landCoverAreas.incrementAge();
 	}
 
 	private boolean isZeroOrNull(Double d) {
@@ -423,10 +352,10 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 	@Override
 	public void interpolateAll(LandUseItem fromItem, LandUseItem toItem, double factor) {		
 		cropFractions = new HashMap<CropType, Double>();
-		landCoverAreas = new HashMap<LandCoverType, LandCoverTile>();
+		landCoverAreas = new LandCoverTile();
 
-		Double fromCropCover = fromItem.landCoverAreas.get(LandCoverType.CROPLAND).getTotalArea();
-		Double toCropCover = toItem.landCoverAreas.get(LandCoverType.CROPLAND).getTotalArea();
+		Double fromCropCover = fromItem.landCoverAreas.getTotalArea(LandCoverType.CROPLAND);
+		Double toCropCover = toItem.landCoverAreas.getTotalArea(LandCoverType.CROPLAND);
 
 		if (!isZeroOrNull(fromCropCover) && isZeroOrNull(toCropCover)) { // if start with crop but end with none, take starting crop fractions
 			cropFractions.putAll(fromItem.cropFractions);
@@ -456,23 +385,16 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 			}
 		}
 		
-		for (LandCoverType landCover : LandCoverType.values()) {
-			for (int age : landCoverAreas.get(landCover).getAgeKeys()) {
-				Double from = fromItem.landCoverAreas.get(landCover).getArea(LandProtectionType.CONVERTIBLE, age);
-				Double to = toItem.landCoverAreas.get(landCover).getArea(LandProtectionType.CONVERTIBLE, age);
-				Double d = Interpolator.interpolate(from, to, factor);
-				LandCoverTile tile = landCoverAreas.computeIfAbsent(landCover, k -> new LandCoverTile());
-				tile.setArea(LandProtectionType.CONVERTIBLE, age, d);
+		for (LandCoverType lcType : LandCoverType.values()) {
+			for (LandProtectionType lpType : LandProtectionType.values()) {
+				for (int age : landCoverAreas.getAgeKeys()) {
+					Double from = fromItem.landCoverAreas.getArea(lcType, lpType, age);
+					Double to = toItem.landCoverAreas.getArea(lcType, lpType, age);
+					Double d = Interpolator.interpolate(from, to, factor);
+					landCoverAreas.setArea(lcType, lpType, age, d);
+				}
 			}
 		}
-		
-		for (LandCoverType landCover : LandCoverType.values()) {
-			Double from = fromItem.landCoverAreas.get(landCover).getTotalArea(LandProtectionType.PROTECTED);
-			Double to = toItem.landCoverAreas.get(landCover).getTotalArea(LandProtectionType.PROTECTED);
-			Double d = Interpolator.interpolate(from, to, factor);
-			LandCoverTile tile = landCoverAreas.computeIfAbsent(landCover, k -> new LandCoverTile());
-			tile.addArea(LandProtectionType.PROTECTED, 0, d);
-		}
 	}
 	
 	public static double getTotalLandCover(Collection<? extends LandUseItem> landUses, LandCoverType landCover) {
@@ -486,18 +408,6 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		}
 		return total;
 	}
-	
-	public static double getTotalLandArea(Collection<? extends LandUseItem> landUses) {
-		double total = 0;
-		for (LandUseItem a : landUses) {
-			if (a!=null) {
-				Double d = a.getTotalLandCoverArea();
-				if (d!=null)
-					total += d;
-			}
-		}
-		return total;
-	}
 
 	public static double getTotalCropArea(Collection<LandUseItem> landUses, CropType crop) {
 		double total = 0;
@@ -512,8 +422,8 @@ public class LandUseItem implements InterpolatingRasterItem<LandUseItem>, Serial
 		Map<LandCoverType, Double> convertibleAreas = new HashMap<LandCoverType, Double>();
 		Map<LandCoverType, Double> protectedAreas = new HashMap<LandCoverType, Double>();
 		for (LandCoverType lcType : LandCoverType.values()) {
-			convertibleAreas.put(lcType, landCoverAreas.get(lcType).getTotalArea(LandProtectionType.CONVERTIBLE));
-			protectedAreas.put(lcType, landCoverAreas.get(lcType).getTotalArea(LandProtectionType.PROTECTED));
+			convertibleAreas.put(lcType, landCoverAreas.getTotalArea(lcType, LandProtectionType.CONVERTIBLE));
+			protectedAreas.put(lcType, landCoverAreas.getTotalArea(lcType, LandProtectionType.PROTECTED));
 		}
 		return "LandUseItem: [landCoverAreas=" + convertibleAreas + ", protectedArea=" + protectedAreas + "]";
 	}
diff --git a/src/ac/ed/lurg/landuse/LandUseReader.java b/src/ac/ed/lurg/landuse/LandUseReader.java
index 48f278d16a4cdf7a994e59ca22f3c7492c587224..ce6d50f0039d944ccf9e3d432562d222f7aea8e3 100644
--- a/src/ac/ed/lurg/landuse/LandUseReader.java
+++ b/src/ac/ed/lurg/landuse/LandUseReader.java
@@ -25,7 +25,7 @@ public class LandUseReader extends AbstractTabularRasterReader<LandUseItem> {
 			luData.setLandCoverArea(cover, LandProtectionType.CONVERTIBLE, getValueForCol(rowValues, cover.getName()));
 		}
 		
-		luData.setProtectedFraction(getValueForCol(rowValues, "protected") / luData.getTotalLandCoverArea());
+		luData.setProtectedFraction(getValueForCol(rowValues, "protected") / luData.getTotalLandArea());
 
 		for (CropType crop : CropType.getAllItems()) {
   			double cropFract = getValueForCol(rowValues, crop.getGamsName() + "_A");
diff --git a/src/ac/ed/lurg/landuse/WoodUsageData.java b/src/ac/ed/lurg/landuse/WoodUsageData.java
index cddab4878b87b12ff86f674346c2ee4bd8599c3f..8dedd1d95c89c596bc486c93cbf0cb0987e586c2 100644
--- a/src/ac/ed/lurg/landuse/WoodUsageData.java
+++ b/src/ac/ed/lurg/landuse/WoodUsageData.java
@@ -8,12 +8,11 @@ public class WoodUsageData implements Serializable {
 	
 	private double rotaHarvest;
 	private double netImport;
-	private double lucHarvest;
+	private double lucHarvest = 0;
 	
-	public WoodUsageData(double harvest, double netImport, double lucHarvest) {
+	public WoodUsageData(double harvest, double netImport) {
 		this.rotaHarvest = harvest;
 		this.netImport = netImport;
-		this.lucHarvest = lucHarvest;
 	}
 
 	public double getHarvest() {
@@ -26,5 +25,9 @@ public class WoodUsageData implements Serializable {
 
 	public double getLucHarvest() {
 		return lucHarvest;
-	}	
+	}
+
+	public void setLucHarvest(double harvest) {
+		lucHarvest = harvest;
+	}
 }
diff --git a/src/ac/ed/lurg/landuse/WoodUsageReader.java b/src/ac/ed/lurg/landuse/WoodUsageReader.java
index efc978947b35d46552218ace27e2815aeb808206..7ad472806ed518a40393f2c88b9fe373689682b8 100644
--- a/src/ac/ed/lurg/landuse/WoodUsageReader.java
+++ b/src/ac/ed/lurg/landuse/WoodUsageReader.java
@@ -4,9 +4,9 @@ import java.io.BufferedReader;
 import java.io.FileReader;
 import java.io.IOException;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
 
-import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.CountryManager;
 import ac.ed.lurg.country.SingleCountry;
 import ac.ed.lurg.ModelConfig;
@@ -14,20 +14,12 @@ import ac.ed.lurg.country.CompositeCountry;
 import ac.ed.lurg.types.WoodType;
 import ac.ed.lurg.utils.LazyHashMap;
 import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.utils.StringTabularReader;
 
 public class WoodUsageReader {
-	private static final int COUNTRY_COL = 0;
-	private static final int YEAR_COL = 1;
-	private static final int WOOD_TYPE_COL = 2;
-	private static final int NET_IMPORT_COL = 4;
-	private static final int HARVEST_COL = 5;
-	
-	private CompositeCountryManager compositeCountryManager;
-	
-	public WoodUsageReader(CompositeCountryManager compositeCountryManager) {
-		this.compositeCountryManager = compositeCountryManager;
-	}
 	
+	public WoodUsageReader() {}
+
 	public Map<CompositeCountry, Map<WoodType, WoodUsageData>> getWoodUsageData() {
 		
 		LazyHashMap<CompositeCountry, Map<WoodType, WoodUsageData>> usageMap = new LazyHashMap<CompositeCountry, Map<WoodType, WoodUsageData>>() {
@@ -40,64 +32,31 @@ public class WoodUsageReader {
 		};
 		
 		String filename = ModelConfig.WOOD_DEMAND_FILE;
-		try {
-			BufferedReader reader = new BufferedReader(new FileReader(filename));
-			String line, countryName, woodTypeName;
-			Integer year;
-			double harvest, netImport;
-			reader.readLine(); // read header
-			
-			while ((line = reader.readLine()) != null) {
-				String[] tokens = line.split(",");
-				
-				if (tokens.length < 6)
-					LogWriter.printlnError("Too few columns in file: " + filename + ", line: " + line);
-				
-				year = Integer.valueOf(tokens[YEAR_COL]);
-				
-				if (year != ModelConfig.BASE_YEAR)
-					continue;
-				
-				countryName = tokens[COUNTRY_COL];
-				woodTypeName = tokens[WOOD_TYPE_COL];
-				harvest = ModelConfig.IS_FORESTRY_ON ? Double.parseDouble(tokens[HARVEST_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR : 0; // m3 to MtC-eq
-				netImport = ModelConfig.IS_FORESTRY_ON ? Double.parseDouble(tokens[NET_IMPORT_COL]) * ModelConfig.WOOD_BIOMASS_CONVERSION_FACTOR : 0; // m3 to MtC-eq
-				
-				SingleCountry country = CountryManager.getForName(countryName);
-				WoodType woodType = WoodType.getForName(woodTypeName);
-				
-				if (country == null) {
-					// LogWriter.printlnWarning("WoodUsageReader can't find single country: " + countryName);
-					continue;
-				}
-				
-				CompositeCountry cc = compositeCountryManager.getForSingleCountry(country);
-				
-				if (cc == null) {
-					// LogWriter.printlnWarning("WoodUsageReader can't find composite country for: " + country.getCountryName());
-					continue;
-				}
-				
-				Map<WoodType, WoodUsageData> countryData = usageMap.lazyGet(cc);
-				WoodUsageData oldData = countryData.get(woodType);
-				
-				// aggregate if multiple countries for cc
-				if (oldData != null) {
-					netImport += oldData.getNetImport();
-					harvest += oldData.getHarvest();
-				}
-				
-				WoodUsageData newData = new WoodUsageData(harvest, netImport, 0);
-				countryData.put(woodType, newData);
+		StringTabularReader reader = new StringTabularReader(
+				",", new String[] {"Iso3", "Type", "Production", "Demand", "NetImports"});
+		reader.read(filename);
+		List<Map<String, String>> rowList = reader.getRowList();
+		for (Map<String, String> row : rowList) {
+			String countryCode = row.get("Iso3");
+			String woodTypeName = row.get("Type");
+			double production = Double.parseDouble(row.get("Production"));
+			double netImport = Double.parseDouble(row.get("NetImports"));
+			SingleCountry country = CountryManager.getInstance().getForCode(countryCode);
+			WoodType woodType = WoodType.getForName(woodTypeName);
+			CompositeCountry cc = CountryManager.getInstance().getForSingleCountry(country);
+
+			Map<WoodType, WoodUsageData> countryData = usageMap.lazyGet(cc);
+			WoodUsageData oldData = countryData.get(woodType);
+
+			// aggregate if multiple countries for cc
+			if (oldData != null) {
+				netImport += oldData.getNetImport();
+				production += oldData.getHarvest();
 			}
-			reader.close();
-			
-		} catch (IOException e) {
-			LogWriter.printlnError("Failed in reading wood usage data");
-			LogWriter.print(e);
+
+			WoodUsageData newData = new WoodUsageData(production, netImport);
+			countryData.put(woodType, newData);
 		}
-		LogWriter.println("Processed " + filename + ", create " + usageMap.size() + " country wood usage map values");
-		
 		return usageMap;
 	}
 }
diff --git a/src/ac/ed/lurg/output/LandUseOutputer.java b/src/ac/ed/lurg/output/LandUseOutputer.java
index ab42c6fc325877d4c40071da70446b15a30a04bc..4af83db6b5e54ecf9d0d331da4acf5bad6211203 100644
--- a/src/ac/ed/lurg/output/LandUseOutputer.java
+++ b/src/ac/ed/lurg/output/LandUseOutputer.java
@@ -57,10 +57,10 @@ public class LandUseOutputer extends AbstractLandUseOutputer {
 			
 				StringBuffer sbData = new StringBuffer(String.format("%.2f %.2f", lon, lat));
 
-				sbData.append(String.format(" %.8f", item.getTotalLandCoverArea()));
+				sbData.append(String.format(" %.8f", item.getTotalLandArea()));
 				sbData.append(String.format(" %.8f", item.getSuitableArea()));
-				sbData.append(String.format(" %.8f", item.getTotalLandCoverArea(LandProtectionType.PROTECTED)));
-				sbData.append(String.format(" %.8f", item.getTotalLandCoverArea(LandProtectionType.PROTECTED)/item.getTotalLandCoverArea()));
+				sbData.append(String.format(" %.8f", item.getTotalLandArea(LandProtectionType.PROTECTED)));
+				sbData.append(String.format(" %.8f", item.getTotalLandArea(LandProtectionType.PROTECTED)/item.getTotalLandArea()));
 				
 				for (LandCoverType cover : LandCoverType.values()) {
 					sbData.append(String.format(" %.8f", item.getLandCoverArea(cover)));
diff --git a/src/ac/ed/lurg/output/LpjgOutputer.java b/src/ac/ed/lurg/output/LpjgOutputer.java
index 725f2976b10b465e23a8c308aff64b3173ac635f..1151414dce616c4a0885fc8b4951737d2e668c36 100644
--- a/src/ac/ed/lurg/output/LpjgOutputer.java
+++ b/src/ac/ed/lurg/output/LpjgOutputer.java
@@ -135,7 +135,7 @@ public class LpjgOutputer extends AbstractLandUseOutputer {
 				double lon = landUseRaster.getYCoordin(key);
 
 				double expectedArea = landUseRaster.getAreaMha(key);
-				double area = item.getTotalLandCoverArea();
+				double area = item.getTotalLandArea();
 				
 				if (area > 0 && Math.abs((expectedArea-area)/expectedArea) > 0.01) {  // zero area, due to data not being found, already reported so ignored here
 					LogWriter.printlnError("Land cover areas look strange " + key + ": expected=" + expectedArea + ", actual=" + area);
@@ -183,7 +183,7 @@ public class LpjgOutputer extends AbstractLandUseOutputer {
 				double lon = landUseRaster.getYCoordin(key);
 
 				double expectedArea = landUseRaster.getAreaMha(key);
-				double area = item.getTotalLandCoverArea();
+				double area = item.getTotalLandArea();
 				
 				if (area > 0 && Math.abs((expectedArea-area)/expectedArea) > 0.01) {  // zero area, due to data not being found, already reported so ignored here
 					LogWriter.printlnError("Land cover areas look strange " + key + ": expected=" + expectedArea + ", actual=" + area);
diff --git a/src/ac/ed/lurg/types/WoodType.java b/src/ac/ed/lurg/types/WoodType.java
index e32c8dee2db391488ecf26208588cffe77f90899..4e4e7c68682b00e433cc2991523eeefd44c6888d 100644
--- a/src/ac/ed/lurg/types/WoodType.java
+++ b/src/ac/ed/lurg/types/WoodType.java
@@ -5,15 +5,19 @@ import java.util.Map;
 import ac.ed.lurg.ModelConfig;
 
 public enum WoodType {
-	IND_ROUNDWOOD("roundwood", ModelConfig.IND_ROUNDWOOD_DEMAND_ELASTICITY),
-	FUELWOOD("fuelwood", ModelConfig.FUELWOOD_DEMAND_ELASTICITY);
+	IND_ROUNDWOOD("roundwood", ModelConfig.IND_ROUNDWOOD_DEMAND_ELASTICITY, 40, ModelConfig.INIT_ROUNDWOOD_PRICE),
+	FUELWOOD("fuelwood", ModelConfig.FUELWOOD_DEMAND_ELASTICITY, 10, ModelConfig.INIT_FUELWOOD_PRICE);
 	
 	private String name;
 	private double incomeDemandElasticity;
+	private int minRotation;
+	private double initialPrice;
 	
-	WoodType(String name, double incomeDemandElasticity) {
+	WoodType(String name, double incomeDemandElasticity, int minRotation, double initialPrice) {
 		this.name = name;
 		this.incomeDemandElasticity = incomeDemandElasticity;
+		this.minRotation = minRotation;
+		this.initialPrice = initialPrice;
 	}
 
 	private static final Map<String, WoodType> nameCache = new HashMap<String, WoodType>();
@@ -40,4 +44,12 @@ public enum WoodType {
 	public double getIncomeDemandElasticity() {
 		return incomeDemandElasticity;
 	}
+
+	public int getMinRotation() {
+		return minRotation;
+	}
+
+	public double getInitialPrice() {
+		return initialPrice;
+	}
 }
diff --git a/src/ac/ed/lurg/utils/KDTree.java b/src/ac/ed/lurg/utils/KDTree.java
new file mode 100644
index 0000000000000000000000000000000000000000..cd150b9f0443f6848dd18df3d191fc0bf5bf8a83
--- /dev/null
+++ b/src/ac/ed/lurg/utils/KDTree.java
@@ -0,0 +1,107 @@
+package ac.ed.lurg.utils;
+
+import ac.sac.raster.RasterKey;
+
+import java.util.*;
+
+// Stores RasterKeys in a K-D Tree structure. Useful for nearest neighbour search.
+public class KDTree {
+    Node root;
+
+    public KDTree(Set<RasterKey> keys) {
+        root = kdTree(new ArrayList<>(keys), 0);
+    }
+    private class Node {
+        private RasterKey location;
+        private Node leftChild;
+        private Node rightChild;
+        private int axis; // Keeps track of which plane is being bisected
+        private int splittingValue;
+
+        public Node(RasterKey location, Node leftChild, Node rightChild, int axis, int splittingValue) {
+            this.location = location;
+            this.leftChild = leftChild;
+            this.rightChild = rightChild;
+            this.axis = axis;
+            this.splittingValue = splittingValue;
+        }
+
+        public RasterKey getLocation() {
+            return location;
+        }
+
+        public Node getLeftChild() {
+            return leftChild;
+        }
+
+        public Node getRightChild() {
+            return rightChild;
+        }
+
+        public int getAxis() {
+            return axis;
+        }
+
+        public int getSplittingValue() {
+            return splittingValue;
+        }
+    }
+
+    private abstract class RasterKeyComparator implements Comparator<RasterKey> {
+        public abstract int compare(RasterKey key1, RasterKey key2);
+    }
+
+    private Node kdTree(List<RasterKey> points, int depth) {
+        if (points.size() == 0) {
+            return null;
+        }
+
+        int axis = depth % 2; // alternate between x and y axis
+
+        RasterKeyComparator comp;
+        if (axis == 0) {
+            comp = new RasterKeyComparator() {
+                @Override
+                public int compare(RasterKey key1, RasterKey key2) {
+                    return Integer.compare(key1.getCol(), key2.getCol());
+                }
+            };
+        } else {
+            comp = new RasterKeyComparator() {
+                @Override
+                public int compare(RasterKey key1, RasterKey key2) {
+                    return Integer.compare(key1.getRow(), key2.getRow());
+                }
+            };
+        }
+        // Sort by axis
+        List<RasterKey> sortedList = new ArrayList<>(points);
+        Collections.sort(sortedList, comp);
+        int m = sortedList.size();
+        int midPoint = m / 2;
+
+        // Here we recursively bisect the plane, alternating between x and y axes
+        return new Node(sortedList.get(midPoint),
+                kdTree(sortedList.subList(0, midPoint), depth + 1),
+                m <= 2 ? null : kdTree(sortedList.subList(midPoint + 1, m - 1), depth + 1),
+                axis, midPoint);
+    }
+
+    private double getSqDistance(RasterKey p1, RasterKey p2) {
+        return Math.pow(p1.getCol() - p2.getCol(), 2) + Math.pow(p1.getRow() - p2.getRow(), 2);
+    }
+
+//    public RasterKey findNearestNeighbour(RasterKey q, Node n, RasterKey p, double wRef) {
+//
+//        if (n.getLeftChild() == null && n.getRightChild() == null) {
+//            double wPrime = getSqDistance(q, n.getLocation());
+//            if (wPrime < wRef) {
+//                wRef = wPrime;
+//                p = n.getLocation();
+//            }
+//        } else {
+//
+//        }
+//
+//    }
+}
diff --git a/src/ac/ed/lurg/utils/StringTabularReader.java b/src/ac/ed/lurg/utils/StringTabularReader.java
index b5785183f91a34a47651d15131df1dd501e1aef5..cb9489a6b1adfdc56068e4ffd3a1a01da99d1546 100644
--- a/src/ac/ed/lurg/utils/StringTabularReader.java
+++ b/src/ac/ed/lurg/utils/StringTabularReader.java
@@ -95,4 +95,8 @@ public class StringTabularReader {
 		
 		return matchingRows;
 	} 
+	
+	public List<Map<String, String>> getRowList() {
+		return rowList;
+	}
 }
diff --git a/src/ac/ed/lurg/utils/cluster/KMeans.java b/src/ac/ed/lurg/utils/cluster/KMeans.java
index 3eb1f40380c7759945f8fcabd00aa554bf785393..7282aaaa1c46973608f6ec23fda746217babf7d2 100644
--- a/src/ac/ed/lurg/utils/cluster/KMeans.java
+++ b/src/ac/ed/lurg/utils/cluster/KMeans.java
@@ -71,19 +71,19 @@ public class KMeans<K, P extends ClusteringPoint<K>> {
 		int clustersWithPoints = 0;
 		for (Cluster<K,P> c : clusters) {
 			i++;
-			LogWriter.println("[Cluster: " + i+"]");
-			LogWriter.println("[Centroid: " + c.getCentroid() + "]");
-			LogWriter.println("[Points: (" + c.getPoints().size() + " points)");
+			LogWriter.println("[Cluster: " + i+"]", 3);
+			LogWriter.println("[Centroid: " + c.getCentroid() + "]", 3);
+			LogWriter.println("[Points: (" + c.getPoints().size() + " points)", 3);
 			for(P p : c.getPoints())
-				LogWriter.println(p.toString());
+				LogWriter.println(p.toString(), 3);
 			
 			if (c.getPoints().size() > 0)
 				clustersWithPoints++;
 			
-			LogWriter.println("]\n");
+			LogWriter.println("]\n", 3);
 		}
 		
-		LogWriter.println(clusters.size() + " clusters, of which " + clustersWithPoints + " have points");
+		LogWriter.println(clusters.size() + " clusters, of which " + clustersWithPoints + " have points", 1);
 	}
 
     public double calculateClusters(int maxIterations, double tolerance) {
diff --git a/src/ac/ed/lurg/yield/YieldResponsesItem.java b/src/ac/ed/lurg/yield/YieldResponsesItem.java
index e6723643f8a33dee46399e90f8eb44eaeef77d1b..95aece9d848a215fde46f3b1edb12a174300c6e9 100644
--- a/src/ac/ed/lurg/yield/YieldResponsesItem.java
+++ b/src/ac/ed/lurg/yield/YieldResponsesItem.java
@@ -49,4 +49,14 @@ public class YieldResponsesItem implements RasterItem {
 	public double getShockRate(CropType crop) {
 		return getYieldResponseForCrop(crop).getShock();
 	}
+
+	public static YieldResponsesItem getDefault() {
+		YieldResponsesItem item = new YieldResponsesItem();
+		for (CropType crop : CropType.values()) {
+			for (YieldType yieldType : YieldType.values()) {
+				item.setYield(yieldType, crop, 0.0);
+			}
+		}
+		return item;
+	}
 }
\ No newline at end of file
diff --git a/src/ac/sac/raster/RasterSet.java b/src/ac/sac/raster/RasterSet.java
index 4e7eea4edeabdd915e50ffd0d338947300527de0..ea7049fcbd95edfc7072f320c806208f51fb52ca 100755
--- a/src/ac/sac/raster/RasterSet.java
+++ b/src/ac/sac/raster/RasterSet.java
@@ -1,9 +1,6 @@
 package ac.sac.raster;
 
-import java.util.ArrayList;
-import java.util.Collection;
-import java.util.HashMap;
-import java.util.Map;
+import java.util.*;
 
 import ac.ed.lurg.utils.LogWriter;