From bf4b2c820fb72e868c51b7f4c67026a315bbee29 Mon Sep 17 00:00:00 2001
From: s1924442 <b.arendarczyk@sms.ed.ac.uk>
Date: Mon, 29 Mar 2021 21:45:13 +0100
Subject: [PATCH] Added carbon market.

---
 GAMS/IntExtOpt.gms                            |   5 +-
 src/ac/ed/lurg/CarbonMarket.java              |  91 +++++++++++++++
 src/ac/ed/lurg/InternationalMarket.java       | 104 ++++++++++++++++--
 src/ac/ed/lurg/ModelConfig.java               |   2 +
 src/ac/ed/lurg/ModelMain.java                 |  12 +-
 .../ed/lurg/country/AbstractCountryAgent.java |  11 +-
 src/ac/ed/lurg/country/CountryAgent.java      |  27 ++++-
 .../ed/lurg/country/CountryAgentManager.java  |  11 +-
 src/ac/ed/lurg/country/GlobalPrice.java       |   8 +-
 .../lurg/country/gams/GamsCountryInput.java   |  15 ++-
 .../country/gams/GamsLocationOptimiser.java   |   8 +-
 .../lurg/country/gams/GamsLocationOutput.java |   9 +-
 .../country/gams/GamsRasterOptimiser.java     |   2 +-
 .../lurg/country/gams/GamsRasterOutput.java   |   9 +-
 src/ac/ed/lurg/landuse/LandUseItem.java       |  12 ++
 15 files changed, 292 insertions(+), 34 deletions(-)
 create mode 100644 src/ac/ed/lurg/CarbonMarket.java

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 665ea1ad..31fa02be 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -138,7 +138,7 @@ $gdxin
 
 
  POSITIVE VARIABLE cropArea, fertI, irrigI, otherIntensity, ruminantFeed, monogastricFeed, importAmount, exportAmount,
-                   agriLandExpansion, cropIncrease, cropDecrease, pastureDecrease, pastureIncrease, totalFeedDM,
+                   totalFeedDM,
                    landCoverArea, landCoverChange, woodHarvest;
 
 * POSITIVE VARIABLE A;
@@ -305,6 +305,7 @@ $gdxin
  parameter netImportCost(all_types);
  parameter feedCostRate(feed_crop);
  parameter productionShock(all_types);
+ scalar netCarbonFlux;
 
 * Production quantities based on smaller area (before unhandledCropArea adjustment applied)
  totalProd(crop) = sum(location, cropArea.l(crop, location) * yield.l(crop, location));
@@ -325,6 +326,8 @@ $gdxin
  netImportCost(import_crop) = importAmount.l(import_crop) * importPrices(import_crop) - exportAmount.l(import_crop) * exportPrices(import_crop);
  netImportAmount(import_crop) = importAmount.l(import_crop) - exportAmount.l(import_crop);
 
+ netCarbonFlux = SUM(location, carbonFlux.L(location));
+
  Scalar totalCostsLessLU;
 
  totalCostsLessLU = sum(crop,totalProdCost(crop))+ sum(import_crop,netImportCost(import_crop));
