diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index e73089b8ae08d81d3b73fd8c876029ecd6cdebae..56a3f4b935609c82453d518d057c4b40f0e4f867 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -11,6 +11,7 @@
SET land_cover / cropland, pasture, timberForest, carbonForest, natural, naturalNew /;
SET forested(land_cover) / timberForest, carbonForest, natural /;
+ SET managed_forest(land_cover) / timberForest, carbonForest /;
SET natural_forest(land_cover) / carbonForest, natural /;
SET non_timber_forest(land_cover) / cropland, pasture, carbonForest, natural, naturalNew /;
SET exc_natural(land_cover) / cropland, pasture, timberForest, naturalNew /;
@@ -72,7 +73,6 @@
SCALAR forestRotationPeriod;
SCALAR timberMaxNetImport;
SCALAR timberMinNetImport;
- SCALAR loggingCost;
*$gdxin "C:\Users\Bart\Documents\PhD\GAMS testing area\_gams_java_gdb1.gdx"
$gdxin %gdxincname%
@@ -82,7 +82,7 @@ $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth, yieldShock
$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
$load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, minDemandPerOilcrop, seedAndWasteRate
$load previousLandCoverArea, carbonFluxRate, woodYieldLuc, carbonPrice, woodPrice
-$load conversionCost, timberForestYield, woodYieldRota, infrastructureCost, netPresentValueFactor, forestRotationPeriod, timberMaxNetImport, timberMinNetImport, loggingCost
+$load conversionCost, timberForestYield, woodYieldRota, infrastructureCost, netPresentValueFactor, forestRotationPeriod, timberMaxNetImport, timberMinNetImport
$gdxin
SCALAR delta "use to smooth power function see 7.5 www.gams.com/dd/docs/solversconopt.pdf" / 0.00000000001 /;
@@ -121,11 +121,12 @@ $gdxin
baseCost('pasture') = 0.04;
otherIntCost('pasture') = 0.8 + otherICost;
- SCALAR forestExpansionMaxRate / 0.0175 /;
+ SCALAR naturalDecreaseMaxRate / 1 /;
+ SCALAR naturalWoodHarvestMaxRate / 0.02 /;
woodYieldLuc('carbonForest', 'naturalNew', location) = woodYieldLuc('carbonForest', 'pasture', location);
woodYieldLuc('natural', 'naturalNew', location) = woodYieldLuc('natural', 'pasture', location);
conversionCost(land_cover, 'naturalNew') = conversionCost(land_cover, 'natural');
- conversionCost('natural', 'naturalNew') = conversionCost('carbonForest', 'naturalNew');
+ conversionCost('natural', 'naturalNew') = 0;
woodYieldLuc('carbonForest', land_cover, location) = 0;
@@ -134,6 +135,10 @@ $gdxin
carbonFluxRate('pasture', 'naturalNew', location) = carbonFluxRate('pasture', 'natural', location);
carbonFluxRate('cropland', 'naturalNew', location) = carbonFluxRate('cropland', 'natural', location);
+ SCALAR forestRotationCost;
+ forestRotationCost = 0.1;
+ SCALAR woodHarvestParam / 1.897 /;
+
VARIABLES
cropArea(crop, location) total area for crops - Mha
fertI(crop, location) fertilizer intensity for each crop - factor between 0 and 1
@@ -152,9 +157,10 @@ $gdxin
totalConversionCost(location) land cover conversion cost - $1000 per ha
woodHarvestLuc(location) wood harvested due to land cover change - Mt C-eq
woodHarvestRota(location) wood harvested from timber forest rotation - Mt C-eq
- totalWoodHarvest total wood harvested from all activities - Mt C-eq
+ totalWoodHarvest(location) total wood harvested from all activities - Mt C-eq
woodExported total wood sold - Mt C-eq
carbonFlux(location) total carbon flux - Mt C
+ loggingCost(location) cost of logging forest
* 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;
@@ -190,20 +196,20 @@ $gdxin
LAND_COVER_CHANGE_CALC(land_cover, location) calc land cover change
LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) can't create primary natural area
-* FOREST_EXPANSION_CONSTRAINT(forested, location) constraint on forest expansion
-* FOREST_CONTRACTION_CONSTRAINT(forested, location) constraint on forest contraction
NATURAL_CONTRACTION_CONSTRAINT(location)
CONVERSION_COST(location) cost of land cover conversion
WOOD_HARVEST_LUC_CALC(location) calc wood harvested from LUC
WOOD_HARVEST_ROTA_CALC(location) calc wood harvested from timber forest rotation
- TOTAL_WOOD_HARVEST_CALC total wood harvested
-* WOOD_MIN_EXPORT_CONSTRAINT constraint on minimum wood export
+ TOTAL_WOOD_HARVEST_CALC(location) total wood harvested
+ NATURAL_WOOD_HARVEST_CONSTRAINT(location)
WOOD_MAX_EXPORT_CONSTRAINT constraint on maximum wood export
-* WOOD_MAX_TRADE_CONSTRAINT
+ WOOD_MIN_TRADE_CONSTRAINT constraint on minimum wood export
+ WOOD_MAX_TRADE_CONSTRAINT
WOOD_ROTA_EXPORT_CONSTRAINT must export all wood produced from timber forest rotation
-* WOOD_ROTA_TO_LUC_CONSTRAINT
+* WOOD_ROTA_TO_LUC_CONSTRAINT
+ LOGGING_COST_CALC(location)
CARBON_FLUX_CALC(location) calc carbon flux
COST_EQ total cost objective function;
@@ -263,35 +269,39 @@ $gdxin
NATURAL_CONVERSION_CONSTRAINT(exc_natural, location) .. landCoverChange(exc_natural, 'natural', location) =E= 0;
-* FOREST_EXPANSION_CONSTRAINT(forested, location) .. landCoverArea(forested, location) =L= previousLandCoverArea(forested, location) + suitableLandArea(location) * forestExpansionMaxRate;
-* FOREST_CONTRACTION_CONSTRAINT(forested, location) .. landCoverArea(forested, location) =G= previousLandCoverArea(forested, location) - suitableLandArea(location) * forestExpansionMaxRate;
-
- NATURAL_CONTRACTION_CONSTRAINT(location) .. landCoverArea('natural', location) =G= previousLandCoverArea('natural', location) * (1 - forestExpansionMaxRate);
+ NATURAL_CONTRACTION_CONSTRAINT(location) .. landCoverArea('natural', location) =G= previousLandCoverArea('natural', location) * (1 - naturalDecreaseMaxRate);
CONVERSION_COST(location) .. totalConversionCost(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * conversionCost(land_cover_before, land_cover_after)) +
sum(exc_natural, landCoverChange('natural', exc_natural, location) * infrastructureCost(location));
* Wood from land cover change
- WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYieldLuc(land_cover_before, land_cover_after, location));
+* WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * woodYieldLuc(land_cover_before, land_cover_after, location));
+
+ WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= woodYieldLuc('natural', 'naturalNew', location) * landCoverChange('natural', 'naturalNew', location);
* Wood from timberForest rotation = harvest from new forest + existing forest
WOOD_HARVEST_ROTA_CALC(location) .. woodHarvestRota(location) =E= [ sum(non_timber_forest, landCoverChange(non_timber_forest, 'timberForest', location) * woodYieldRota(non_timber_forest, 'timberForest', location)) +
landCoverChange('timberForest', 'timberForest', location) * timberForestYield(location)
] / forestRotationPeriod;
- TOTAL_WOOD_HARVEST_CALC .. totalWoodHarvest =E= sum(location, woodHarvestLuc(location) + woodHarvestRota(location));
+ TOTAL_WOOD_HARVEST_CALC(location) .. totalWoodHarvest(location) =E= woodHarvestLuc(location) + woodHarvestRota(location);
+
+ NATURAL_WOOD_HARVEST_CONSTRAINT(location) .. woodHarvestLuc(location) =L= previousLandCoverArea('natural', location) * woodYieldLuc('natural', 'naturalNew', location) * naturalWoodHarvestMaxRate;
-* WOOD_MIN_EXPORT_CONSTRAINT .. 0 - woodExported =L= timberMaxNetImport;
+ WOOD_MAX_EXPORT_CONSTRAINT .. woodExported =L= sum(location, totalWoodHarvest(location));
- WOOD_MAX_EXPORT_CONSTRAINT .. woodExported =L= totalWoodHarvest;
+ WOOD_MIN_TRADE_CONSTRAINT .. 0 - woodExported =L= timberMaxNetImport;
-* WOOD_MAX_TRADE_CONSTRAINT .. 0 - woodExported =G= timberMinNetImport;
+ WOOD_MAX_TRADE_CONSTRAINT .. 0 - woodExported =G= timberMinNetImport;
* Must sell timber producted in timberForest
WOOD_ROTA_EXPORT_CONSTRAINT .. woodExported =G= sum(location, woodHarvestRota(location));
* WOOD_ROTA_TO_LUC_CONSTRAINT .. sum(location, woodHarvestLuc(location)) =L= sum(location, woodHarvestRota(location)) * 2;
+ LOGGING_COST_CALC(location) .. loggingCost(location) =E= 0.04 + 0.06 * landCoverChange('natural', 'naturalNew', location) / [previousLandCoverArea('natural', location) + delta] *
+ landCoverChange('natural', 'naturalNew', location) * woodYieldLuc('natural', 'naturalNew', location);
+
CARBON_FLUX_CALC(location) .. carbonFlux(location) =E= sum((land_cover_before, land_cover_after), landCoverChange(land_cover_before, land_cover_after, location) * carbonFluxRate(land_cover_before, land_cover_after, location));
@@ -301,9 +311,9 @@ $gdxin
SUM(import_crop, importAmount(import_crop) * importPrices(import_crop) - exportAmount(import_crop) * exportPrices(import_crop)) +
- SUM(location, totalConversionCost(location)) -
+ SUM(location, totalConversionCost(location)) +
- woodExported * woodPrice +
+ [SUM((managed_forest, location), landCoverArea(managed_forest, location) * forestRotationCost) - woodExported * woodPrice] + SUM(location, loggingCost(location)) +
SUM(location, carbonFlux(location)) * carbonPrice
);
@@ -320,7 +330,6 @@ $gdxin
importAmount.L(all_types) = previousImportAmount(all_types);
exportAmount.L(all_types) = previousExportAmount(all_types);
woodHarvestRota.L(location) = previousLandCoverArea('timberForest', location) * timberForestYield(location) / forestRotationPeriod;
- totalWoodHarvest.L = sum(location, woodHarvestRota.L(location));
totalConversionCost.L(location) = previousLandCoverArea('timberForest', location) * conversionCost('timberForest', 'timberForest');
LAND_USE.OptFile = 1;
diff --git a/data/conversion_costs.csv b/data/conversion_costs.csv
index 5e44b81153f28bdba389d806679d0e72198b80f4..353b1a7b8a3dd375f7b6058f1df5cbf9dff3cd0e 100644
--- a/data/conversion_costs.csv
+++ b/data/conversion_costs.csv
@@ -1,26 +1,21 @@
fromLC,toLC,Cost
-pasture,cropland,0.21
-pasture,pasture,0
+pasture,cropland,0.3
pasture,timberForest,0.2
pasture,carbonForest,0.2
pasture,natural,0.2
-cropland,cropland,0
cropland,pasture,0.4
cropland,timberForest,0.33
cropland,carbonForest,0.33
cropland,natural,0.33
-timberForest,cropland,0.01
-timberForest,pasture,0.07
-timberForest,timberForest,0
-timberForest,carbonForest,0
-timberForest,natural,0
-carbonForest,cropland,0.01
-carbonForest,pasture,0.07
-carbonForest,timberForest,0
-carbonForest,carbonForest,0
-carbonForest,natural,0
-natural,cropland,0.01
-natural,pasture,0.07
-natural,timberForest,0
-natural,carbonForest,0
-natural,natural,0
+timberForest,cropland,0.206
+timberForest,pasture,0.176
+timberForest,carbonForest,0.05
+timberForest,natural,0.07
+carbonForest,cropland,0.206
+carbonForest,pasture,0.176
+carbonForest,timberForest,0.05
+carbonForest,natural,0.07
+natural,cropland,0.106
+natural,pasture,0.076
+natural,timberForest,0.07
+natural,carbonForest,0.07
diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index cbd77efa7541fb3b9ff706e1d78bd81ac6730f90..2ebc7bf393d5b3818b77f50c7dadf2d87163f6cb 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -41,8 +41,8 @@ public class InternationalMarket {
worldPrices.put(crop, GlobalPrice.createInitial(crop.getInitialPrice(), initialStock));
}
- carbonPrice = GlobalPrice.createInitial(ModelConfig.CARBON_PRICE, 0);
- timberPrice = GlobalPrice.createInitial(ModelConfig.WOOD_PRICE, 1000);
+ carbonPrice = GlobalPrice.createInitial(ModelConfig.INIT_CARBON_PRICE, 0);
+ timberPrice = GlobalPrice.createInitial(ModelConfig.INIT_WOOD_PRICE, ModelConfig.INIT_WOOD_STOCK);
}
else {
List<Object> deserializedPrices = deserializeGlobalPrice();
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 9788115077d73f5ca277192ebc42d8f920f2d578..8997c16f5ac2ba108634a2138c22f30657ff855c 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -422,18 +422,21 @@ public class ModelConfig {
public static final int DIET_CHANGE_END_YEAR = getIntProperty("DIET_CHANGE_END_YEAR", 2040);
public static final String TARGET_DIET_FILE = getProperty("TARGET_DIET_FILE", DATA_DIR + File.separator + "TargetDiet.txt");
- // Forestry parameters
+ // Forestry and carbon market parameters
public static final String CONVERSION_COST_FILE = DATA_DIR + File.separator + "conversion_costs.csv"; // cost of converting from one land cover to another, $1000/ha
public static final String SERIALIZED_CARBON_MARKET_FILE = CALIB_DIR + File.separator + "CarbonMarket.ser";
public static final String CARBON_DEMAND_FILENAME = getProperty("CARBON_DEMAND_FILENAME", "carbon_demand.csv");
public static final String CARBON_DEMAND_FILE = getProperty("CARBON_DEMAND_FILE", DATA_DIR + File.separator + CARBON_DEMAND_FILENAME);
- public static final double CARBON_PRICE = getDoubleProperty("CARBON_PRICE", 0.02); // $1000/tC-eq
+ public static final double INIT_CARBON_PRICE = getDoubleProperty("INIT_CARBON_PRICE", 0.02); // $1000/tC-eq
public static final String TIMBER_DEMAND_FILENAME = getProperty("TIMBER_DEMAND_FILENAME", "timber_demand.csv");
public static final String TIMBER_DEMAND_FILE = getProperty("TIMBER_DEMAND_FILE", DATA_DIR + File.separator + TIMBER_DEMAND_FILENAME);
- public static final double WOOD_PRICE = getDoubleProperty("WOOD_PRICE", 0.05); // $1000/tC-eq
- public static final int FOREST_ROTATION_PERIOD = getIntProperty("FOREST_LOCKIN_PERIOD", 20); // cannot convert forest after planting for this many years
+ public static final double INIT_WOOD_PRICE = getDoubleProperty("INIT_WOOD_PRICE", 0.05); // $1000/tC-eq
+ public static final double INIT_WOOD_STOCK = getDoubleProperty("INIT_WOOD_STOCK", 0.0); // $1000/tC-eq
+ public static final int FOREST_ROTATION_PERIOD = getIntProperty("FOREST_ROTATION_PERIOD", 30); // cannot convert forest after planting for this many years
public static final int NATURAL_REGROWTH_TIME = getIntProperty("NATURAL_REGROWTH_TIME", 50); // years
+ public static final double NATURAL_REGROWTH_PARAM_P = getDoubleProperty("NATURAL_REGROWTH_PARAM_P", 1.0); // Chapman-Richards paramater. Controls gradient of growth curve
+ public static final double NATURAL_REGROWTH_PARAM_K = -Math.log(1 - Math.pow(0.95, (1 / NATURAL_REGROWTH_PARAM_P))) / NATURAL_REGROWTH_TIME; // Calc Chapman-Richards paramater. 95% yield reached at NATURAL_REGROWTH_TIME
public static final double DISCOUNT_RATE = getDoubleProperty("DISCOUNT_RATE", 0.05);
- public static final double TREE_PLANTING_COST = getDoubleProperty("TREE_PLANTING_COST", 1.0); // $1000/ha
+ public static final double TREE_PLANTING_COST = getDoubleProperty("TREE_PLANTING_COST", 1.2); // $1000/ha
public static final double TREE_LOGGING_COST = getDoubleProperty("TREE_LOGGING_COST", 0.02); // $1000/tC-eq
}
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 277042c354a645afe95e5f4a3170797293fec20c..6ed73fec233ba2a9f9079a326828838f24580d29 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -233,14 +233,7 @@ public class CountryAgent extends AbstractCountryAgent {
// Timber import/export constraints
TradeConstraint timberTradeConstraint;
{
- double baseTrade;
- if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.isInitialTimestep()) {
- double forestedAndNaturalArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.NATURAL) +
- LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.TIMBER_FOREST);
- baseTrade = -timberDemandIncrease * forestedAndNaturalArea / 4000;
- } else {
- baseTrade = -getTimberProduction(); // net imports negative as all exported
- }
+ double baseTrade = -getTimberProduction(); // net imports negative as all exported
double countryArea = LandUseItem.getTotalLandArea(previousGamsRasterOutput.getLandUses().values());
double changeUp = 0.0;
double changeDown = 0.0;
@@ -249,7 +242,14 @@ public class CountryAgent extends AbstractCountryAgent {
changeUp = baseTrade * allowedImportChange;
}
// TODO add initial timber production
- changeDown = Math.min(baseTrade * allowedImportChange, -(timberDemandIncrease * countryArea / 15000 * 2));
+ if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.isInitialTimestep()) {
+ double forestedAndNaturalArea = LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.NATURAL) +
+ LandUseItem.getTotalLandCover(previousGamsRasterOutput.getLandUses().values(), LandCoverType.TIMBER_FOREST);
+ changeDown = -timberDemandIncrease * forestedAndNaturalArea / 4000;
+ if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {changeDown = -timberDemandIncrease;}
+ } else {
+ changeDown = Math.min(baseTrade * allowedImportChange, -(timberDemandIncrease * countryArea / 15000 * 2));
+ }
changeUp = baseTrade * allowedImportChange;
timberTradeConstraint = new TradeConstraint(baseTrade + changeDown, baseTrade - changeUp);
diff --git a/src/ac/ed/lurg/country/ForestManager.java b/src/ac/ed/lurg/country/ForestManager.java
index 5d8b4bf296252e5e0e4bc39197d4f1135f7e031e..1e582d7c1da5c6e9d48a231098b81f66841e29ff 100644
--- a/src/ac/ed/lurg/country/ForestManager.java
+++ b/src/ac/ed/lurg/country/ForestManager.java
@@ -72,11 +72,6 @@ public class ForestManager {
otherNaturalArea.put(clusterId, otherNaturalAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.NATURAL) * (1 - naturalForestFraction)); // Natural is other natural + unmanaged forest so need to subtract the latter
managedForestArea.put(clusterId, timberForestAreaSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST));
timberYieldTimesArea.put(clusterId, timberYieldSoFar + luItem.getUnprotectedLandCoverArea(LandCoverType.TIMBER_FOREST) * timberYieldThisTime);
-
- /*
- if (clusterId == 3) {
- LogWriter.println(""+timberYieldThisTime+","+timberYieldSoFar+","+luItem.getLandCoverArea(LandCoverType.TIMBER_FOREST)+","+timberForestAreaSoFar);
- }*/
} catch (Exception e) {
LogWriter.print(e);
@@ -120,7 +115,9 @@ public class ForestManager {
// Calculates wood yield factor for natural areas based on age
private double naturalWoodYieldFunction(double age) {
- return Math.min(age, ModelConfig.NATURAL_REGROWTH_TIME) / ModelConfig.NATURAL_REGROWTH_TIME;
+ //return Math.min(age, ModelConfig.NATURAL_REGROWTH_TIME) / ModelConfig.NATURAL_REGROWTH_TIME;
+ // Chapman-Richards function
+ return Math.pow((1 - Math.exp(-ModelConfig.NATURAL_REGROWTH_PARAM_K * age)), ModelConfig.NATURAL_REGROWTH_PARAM_P);
}
public Map<Integer, Double> getNaturalYieldFactorForTimestep(Timestep timestep) {
@@ -143,13 +140,15 @@ public class ForestManager {
public void updateForestHistory(ArrayList<LandCoverChangeItem> gamsLandCoverChanges, Map<Integer, Double> baseTimberYield) {
- // Update timber forest yield
- for (Integer locId : timberForestYield.keySet()) {
- double area = timberForestArea.get(locId);
- double currentYield = timberForestYield.get(locId);
- double baseYield = baseTimberYield.get(locId);
- double updatedYield = updateTimberForestYield(area, currentYield, baseYield);
- timberForestYield.put(locId, updatedYield);
+ if (!ModelConfig.IS_CALIBRATION_RUN) {
+ // Update timber forest yield - the yield effect of previous land use gradually decreases
+ for (Integer locId : timberForestYield.keySet()) {
+ double area = timberForestArea.get(locId);
+ double currentYield = timberForestYield.get(locId);
+ double baseYield = baseTimberYield.get(locId);
+ double updatedYield = updateTimberForestYield(area, currentYield, baseYield);
+ timberForestYield.put(locId, updatedYield);
+ }
}
// Removed timber forest
@@ -161,6 +160,7 @@ public class ForestManager {
timberForestArea.put(locId, currentArea - removedArea);
}
}
+
// New timber forest
for (LandCoverChangeItem item : gamsLandCoverChanges) {
if (item.getToLandCover().equals(LandCoverType.TIMBER_FOREST.getName()) && !item.getFromLandCover().equals(LandCoverType.TIMBER_FOREST.getName())) {
@@ -182,6 +182,10 @@ public class ForestManager {
// Safety check
timberForestArea.forEach((k, v) -> {if (v < 0) LogWriter.printlnError("ForestManager.updateForestHistory(): timberForestArea < 0");});
+ // Don't update natural areas during calibration run
+ if (ModelConfig.IS_CALIBRATION_RUN)
+ return;
+
// Update natural areas
Map<Integer, Double> newNaturalArea = new HashMap<Integer, Double>();
Map<Integer, Double> clearedNaturalArea = new HashMap<Integer, Double>();
@@ -193,10 +197,6 @@ public class ForestManager {
String toLC = item.getToLandCover();
double area = item.getArea();
- if (locId == 13) {
- int foo = 1;
- }
-
// new natural area
if (!fromLC.equals("natural") && toLC.equals("naturalNew")) {
double areaSoFar = (newNaturalArea.containsKey(locId)) ? newNaturalArea.get(locId) : 0;
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 0d44762fd7144c34bcb3c4598dc5a56ea20888bb..b5c5e20309dd64d6e952406be88a41e3e9b93906 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -180,7 +180,7 @@ public class GamsLocationOptimiser {
}
// Infrastructure cost
- double infrastructureCost = landUseItem.getLandCoverFract(LandCoverType.NATURAL) * 0.09 + 0.01;
+ double infrastructureCost = landUseItem.getLandCoverFract(LandCoverType.NATURAL) * 0.1;
setGamsParamValue(infrastructureCostP.addRecord(locString), infrastructureCost, 5);
}
@@ -431,7 +431,7 @@ public class GamsLocationOptimiser {
// Land cover conversion cost
GAMSParameter conversionCostP = inDB.addParameter("conversionCost", 2);
- /*
+
DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = inputData.getConversionCosts();
for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap : conversionCosts.getMap().entrySet()) {
@@ -445,7 +445,7 @@ public class GamsLocationOptimiser {
setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
}
}
- */
+ /*
for (LandCoverType fromLc : LandCoverType.getConvertibleTypes()) {
for (LandCoverType toLc : LandCoverType.getConvertibleTypes()) {
double cost = 0;
@@ -477,7 +477,7 @@ public class GamsLocationOptimiser {
cost += ModelConfig.CROP_INCREASE_COST;
break;
case TIMBER_FOREST:
- cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
+ cost += ModelConfig.TREE_PLANTING_COST;
break;
case CARBON_FOREST:
cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
@@ -490,13 +490,22 @@ public class GamsLocationOptimiser {
if (toLc.equals(fromLc))
cost = 0;
-
+
+ if (toLc.equals(LandCoverType.NATURAL) && fromLc.isManagedForest())
+ cost = 0;
+
// Forest rotation cost
if (toLc.isManagedForest() && fromLc.isManagedForest() & toLc.equals(fromLc))
cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD;
- if (!toLc.isNatural() && fromLc.isNatural())
- cost += ModelConfig.AGRI_EXPANSION_COST_BASE;
+ if (toLc.equals(LandCoverType.NATURAL) && fromLc.equals(LandCoverType.TIMBER_FOREST))
+ cost += ModelConfig.TREE_PLANTING_COST / ModelConfig.FOREST_ROTATION_PERIOD * 0.9;
+
+ //if (toLc.equals(LandCoverType.TIMBER_FOREST))
+ // cost = cost / ModelConfig.FOREST_ROTATION_PERIOD;
+
+
+ cost = 0;
Vector<String> v = new Vector<String>();
v.add(fromLc.getName());
@@ -504,7 +513,7 @@ public class GamsLocationOptimiser {
setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
}
}
-
+ */
// Net present value factor for timber i.e. multiply by value of potential timber harvest to get the future discounted value over the forest rotation period
double npvFactor = 0;