diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index eccb3b92c69a92622cb724ca144eac08af6f7861..4dc891f2e971fea8f0fa771d576abff138bd8680 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -59,7 +59,7 @@
PARAMETER woodYieldMaxParam(location);
PARAMETER woodYieldSlopeParam(location);
PARAMETER woodYieldShapeParam(location);
- PARAMETER previousRotationIntensity(location);
+ PARAMETER previousRotationIntensity(protection, location);
PARAMETER maxLandCoverChange(land_cover_before, land_cover_after, protection, location);
PARAMETER maxLandCoverArea(land_cover, protection, location);
@@ -85,6 +85,8 @@
SCALAR fertLimConventional;
SCALAR irrigLimConventional;
SCALAR otherLimConventional;
+ SCALAR timberLimProtected;
+ SCALAR timberLimUnprotected;
SCALAR forestManagementCost cost $1000 per ha;
SCALAR carbonForestMaxProportion maximum proportion of land cover as carbon forest;
@@ -103,7 +105,7 @@ $load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImpo
$load meatEfficency, otherICost, irrigCost, irrigConstraintLocation, irrigConstraintFarmingType, irrigMaxRate, fertiliserUnitCost, domesticPriceMarkup, seedAndWasteRate
$load previousLandCoverArea, carbonCreditRate, conversionCost, woodYieldMaxParam, woodYieldSlopeParam, woodYieldShapeParam, tradeAdjustmentCostRate
$load forestManagementCost, carbonForestMaxProportion, previousRotationIntensity, maxLandCoverChange, maxLandCoverArea
-$load fertLimRestricted, otherLimRestricted, fertLimConventional, otherLimConventional
+$load fertLimRestricted, otherLimRestricted, fertLimConventional, otherLimConventional, timberLimProtected, timberLimUnprotected
$load solarEnergyDensity, agrivoltaicsYieldFactor, prodCost, animalProdCost, solarCostRate
$load discountRate, energyPrice, monoculturePenalty, agrivoltaicsCropCostFactor, previousProduction, pastureIntensityPenalty
$gdxin
@@ -162,8 +164,8 @@ $gdxin
landCoverChange(land_cover_before, land_cover_after, protection, location) land cover change in Mha
totalConversionCost(location) land cover conversion cost - $1000 per ha
- woodYieldRota(location) wood yield m3 per ha
- rotationIntensity(location) equivalent to 1 divided by rotation period
+ woodYieldRota(protection, location) wood yield m3 per ha
+ rotationIntensity(protection, location) equivalent to 1 divided by rotation period
woodSupply total wood supply - million m3
carbonCredits(location) carbon credits generated - Mt C
@@ -205,7 +207,7 @@ $gdxin
CONVERSION_COST_CALC(location) cost of land cover conversion
MAX_CROPLAND_CONSTRAINT(protection, location) limits on maximum cropland due to slope
- WOOD_YIELD_CALC(location) wood yield at rotation - tC per ha
+ WOOD_YIELD_CALC(protection, location) wood yield at rotation - tC per ha
WOOD_SUPPLY_CALC wood supply calculation
WOOD_DEMAND_CONSTRAINT wood demand constraint
@@ -310,13 +312,15 @@ $gdxin
************* Forestry ***********************************
- WOOD_YIELD_CALC(location) .. woodYieldRota(location) =E= woodYieldMaxParam(location) * rpower(1 - exp(woodYieldSlopeParam(location) / rotationIntensity(location)), woodYieldShapeParam(location));
+ WOOD_YIELD_CALC(protection, location) .. woodYieldRota(protection, location) =E= woodYieldMaxParam(location) * rpower(1 - exp(woodYieldSlopeParam(location) / rotationIntensity(protection, location)), woodYieldShapeParam(location));
* Upper and lower bounds on rotation intensity
- rotationIntensity.UP(location) = 1 / 10;
- rotationIntensity.LO(location) = 1 / 160;
-
- WOOD_SUPPLY_CALC .. woodSupply =E= sum(location, woodYieldRota(location) * rotationIntensity(location) * landCoverArea('timberForest', 'unprotected', location));
+
+ rotationIntensity.LO(protection, location) = 1 / 160;
+ rotationIntensity.UP('unprotected', location) = timberLimUnprotected;
+ rotationIntensity.UP('protected', location) = timberLimProtected;
+
+ WOOD_SUPPLY_CALC .. woodSupply =E= sum((protection, location), woodYieldRota(protection, location) * rotationIntensity(protection, location) * landCoverArea('timberForest', protection, location));
WOOD_DEMAND_CONSTRAINT .. woodSupply =E= demand('wood') + exportAmount('wood') - importAmount('wood');
@@ -366,7 +370,7 @@ $gdxin
* Cost of imports
- sum(import_types, importPrices(import_types) * importAmount(import_types))
* Cost of wood production
- - sum(location, forestManagementCost * rotationIntensity(location) * landCoverArea('timberForest', 'unprotected', location))
+ - sum((protection, location), forestManagementCost * rotationIntensity(protection, location) * landCoverArea('timberForest', protection, location))
* Cost of solar energy
- sum((solar_land, location), solarCostRate(solar_land, location) * landCoverArea(solar_land, 'unprotected', location))
* Land cover conversion costs
@@ -400,9 +404,9 @@ $gdxin
importAmount.L(all_types) = previousImportAmount(all_types);
exportAmount.L(all_types) = previousExportAmount(all_types);
netImportAmount.L(import_types) = previousImportAmount(import_types) - previousExportAmount(import_types);
- rotationIntensity.L(location) = previousRotationIntensity(location);
- rotationIntensity.L(location) $ (previousRotationIntensity(location) = 0) = 1 / 160;
- woodYieldRota.L(location) = woodYieldMaxParam(location) * rpower(1 - exp(woodYieldSlopeParam(location) / rotationIntensity.L(location)), woodYieldShapeParam(location));
+ rotationIntensity.L(protection, location) = previousRotationIntensity(protection, location);
+ rotationIntensity.L(protection, location) $ (previousRotationIntensity(protection, location) = 0) = 1 / 160;
+ woodYieldRota.L(protection, location) = woodYieldMaxParam(location) * rpower(1 - exp(woodYieldSlopeParam(location) / rotationIntensity.L(protection, location)), woodYieldShapeParam(location));
totalProduction.L(agri_types) = previousProduction(agri_types);
option nlp = conopt4;
@@ -447,7 +451,7 @@ $gdxin
+ (1 - importFraction(feed_crop)) * prodCostRate(feed_crop)) + totalProduction.l(animal) * animalProdCost(animal);
If(woodSupply.L > 0,
- woodProdCost = sum(location, forestManagementCost * rotationIntensity.L(location) * landCoverArea.L('timberForest', 'unprotected', location));
+ woodProdCost = sum((protection, location), forestManagementCost * rotationIntensity.L(protection, location) * landCoverArea.L('timberForest', protection, location));
);
netCarbonCredits = SUM(location, carbonCredits.L(location));
@@ -463,7 +467,7 @@ $gdxin
* Diagnostics
parameter totalLandCoverChange(land_cover, land_cover_after);
parameter landCoverDiff(land_cover);
- parameter rotationPeriod(location);
+ parameter rotationPeriod(protection, location);
parameter irrigWaterWithdrawn(crop, location);
parameter totalIrrig;
parameter fertiliserUsed(crop);
@@ -474,7 +478,7 @@ $gdxin
parameter woodProd(location);
totalLandCoverChange(land_cover, land_cover_after) = sum((protection, location), landCoverChange.L(land_cover, land_cover_after, protection, location));
landCoverDiff(land_cover) = sum((protection, location), landCoverArea.L(land_cover, protection, location)) - sum((protection, location), previousLandCoverArea(land_cover, protection, location));
- rotationPeriod(location) = 1 / rotationIntensity.L(location);
+ rotationPeriod(protection, location) = 1 / rotationIntensity.L(protection, location);
irrigWaterWithdrawn(crop, location) = sum(farming_type, irrigMaxRate(crop, location) * irrigI.L(crop, farming_type, location) * cropArea.L(crop, farming_type, location) * 0.01);
totalIrrig = sum((crop, location), irrigWaterWithdrawn(crop, location));
fertiliserUsed(crop) = sum((farming_type, location), fertI.L(crop, farming_type, location) * cropArea.L(crop, farming_type, location));
@@ -482,4 +486,4 @@ $gdxin
abandonedArea(crop) = sum((farming_type, location)$(otherIntensity.L(crop, farming_type, location) = 0), cropArea.L(crop, farming_type, location));
otherIntUsed(crop, location) = sum(farming_type, otherIntensity.L(crop, farming_type, location) * cropArea.L(crop, farming_type, location));
totalOtherI = sum((crop, location), otherIntUsed(crop, location));
- woodProd(location) = woodYieldRota.L(location) * landCoverArea.L('timberForest', 'unprotected', location) * rotationIntensity.L(location);
+ woodProd(location) = sum(protection, woodYieldRota.L(protection, location) * landCoverArea.L('timberForest', protection, location) * rotationIntensity.L(protection, location));
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 7cdd8324aec2396210c20e9f718ac0ed3888d584..fb0edc718ed3f5b16bcde24b24f46704bf1f976e 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -573,6 +573,8 @@ public class ModelConfig {
public static final double WOOD_TRADE_BARRIER = getDoubleProperty("WOOD_TRADE_BARRIER", 0.2); //$1000/m3
public static final double INIT_WOOD_PRICE = IS_FORESTRY_ON ? getDoubleProperty("INIT_WOOD_PRICE", 0.03) : 0.0; // $1000/m3
public static final int CARBON_WOOD_MAX_TIME = getIntProperty("CARBON_WOOD_MAX_TIME", 160); // upper data limit, years
+ public static final double PROTECTED_TIMBER_INTENSITY_LIMIT = getDoubleProperty("PROTECTED_TIMBER_INTENSITY_LIMIT", 0.02); // = 1/rotation_period
+ public static final double UNPROTECTED_TIMBER_INTENSITY_LIMIT = getDoubleProperty("UNPROTECTED_TIMBER_INTENSITY_LIMIT", 0.1);
// Carbon
public static final boolean IS_CARBON_ON = !IS_CALIBRATION_RUN && getBooleanProperty("IS_CARBON_ON", false);
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 3a321d8d60d795fc83e912938fad9acffac1d65f..4ec7fa77b190db4ca30dd611a50aca7f96d57fe2 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -120,7 +120,7 @@ public class GamsLocationOptimiser {
GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 3);
GAMSParameter maxLandCoverAreaP = inDB.addParameter("maxLandCoverArea", 3);
- GAMSParameter prevRotaIntP = inDB.addParameter("previousRotationIntensity", 1);
+ GAMSParameter prevRotaIntP = inDB.addParameter("previousRotationIntensity", 2);
@SuppressWarnings("unused")
double totalAgriLand = 0;
@@ -207,8 +207,15 @@ public class GamsLocationOptimiser {
}
// Forestry
- ForestDetail forestDetail = (ForestDetail) landUseItem.getLandCoverDetail(new LandKey(LandCoverType.TIMBER_FOREST, LandProtectionType.UNPROTECTED));
- setGamsParamValue(prevRotaIntP.addRecord(locString), forestDetail != null ? forestDetail.getRotationIntensity() : 0.02, -1);
+ for (LandKey landKey : landUseItem.getLandKeysForLandCover(LandCoverType.TIMBER_FOREST)) {
+ ForestDetail forestDetail = (ForestDetail) landUseItem.getLandCoverDetail(landKey);
+ Vector<String> v = new Vector<>();
+ v.add(landKey.getLpType().getGamsName());
+ v.add(locString);
+
+ setGamsParamValue(prevRotaIntP.addRecord(v), forestDetail != null ? forestDetail.getRotationIntensity() : 0.02, -1);
+ }
+
}
if (DEBUG) LogWriter.println(String.format(" Total agricultural %.1f,\t suitable %.1f", totalAgriLand, totalSuitable));
@@ -465,6 +472,8 @@ public class GamsLocationOptimiser {
setGamsParamValue(importPriceP.addRecord("wood"), countryInput.getWoodPrice().getImportPrice() , -1);
setGamsParamValue(exportPriceP.addRecord("wood"), countryInput.getWoodPrice().getExportPrice(), -1);
+ addScalar(inDB, "timberLimProtected", ModelConfig.PROTECTED_TIMBER_INTENSITY_LIMIT, -1);
+ addScalar(inDB, "timberLimUnprotected", ModelConfig.UNPROTECTED_TIMBER_INTENSITY_LIMIT, -1);
addScalar(inDB, "forestManagementCost",
ModelConfig.getAdjParam("FOREST_MANAGEMENT_COST") * calibrationItem.getForestryCostFactor(), -1);
@@ -835,12 +844,14 @@ public class GamsLocationOptimiser {
}
// Rotation intensity
- Map<Integer, Double> rotationIntensities = new HashMap<>();
+ Map<LandProtectionType, Map<Integer, Double>> rotationIntensities = new HashMap<>();
GAMSVariable varRota = outDB.getVariable("rotationIntensity");
for (GAMSVariableRecord rec : varRota) {
- int loc = Integer.parseInt(rec.getKey(0));
+ LandProtectionType lpType = LandProtectionType.getForGamsName(rec.getKey(0));
+ int loc = Integer.parseInt(rec.getKey(1));
double intensity = rec.getLevel();
- rotationIntensities.put(loc, intensity);
+ Map<Integer, Double> locMap = rotationIntensities.computeIfAbsent(lpType, k -> new HashMap<>());
+ locMap.put(loc, intensity);
}
// Carbon
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
index 8de888c98784cb0df2deebc00308f1850858d0e3..46c51c92dd205f539fec84d1b4866ad9ac2da587 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOutput.java
@@ -4,6 +4,7 @@ import java.util.List;
import java.util.Map;
import ac.ed.lurg.landuse.*;
+import ac.ed.lurg.types.LandProtectionType;
import com.gams.api.GAMSGlobals.ModelStat;
import ac.ed.lurg.types.CropType;
@@ -17,14 +18,14 @@ public class GamsLocationOutput {
Map<Integer, List<LandCoverChangeItem>> landCoverChanges;
private CarbonUsageData carbonUsageData;
private WoodUsageData woodUsageData;
- private Map<Integer, Double> rotationIntensities;
+ private Map<LandProtectionType, Map<Integer, Double>> rotationIntensities;
public GamsLocationOutput(ModelStat status,
Map<Integer, LandUseItem> landUses,
Map<CropType, CropUsageData> cropUsageData,
Map<Integer, List<LandCoverChangeItem>> landCoverChanges,
CarbonUsageData carbonUsageData, WoodUsageData woodUsageData,
- Map<Integer, Double> rotationIntensities) {
+ Map<LandProtectionType, Map<Integer, Double>> rotationIntensities) {
super();
this.status = status;
this.landUses = landUses;
@@ -58,7 +59,7 @@ public class GamsLocationOutput {
return woodUsageData;
}
- public Map<Integer, Double> getRotationIntensities() {
+ public Map<LandProtectionType, Map<Integer, Double>> getRotationIntensities() {
return rotationIntensities;
}
}
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 4d7495fb6c79c0c55ef65f15af2786076cc44e96..cbfb7b77a111eac708bd6b4d9eece42a57e29311 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -89,13 +89,15 @@ public class GamsRasterOptimiser {
}
// Forest rotation intensity and yields
- Map<Integer, Double> rotaIntensities = gamsOutput.getRotationIntensities();
+ Map<LandProtectionType, Map<Integer, Double>> rotaIntensities = gamsOutput.getRotationIntensities();
for (RasterKey key : newLandUseRaster.keySet()) {
LandUseItem luItem = newLandUseRaster.get(key);
Integer locId = mapping.get(key).getInt();
- double rotaInt = rotaIntensities.getOrDefault(locId, 0.0);
- ForestDetail detail = (ForestDetail) luItem.getLandCoverDetail(new LandKey(LandCoverType.TIMBER_FOREST, LandProtectionType.UNPROTECTED));
- detail.setRotationIntensity(rotaInt);
+ for (LandProtectionType lpType : LandProtectionType.values()) {
+ double rotaInt = rotaIntensities.containsKey(lpType) ? rotaIntensities.get(lpType).getOrDefault(locId, 0.0) : 0.0;
+ ForestDetail detail = (ForestDetail) luItem.getLandCoverDetail(new LandKey(LandCoverType.TIMBER_FOREST, lpType));
+ detail.setRotationIntensity(rotaInt);
+ }
}
// Crop areas and intensity
diff --git a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
index 74a0337a12cafb2f6636b6927e96221612e1f4cc..2a248ff6eddbce7f17a22b2e9af019c808959bd6 100644
--- a/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
+++ b/src/ac/ed/lurg/forestry/ForestryDataOutputer.java
@@ -32,7 +32,7 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
BufferedWriter outputWriter = null;
try {
outputWriter = new BufferedWriter(new FileWriter(outputFileName, false));
- outputWriter.write("Lon Lat RotationIntensity TimberForestArea Yield");
+ outputWriter.write("Lon Lat Protection RotationIntensity TimberForestArea Yield");
outputWriter.newLine();
for (Map.Entry<RasterKey, LandUseItem> entry : landUseRaster.entrySet()) {
@@ -46,18 +46,19 @@ public class ForestryDataOutputer extends AbstractLandUseOutputer {
double lat = landUseRaster.getXCoordin(key);
double lon = landUseRaster.getYCoordin(key);
- LandKey landKey = new LandKey(LandCoverType.TIMBER_FOREST, LandProtectionType.UNPROTECTED);
- ForestDetail forestDetail = (ForestDetail) luItem.getLandCoverDetail(landKey);
- double rotaInt = forestDetail.getRotationIntensity();
- double area = luItem.getLandCoverArea(landKey);
- if (area == 0) { // skip if no forest
- continue;
- }
- double yield = wyItem.getYield(rotaInt);
- outputWriter.write(String.format("%.2f %.2f %.14f %.14f %.14f", lat, lon, rotaInt, area, yield));
- outputWriter.newLine();
-
+ for (LandKey landKey : luItem.getLandKeysForLandCover(LandCoverType.TIMBER_FOREST)) {
+ ForestDetail forestDetail = (ForestDetail) luItem.getLandCoverDetail(landKey);
+ double rotaInt = forestDetail.getRotationIntensity();
+ double area = luItem.getLandCoverArea(landKey);
+ if (area == 0) { // skip if no forest
+ continue;
+ }
+ double yield = wyItem.getYield(rotaInt);
+
+ outputWriter.write(String.format("%.2f %.2f %s %.14f %.14f %.14f", lat, lon, landKey.getLpType().getGamsName(), rotaInt, area, yield));
+ outputWriter.newLine();
+ }
}
}
catch (IOException e) {