diff --git a/src/ac/ed/lurg/CarbonMarket.java b/src/ac/ed/lurg/CarbonMarket.java
new file mode 100644
index 00000000..735a277f
--- /dev/null
+++ b/src/ac/ed/lurg/CarbonMarket.java
@@ -0,0 +1,91 @@
+package ac.ed.lurg;
+
+import java.io.FileInputStream;
+import java.io.FileOutputStream;
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+import java.util.Collection;
+import java.util.Map;
+import java.util.Map.Entry;
+
+import ac.ed.lurg.country.AbstractCountryAgent;
+import ac.ed.lurg.country.CountryAgent;
+import ac.ed.lurg.country.GlobalPrice;
+import ac.ed.lurg.landuse.CropUsageData;
+import ac.ed.lurg.types.CropToDoubleMap;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.LogWriter;
+
+public class CarbonMarket {
+	private GlobalPrice carbonPrice;
+	
+	public CarbonMarket() {
+		if (ModelConfig.IS_CALIBRATION_RUN) {
+			carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
+		} else {
+			carbonPrice = deserializeCarbonPrice();
+		}
+	}
+	
+	public GlobalPrice getCarbonPrice() {
+		return carbonPrice;
+	}
+	
+	/*
+	public void updateCarbonPrice(Collection<CountryAgent> countryAgents, double carbonDemand, Timestep timestep) {
+		double globalNetCarbonFlux;
+		for (CountryAgent ca : countryAgents) {
+			double netCarbonFlux = ca.getNetCarbonFlux();
+			globalNetCarbonFlux += netCarbonFlux;
+		}
+
+			GlobalPrice prevPrice = carbonPrice;
+			double imports = carbonDemand;
+			double exportsBeforeTransportLosses = globalNetCarbonFlux;
+			LogWriter.println(timestep.getYear() + " Updating carbon prices");
+			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exportsBeforeTransportLosses, timestep, globalNetCarbonFlux);
+			LogWriter.println( String.format("Price for carbon updated from %s to %s \n", prevPrice, adjustedPrice));
+			if (adjustedPrice.getStockLevel() < 0)
+				LogWriter.println("Global stocks are below zero carbon, " + timestep.getYear());
+
+			carbonPrice = adjustedPrice;
+	}
+	*/
+	
+	public void serializeGlobalPrices() {
+		try {
+			LogWriter.println("Starting serializing CarbonPrice to " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			ObjectOutputStream out = new ObjectOutputStream(fileOut);
+			out.writeObject(carbonPrice);
+			out.close();
+			fileOut.close();
+			LogWriter.println("Serialized carbon price is saved");
+		} catch (IOException i) {
+			i.printStackTrace();
+		}
+	}
+	
+	@SuppressWarnings("unchecked")
+	private GlobalPrice deserializeCarbonPrice() {
+		try {
+			GlobalPrice initCarbonPrice;
+			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			ObjectInputStream in = new ObjectInputStream(fileIn);
+			initCarbonPrice = (GlobalPrice) in.readObject();
+			in.close();
+			fileIn.close();
+			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			return initCarbonPrice;
+		} catch (IOException i) {
+			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			LogWriter.print(i);
+			return null;
+		} catch (ClassNotFoundException c) {
+			LogWriter.printlnError("GlobalPrice class not found");
+			c.printStackTrace();
+			return null;
+		}
+	}
+}
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 5f3d517f..c069e517 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -6,8 +6,10 @@ import java.io.FileOutputStream;
 import java.io.IOException;
 import java.io.ObjectInputStream;
 import java.io.ObjectOutputStream;
+import java.util.ArrayList;
 import java.util.Collection;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
 import java.util.Map.Entry;
 
@@ -23,22 +25,28 @@ import ac.ed.lurg.utils.LogWriter;
 public class InternationalMarket {
 
 	private Map<CropType, GlobalPrice> worldPrices;
+	private GlobalPrice carbonPrice;
+	private GlobalPrice timberPrice;
 	private PriceShockManager priceShockManager;
 	
+	@SuppressWarnings("unchecked")
 	public InternationalMarket() {
 		if(ModelConfig.IS_CALIBRATION_RUN) {
 			worldPrices = new HashMap<CropType, GlobalPrice>();
 
-
 			Map<CropType, Double> stockLevel = getInitialStockLevels();
 
 			for (CropType crop : CropType.getImportedTypes()) {
 				Double initialStock = stockLevel.get(crop);
 				worldPrices.put(crop, GlobalPrice.createInitial(crop.getInitialPrice(), initialStock));
 			}
+			
+			carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
 		}
 		else {
-			worldPrices = deserializeGlobalPrice();
+			List<Object> deserializedPrices = deserializeGlobalPrice();
+			worldPrices = (Map<CropType, GlobalPrice>) deserializedPrices.get(0);
+			carbonPrice = (GlobalPrice) deserializedPrices.get(1);
 		}
 		priceShockManager = new PriceShockManager();
 	}
@@ -46,6 +54,10 @@ public class InternationalMarket {
 	public  Map<CropType, GlobalPrice> getWorldPrices() {
 		return worldPrices;
 	}
+	
+	public GlobalPrice getCarbonPrice() {
+		return carbonPrice;
+	}
 
 	private Map<CropType, Double> getInitialStockLevels() {
 		Map<CropType, Double> initialStockLevels = new StockReader().read();
@@ -53,7 +65,7 @@ public class InternationalMarket {
 		return initialStockLevels;
 	}
 
-	void determineInternationalTrade(Collection<AbstractCountryAgent> countryAgents, double gen2EcDDemand, Timestep timestep) {
+	void determineInternationalTrade(Collection<AbstractCountryAgent> countryAgents, double gen2EcDDemand, double carbonDemand, Timestep timestep) {
 		CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
 		CropToDoubleMap totalExportCommodities = new CropToDoubleMap();
 		CropToDoubleMap totalProduction = new CropToDoubleMap();
@@ -68,6 +80,7 @@ public class InternationalMarket {
 				double countryNetImports = entry.getValue().getShockedNetImports();
 				totalProduction.incrementValue(c, (entry.getValue().getProductionExpected()-entry.getValue().getProductionShock()));
 
+				//TODO this doesn't work when both imports and exports > 0
 				if (countryNetImports > 0)
 					totalImportCommodities.incrementValue(c, countryNetImports);
 				else
@@ -85,13 +98,34 @@ public class InternationalMarket {
 			double imports = totalImportCommodities.containsKey(crop) ? totalImportCommodities.get(crop) : 0.0;
 			double exportsBeforeTransportLosses = totalExportCommodities.containsKey(crop) ? totalExportCommodities.get(crop) : 0.0;
 			LogWriter.println(timestep.getYear() + " Updating " + crop.getGamsName() + " prices");
-			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exportsBeforeTransportLosses, timestep, totalProduction.get(crop));
+			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exportsBeforeTransportLosses, timestep, totalProduction.get(crop), true);
 			LogWriter.println( String.format("Price for %s updated from %s to %s \n", crop.getGamsName(), prevPrice, adjustedPrice));
 			if (adjustedPrice.getStockLevel() < 0)
 				LogWriter.println("Global stocks are below zero" + crop.getGamsName() + ", " + timestep.getYear());
 
 			worldPrices.put(crop, adjustedPrice);
 		}
+		
+		// Update carbon price
+		double totalCarbonSequestered = 0;
+		double totalCarbonEmitted = 0;
+		for (AbstractCountryAgent ca : countryAgents) {
+			double netCarbonFlux = ca.getNetCarbonFlux();
+			if (netCarbonFlux >= 0)
+				totalCarbonEmitted += netCarbonFlux;
+			else
+				totalCarbonSequestered -= netCarbonFlux;
+		}
+		
+		totalCarbonSequestered = Math.max(totalCarbonSequestered, 0.0000001);
+		GlobalPrice prevPrice = carbonPrice;
+				LogWriter.println(timestep.getYear() + " Updating carbon prices");
+		GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(carbonDemand, totalCarbonSequestered, timestep, totalCarbonSequestered, false);
+		LogWriter.println( String.format("Price for carbon updated from %s to %s \n", prevPrice, adjustedPrice));
+		if (adjustedPrice.getStockLevel() < 0)
+			LogWriter.println("Global stocks are below zero carbon, " + timestep.getYear());
+
+		carbonPrice = adjustedPrice;
 	}
 
 	void writeGlobalMarketFile(Timestep timestep, BufferedWriter outputFile) throws IOException {
@@ -105,14 +139,27 @@ public class InternationalMarket {
 			outputFile.write(sbData.toString());
 			outputFile.newLine();
 		}
+		// Carbon price
+		StringBuffer sbData = new StringBuffer();
+		sbData.append(String.format("%d,%s", timestep.getYear(), "carbon"));
+		sbData.append(String.format(",%.1f,%.1f", carbonPrice.getImportAmount(), carbonPrice.getExportsAfterTransportLosses()));
+		sbData.append(String.format(",%.8f,%.3f", carbonPrice.getExportPrice(), carbonPrice.getStockLevel()));
+
+		outputFile.write(sbData.toString());
+		outputFile.newLine();
 	}
 
 	public void serializeGlobalPrices() {
+		
+		List<Object> pricesToSerialize = new ArrayList<Object>();
+		pricesToSerialize.add(worldPrices);
+		pricesToSerialize.add(carbonPrice);
+		
 		try {
 			LogWriter.println("Starting serializing GlobalPrice to " + ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE);
 			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE);
 			ObjectOutputStream out = new ObjectOutputStream(fileOut);
-			out.writeObject(worldPrices);
+			out.writeObject(pricesToSerialize);
 			out.close();
 			fileOut.close();
 			LogWriter.println("Serialized international market data is saved");
@@ -120,14 +167,29 @@ public class InternationalMarket {
 			i.printStackTrace();
 		}
 	}
-
+	
+/*
+	public void serializeCarbonPrices() {
+		try {
+			LogWriter.println("Starting serializing CarbonPrice to " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			ObjectOutputStream out = new ObjectOutputStream(fileOut);
+			out.writeObject(carbonPrice);
+			out.close();
+			fileOut.close();
+			LogWriter.println("Serialized carbon price is saved");
+		} catch (IOException i) {
+			i.printStackTrace();
+		}
+	}
+*/
 	@SuppressWarnings("unchecked")
-	private Map<CropType, GlobalPrice> deserializeGlobalPrice() {
+	private List<Object> deserializeGlobalPrice() {
 		try {
-			Map<CropType, GlobalPrice> initGlobalPrices;
+			List<Object> initGlobalPrices;
 			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE);
 			ObjectInputStream in = new ObjectInputStream(fileIn);
-			initGlobalPrices = (Map<CropType, GlobalPrice>) in.readObject();
+			initGlobalPrices = (List<Object>) in.readObject();
 			in.close();
 			fileIn.close();
 			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_INTERNATIONAL_MARKET_FILE);
@@ -142,6 +204,30 @@ public class InternationalMarket {
 			return null;
 		}
 	}
+	
+	/*
+	@SuppressWarnings("unchecked")
+	private GlobalPrice deserializeCarbonPrice() {
+		try {
+			GlobalPrice initCarbonPrice;
+			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			ObjectInputStream in = new ObjectInputStream(fileIn);
+			initCarbonPrice = (GlobalPrice) in.readObject();
+			in.close();
+			fileIn.close();
+			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			return initCarbonPrice;
+		} catch (IOException i) {
+			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_CARBON_MARKET_FILE);
+			LogWriter.print(i);
+			return null;
+		} catch (ClassNotFoundException c) {
+			LogWriter.printlnError("GlobalPrice class not found");
+			c.printStackTrace();
+			return null;
+		}
+	}
+	*/
 
 	public void applyPriceShocks(Timestep timestep) {
 		Map<CropType, Double> shocks = priceShockManager.getShocks(timestep);
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 5c172641..77dddcf6 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -7,6 +7,7 @@ import java.util.Enumeration;
 import java.util.HashMap;
 import java.util.Map;
 import java.util.Properties;
+import java.lang.Long;
 
 import ac.ed.lurg.types.ModelFitType;
 import ac.ed.lurg.types.Parameter;
@@ -424,6 +425,7 @@ public class ModelConfig {
 	
 	// Forestry parameters
 	public static final String CONVERSION_COST_FILE = DATA_DIR + File.separator + "conversion_costs.csv";
+	public static final String SERIALIZED_CARBON_MARKET_FILE = CALIB_DIR + File.separator +  "CarbonMarket.ser";
 	public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 0.02); // $1000/tC-eq
 	public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.05); // $1000/tC-eq
 	public static final int FOREST_LOCKIN_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 30); // cannot convert forest after planting for this many years
diff --git a/src/ac/ed/lurg/ModelMain.java b/src/ac/ed/lurg/ModelMain.java
index 8285834d..6d8d436f 100644
--- a/src/ac/ed/lurg/ModelMain.java
+++ b/src/ac/ed/lurg/ModelMain.java
@@ -74,6 +74,7 @@ public class ModelMain {
 	private RasterHeaderDetails desiredProjection;
 
 	private InternationalMarket internationalMarket;
+	//private CarbonMarket carbonMarket;
 	private IrrigationRasterSet currentIrrigationData;
 	private RasterSet<LandUseItem> globalLandUseRaster;
 	private RasterSet<IntegerRasterItem> clusterIdRaster;
@@ -110,6 +111,7 @@ public class ModelMain {
 
 		globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
 		internationalMarket = new InternationalMarket();
+		//carbonMarket = new CarbonMarket();
 		createCountryAgents(compositeCountryManager.getAll());
 	}
 
@@ -141,6 +143,10 @@ public class ModelMain {
 		double previousGen2EcDDemand = (timestep.isInitialTimestep() || ModelConfig.IS_CALIBRATION_RUN ) ? 0: demandManager.getSecondGenBioenergyDemand(timestep.getPreviousTimestep());
 		double gen2EcDDemand = demandManager.getSecondGenBioenergyDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
 		double gen2Increase = (gen2EcDDemand>previousGen2EcDDemand) ? gen2EcDDemand - previousGen2EcDDemand : 0.0;
+		double previousCarbonDemand = (timestep.isInitialTimestep() || ModelConfig.IS_CALIBRATION_RUN ) ? 0: 0.0;
+		double carbonDemand = 1000; //TODO read from file
+		double carbonDemandIncrease = (carbonDemand > previousCarbonDemand) ? carbonDemand - previousCarbonDemand: 0;
+		
 		
 		CarbonFluxRasterSet currentCarbonFluxData = getCarbonFluxData(timestep);
 		
@@ -148,15 +154,15 @@ public class ModelMain {
 		
 		DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = new ConversionCostReader().read();
 		
-		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData, conversionCosts);
-		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, timestep);
+		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase, currentCarbonFluxData, currentWoodYieldData, conversionCosts, carbonDemandIncrease);
+		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep);
 		
 		// loop for iterations.  Could check within a tolerance using internationalMarket.findMaxPriceDiff, not doing so as doesn't find a solution due to inelastic consumption
 		for (int i=0; i < ModelConfig.DEMAND_RECALC_MAX_ITERATIONS; i++) {
 			LogWriter.println("\n++ Re-estimating prices and demand " + i);
 			countryAgents.recalculateDemandForAll(); // recalculate demand from new prices
 			countryAgents.updateNetImportsForAll(); // calculate imports and exports
-			internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, timestep); // calculate prices
+			internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, carbonDemand, timestep); // calculate prices
 		}
 		internationalMarket.applyPriceShocks(timestep);
 		
diff --git a/src/ac/ed/lurg/country/AbstractCountryAgent.java b/src/ac/ed/lurg/country/AbstractCountryAgent.java
index c98997ab..87937387 100644
--- a/src/ac/ed/lurg/country/AbstractCountryAgent.java
+++ b/src/ac/ed/lurg/country/AbstractCountryAgent.java
@@ -18,6 +18,7 @@ public abstract class AbstractCountryAgent {
 	protected Map<CropType, Double> tradeBarriers;
 	protected Map<CommodityType, Double> currentProjectedDemand;
 	protected Map<CropType, CountryPrice> currentCountryPrices;
+	protected CountryPrice currentCarbonPrice;
 	protected Map<CommodityType, Map<CropType, Double>> baseDemandFact;
 	protected Timestep currentTimestep;
 	protected Map<CommodityType, Map<CropType, Double>> currentMinDemandFract;
@@ -39,7 +40,7 @@ public abstract class AbstractCountryAgent {
 		currentMinDemandFract = getMinDemandFraction();
 	}
 
-	private void calculateCountryPrices(Map<CropType, GlobalPrice> worldPrices){
+	private void calculateCountryPrices(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice){
 		Map<CropType, CountryPrice> countryPrices = new HashMap <CropType, CountryPrice>();
 
 		for (CropType c : CropType.getImportedTypes()) {
@@ -51,6 +52,7 @@ public abstract class AbstractCountryAgent {
 		}
 
 		currentCountryPrices = countryPrices;
+		currentCarbonPrice = new CountryPrice(0.0, carbonPrice.getExportPrice());
 	}
 	
 	private void calculateDemand(boolean outputGamsDemand) {
@@ -58,8 +60,8 @@ public abstract class AbstractCountryAgent {
 		currentProjectedDemand = demandManager.getDemand(country, currentTimestep.getYear(), producerPrices, outputGamsDemand);
 	}
 	
-	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, boolean outputGamsDemand) {
-		calculateCountryPrices(worldPrices);
+	protected void calculateCountryPricesAndDemand(Map<CropType, GlobalPrice> worldPrices, GlobalPrice carbonPrice, boolean outputGamsDemand) {
+		calculateCountryPrices(worldPrices, carbonPrice);
 		calculateDemand(outputGamsDemand);
 	}
 	
@@ -169,4 +171,7 @@ public abstract class AbstractCountryAgent {
 	}
 	
 	abstract public Map<CropType, CropUsageData> getCropUsageData();
+	
+	abstract public double getNetCarbonFlux();
+	
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index fb713fec..8478526c 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -114,10 +114,11 @@ public class CountryAgent extends AbstractCountryAgent {
 	
 	public GamsRasterOutput determineProduction(YieldRaster countryYieldSurfaces, RasterSet<IrrigationItem> irrigData, 
 			Map<CropType, GlobalPrice> worldPrices, double globalGen2EcIncrease, RasterSet<CarbonFluxItem> carbonFluxData, 
-			RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts) {
+			RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
+			double carbonDemandIncrease, GlobalPrice carbonPrice) {
 
 		// project demand
-		calculateCountryPricesAndDemand(worldPrices, false);
+		calculateCountryPricesAndDemand(worldPrices, carbonPrice, false);
 		
 		if (saveGamsGdxFiles && ModelConfig.PRICE_ELASTIC_DEMAND)
 			saveGDXFile("demand");
@@ -135,7 +136,7 @@ public class CountryAgent extends AbstractCountryAgent {
 				yieldClusters = calcYieldClusters(irrigData, countryYieldSurfaces);  // this should only be on the first timestep
 			
 			// optimize areas and intensity 
-			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, woodYieldData, conversionCosts);
+			GamsRasterInput input = getGamsRasterInput(irrigData, countryYieldSurfaces, globalGen2EcIncrease, carbonFluxData, woodYieldData, conversionCosts, carbonDemandIncrease);
 			GamsRasterOptimiser opti = new GamsRasterOptimiser(input, yieldClusters);
 			LogWriter.println("Running " + country.getName() + ", currentTimestep " + currentTimestep);
 
@@ -173,7 +174,8 @@ public class CountryAgent extends AbstractCountryAgent {
 	}
 
 	private GamsRasterInput getGamsRasterInput(RasterSet<IrrigationItem> irrigData, YieldRaster countryYieldSurfaces, double gen2EcIncrease,
-			RasterSet<CarbonFluxItem> carbonFluxData, RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts) {
+			RasterSet<CarbonFluxItem> carbonFluxData, RasterSet<WoodYieldItem> woodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
+			double carbonDemandIncrease) {
 		double allowedImportChange;
 
 		if (currentTimestep.isInitialTimestep() || (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION)) {  // initialisation time-step
@@ -203,7 +205,7 @@ public class CountryAgent extends AbstractCountryAgent {
 			// avoid getting stuck at 0
 			if (CropType.ENERGY_CROPS.equals(crop)) {
 				double croplandArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.CROPLAND);
-				double ecMaxExport = gen2EcIncrease * croplandArea / 1500.0 * 2; // twice the share
+				double ecMaxExport = gen2EcIncrease * croplandArea / 1500.0 * 2; // gen2EcIncrease * croplandArea / totalArea * sharePerCountry
 				if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
 					baseTrade = 0.0;
 					ecMaxExport = 0.0; 
@@ -217,6 +219,15 @@ public class CountryAgent extends AbstractCountryAgent {
 			importConstraints.put(crop, new TradeConstraint(baseTrade - changeDown, baseTrade + changeUp));
 		}
 		
+		// Carbon import/export constraints
+		double baseTrade = getNetCarbonFlux();
+		double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
+		double changeUp = 0.0;
+		if (allowedImportChange > 0.0)
+			changeUp = Math.max(baseTrade * allowedImportChange, carbonDemandIncrease * countryArea / 1500 * 2);
+			
+		TradeConstraint carbonTradeConstraint = new TradeConstraint(0, baseTrade + changeUp);
+		
 		// Aggregate minimum land cover for current timestep
 		DoubleMap<Integer, LandCoverType, Double> minimumLandCoverForTimestep = new DoubleMap<Integer, LandCoverType, Double>(); //<Location, LandCoverType, Area>
 		
@@ -231,7 +242,7 @@ public class CountryAgent extends AbstractCountryAgent {
 		}
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
-				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates);	
+				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates, currentCarbonPrice, carbonTradeConstraint);	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, 
 				carbonFluxData, woodYieldData, minimumLandCoverForTimestep, conversionCosts, countryLevelInputs);
 
@@ -392,4 +403,8 @@ public class CountryAgent extends AbstractCountryAgent {
 			 } 
 		 }
 	}
+	
+	public double getNetCarbonFlux() {
+		return previousGamsRasterOutput.getNetCarbonFlux();
+	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/CountryAgentManager.java b/src/ac/ed/lurg/country/CountryAgentManager.java
index a63d938e..1012723d 100644
--- a/src/ac/ed/lurg/country/CountryAgentManager.java
+++ b/src/ac/ed/lurg/country/CountryAgentManager.java
@@ -9,6 +9,7 @@ import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
 
+import ac.ed.lurg.CarbonMarket;
 import ac.ed.lurg.InternationalMarket;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.Timestep;
@@ -47,7 +48,7 @@ public class CountryAgentManager {
 	private CraftyProdManager craftyManager;
 	
 	public CountryAgentManager (CompositeCountryManager compositeCountryManager, AbstractDemandManager demandManager, 
-			CountryBoundaryRaster countryBoundaryRaster, InternationalMarket internationalMarket, 
+			CountryBoundaryRaster countryBoundaryRaster, InternationalMarket internationalMarket,
 			RasterSet<IntegerRasterItem> clusterIdRaster, RasterSet<LandUseItem> globalLandUseRaster) {
 		this.demandManager = demandManager;
 		this.countryBoundaryRaster = countryBoundaryRaster;
@@ -103,7 +104,8 @@ public class CountryAgentManager {
 	}
 
 	public void determineProductionForAll(Timestep timestep, YieldRaster yieldSurfaces, IrrigationRasterSet currentIrrigationData, double gen2Increase,
-			CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts) {
+			CarbonFluxRasterSet currentCarbonFluxData, WoodYieldRasterSet currentWoodYieldData, DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts,
+			double carbonDemandIncrease) {
 		
 		for (AbstractCountryAgent aca : countryAgents) {		
 			aca.setCurrentTimestep(timestep);
@@ -118,7 +120,8 @@ public class CountryAgentManager {
 			RasterSet<WoodYieldItem> woodYieldData = currentWoodYieldData.createSubsetForKeys(countryKeys);
 			
 			try {
-				ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase, carbonFluxData, woodYieldData, conversionCosts);
+				ca.determineProduction(countryYieldSurfaces, irrigData, internationalMarket.getWorldPrices(), gen2Increase, carbonFluxData, woodYieldData, 
+						conversionCosts, carbonDemandIncrease, internationalMarket.getCarbonPrice());
 				
 				// update global rasters
 				globalLandUseRaster.putAll(ca.getLandUses());
@@ -147,7 +150,7 @@ public class CountryAgentManager {
 	public void recalculateDemandForAll() {
 		for (CountryAgent ca : gamsCountryAgents) {
 			ca.shockGDP();
-			ca.calculateCountryPricesAndDemand(internationalMarket.getWorldPrices(), true);
+			ca.calculateCountryPricesAndDemand(internationalMarket.getWorldPrices(), internationalMarket.getCarbonPrice(), true);
 		}
 	}
 	
diff --git a/src/ac/ed/lurg/country/GlobalPrice.java b/src/ac/ed/lurg/country/GlobalPrice.java
index eb54a680..163eb9da 100644
--- a/src/ac/ed/lurg/country/GlobalPrice.java
+++ b/src/ac/ed/lurg/country/GlobalPrice.java
@@ -66,11 +66,15 @@ public class GlobalPrice implements Serializable {
 		return transportLosses;
 	}
 
-	public GlobalPrice createWithUpdatedMarketPrices(double newImports, double newExportAmountBeforeLoss, Timestep timestep, double production) {
+	public GlobalPrice createWithUpdatedMarketPrices(double newImports, double newExportAmountBeforeLoss, Timestep timestep, double production, boolean hasTransportLosses) {
 
 		if (newImports > 0 || newExportAmountBeforeLoss > 0) {
 			double oldDiff = timestep.getTimestep() != timestepId ? 0.0 : exportAmountBeforeLoss - transportLosses - importAmount; // if recomputing for same year need to back our previous adjustment
-			double transportLossRate = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.TRANSPORT_LOSSES : ModelConfig.getParameter(Parameter.TRANSPORT_LOSSES, timestep.getYear());
+			double transportLossRate;
+			if (hasTransportLosses)
+				transportLossRate = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.TRANSPORT_LOSSES : ModelConfig.getParameter(Parameter.TRANSPORT_LOSSES, timestep.getYear());
+			else
+				transportLossRate = 0;
 			double newTransportLosses = newExportAmountBeforeLoss * transportLossRate;
 			double stockChange = newExportAmountBeforeLoss - newTransportLosses - newImports - oldDiff;
 			
diff --git a/src/ac/ed/lurg/country/gams/GamsCountryInput.java b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
index 56e8ef9e..f0fec4b8 100644
--- a/src/ac/ed/lurg/country/gams/GamsCountryInput.java
+++ b/src/ac/ed/lurg/country/gams/GamsCountryInput.java
@@ -20,10 +20,13 @@ public class GamsCountryInput {
 	private Map<CropType, CropUsageData> previousCropUsageData;
 	private Map<CommodityType, Map<CropType, Double>> minDemandFractions;
 	private Map<CropType, Double> subsidyRates;
+	private CountryPrice carbonPrice;
+	private TradeConstraint carbonTradeConstraint;
 
 	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) {
+			Map<CommodityType, Map<CropType, Double>> minDemandFracts, Map<CropType, Double> subsidyRates,
+			CountryPrice carbonPrice, TradeConstraint carbonTradeConstraint) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
@@ -32,6 +35,8 @@ public class GamsCountryInput {
 		this.previousCropUsageData = previousCropUsageData;
 		this.minDemandFractions = minDemandFracts;
 		this.subsidyRates = subsidyRates;
+		this.carbonPrice = carbonPrice;
+		this.carbonTradeConstraint = carbonTradeConstraint;
 	}
 		
 	public CompositeCountry getCountry() { 
@@ -83,4 +88,12 @@ public class GamsCountryInput {
 	public Map<CropType, Double> getSubsidyRates() {
 		return subsidyRates;
 	}
+	
+	public CountryPrice getCarbonPrice() {
+		return carbonPrice;
+	}
+	
+	public TradeConstraint getCarbonTradeConstraint() {
+		return carbonTradeConstraint;
+	}
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 51ee8d37..047774d8 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -313,6 +313,7 @@ public class GamsLocationOptimiser {
 			setGamsParamValue(previousMonogastricFeedP.addRecord(crop.getGamsName()), monogastricFeed, 3);
 			setGamsParamValue(seedAndWasteRateP.addRecord(crop.getGamsName()), seedAndWasteRate, 3);
 		}
+		
 
 		//below in case running without shocks to use original values in model config 
 		
@@ -335,7 +336,7 @@ public class GamsLocationOptimiser {
 		}
 		addScalar(inDB, "maxLandExpansionRate", maxExpansion, 3);
 		
-		addScalar(inDB, "carbonPrice", ModelConfig.CARBON_PRICE, 3);
+		addScalar(inDB, "carbonPrice", countryInput.getCarbonPrice().getExportPrice(), 3);
 		addScalar(inDB, "woodPrice", ModelConfig.WOOD_PRICE, 3);
 		
 		// Carbon fluxes
@@ -591,8 +592,11 @@ public class GamsLocationOptimiser {
 				reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);
 			}							
 		}
+		
+		// Carbon flux
+		double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue();
 
-		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, minimumLandCoverAdditions, reservedDeforested);
+		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, minimumLandCoverAdditions, reservedDeforested, netCarbonFlux);
 		return results;
 	}
 	
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 05a7a54f..7d579e79 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -19,13 +19,15 @@ public class GamsLocationOutput {
 	TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges;
 	DoubleMap<Integer, LandCoverType, Double> minimumLandCover;
 	DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
+	double netCarbonFlux;
 	
 	public GamsLocationOutput(ModelStat status, 
 			Map<Integer, LandUseItem> landUses, 
 			Map<CropType, CropUsageData> cropUsageData,
 			TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChange,
 			DoubleMap<Integer, LandCoverType, Double> minimumLandCover,
-			DoubleMap<Integer, LandCoverType, Double> reservedDeforested) {
+			DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
+			double netCarbonFlux) {
 		super();
 		this.status = status;
 		this.landUses = landUses;
@@ -33,6 +35,7 @@ public class GamsLocationOutput {
 		this.landCoverChanges = landCoverChange;
 		this.minimumLandCover = minimumLandCover;
 		this.reservedDeforested = reservedDeforested;
+		this.netCarbonFlux = netCarbonFlux;
 	}
 	
 	public ModelStat getStatus() {
@@ -57,4 +60,8 @@ public class GamsLocationOutput {
 	public DoubleMap<Integer, LandCoverType, Double> getReservedDeforested() {
 		return reservedDeforested;
 	}
+	
+	public double getNetCarbonFlux() {
+		return netCarbonFlux;
+	}
 }
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 7c3e196e..8afd9b65 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -53,7 +53,7 @@ public class GamsRasterOptimiser {
 	private GamsRasterOutput convertToRaster(GamsLocationInput gamsInput, GamsLocationOutput gamsOutput) {		
 		RasterSet<LandUseItem> newIntensityRaster = allocAreas(gamsInput.getPreviousLandUse(), gamsOutput, gamsInput.getTimestep().getYear());
 
-		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getMinimumLandCover(), gamsOutput.getReservedDeforested());
+		return new GamsRasterOutput(gamsOutput.getStatus(), newIntensityRaster, gamsOutput.getCommoditiesData(), gamsOutput.getMinimumLandCover(), gamsOutput.getReservedDeforested(), gamsOutput.netCarbonFlux);
 	}
 
 	private RasterSet<LandUseItem> createWithSameLandCovers(RasterSet<LandUseItem> toCopy) {
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
index b141c061..c16a6b80 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOutput.java
@@ -18,6 +18,7 @@ public class GamsRasterOutput {
 	private Map<CropType, CropUsageData> cropUsageData;
 	private DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions;
 	DoubleMap<Integer, LandCoverType, Double> reservedDeforested;
+	double netCarbonFlux;
 
 	public GamsRasterOutput(RasterSet<LandUseItem> landUses, Map<CropType, CropUsageData> cropUsageData) {
 		super();
@@ -26,11 +27,13 @@ public class GamsRasterOutput {
 	}
 
 	public GamsRasterOutput(ModelStat status, RasterSet<LandUseItem> intensityRaster, Map<CropType, CropUsageData> cropUsageData,
-			DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions, DoubleMap<Integer, LandCoverType, Double> reservedDeforested) {
+			DoubleMap<Integer, LandCoverType, Double> minimumLandCoverAdditions, DoubleMap<Integer, LandCoverType, Double> reservedDeforested,
+			double netCarbonFlux) {
 		this(intensityRaster, cropUsageData);
 		this.status = status;
 		this.minimumLandCoverAdditions = minimumLandCoverAdditions;
 		this.reservedDeforested = reservedDeforested;
+		this.netCarbonFlux = netCarbonFlux;
 	}
 	
 	public ModelStat getStatus() {
@@ -52,4 +55,8 @@ public class GamsRasterOutput {
 	public DoubleMap<Integer, LandCoverType, Double> getReservedDeforested() {
 		return reservedDeforested;
 	}
+	
+	public double getNetCarbonFlux() {
+		return netCarbonFlux;
+	}
 }
diff --git a/src/ac/ed/lurg/landuse/LandUseItem.java b/src/ac/ed/lurg/landuse/LandUseItem.java
index 1e6c5cb6..56ca4334 100644
--- a/src/ac/ed/lurg/landuse/LandUseItem.java
+++ b/src/ac/ed/lurg/landuse/LandUseItem.java
@@ -455,6 +455,18 @@ 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;
-- 
GitLab