Annualised wood harvest.

GAMS/IntExtOpt.gms

@@ -154,7 +154,7 @@ $gdxin
        landCoverArea(land_cover, location)  land cover area in Mha
        landCoverChange(land_cover_before, land_cover_after, location) land cover change in Mha
        totalConversionCost(location)      land cover conversion cost - $1000 per ha
-       newForestPotentialHarvest(location)
+       forestRotaHarvest(location)
        naturalDeforestedArea(location)
        woodHarvestNat(location)           
        woodHarvestLuc(location)
@@ -164,7 +164,6 @@ $gdxin
        woodSold
        carbonFlux(location)               total carbon flux  - Mt C
        supply(all_types)
-       npv
 *       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;
 
@@ -202,7 +201,7 @@ $gdxin
        LAND_COVER_CHANGE_CONSTRAINT(land_cover, location) conservation of land area
        CONVERSION_COST(location)                          cost of land cover conversion
 
-       NEW_FOREST_POTENTIAL_HARVEST(location)
+       FOREST_ROTA_HARVEST_CALC(location)
        WOOD_HARVEST_NATURAL_CALC(location)
        WOOD_HARVEST_NATURAL_CONSTRAINT(location)
        WOOD_HARVEST_LUC_CALC(location)
@@ -214,7 +213,7 @@ $gdxin
  
        CARBON_FLUX_CALC(location)                         calc carbon flux      
 
-       NPV_CALC         Net Present Value function;
+       COST_CALC         Net Present Value function;
 
 **************** Crop cost **************************
 
@@ -279,7 +278,7 @@ $gdxin
 
 ************* Forestry ***********************************
 
- NEW_FOREST_POTENTIAL_HARVEST(location) .. newForestPotentialHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location));
+ FOREST_ROTA_HARVEST_CALC(location) .. forestRotaHarvest(location) =E= sum(land_cover, landCoverChange(land_cover, 'timberForest', location) * woodYieldRota(land_cover, 'timberForest', location)) / forestRotationPeriod(location);
 
  WOOD_HARVEST_NATURAL_CALC(location) .. woodHarvestNat(location) =E= naturalDeforestedArea(location) * woodYieldLuc('natural', 'natural', location);
  
@@ -287,7 +286,7 @@ $gdxin
                                                           
  WOOD_HARVEST_LUC_CALC(location) .. woodHarvestLuc(location) =E= sum((forested, exc_natural), landCoverChange(forested, exc_natural, location) * woodYieldLuc(forested, exc_natural, location));
 
- WOOD_SUPPLY_CALC .. woodSupply =E= woodFromRotation + sum(location, woodHarvestNat(location) + woodHarvestLuc(location));
+ WOOD_SUPPLY_CALC .. woodSupply =E= sum(location, forestRotaHarvest(location) + woodHarvestNat(location) + woodHarvestLuc(location));
  
  WOOD_MARKET_CONSTRAINT .. woodSold =L= woodDemand + woodExported - woodImported;
  
@@ -305,26 +304,13 @@ $gdxin
 
 ************ Net Present Value ******************************
 
- NPV_CALC .. npv =E= {
-                        [
-                           sum(import_crop, exportPrices(import_crop) * supply(import_crop)) +
-                           sum(location, landCoverArea('natural', location)) * naturalAreaValue
-                        ] * exp(-discountRate) -
-                        sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop)))
-                     } / (1 - exp(-discountRate)) +
-                     
-                     sum(location,
-                        [
-                            (woodExportPrice - vegClearingCost) * newForestPotentialHarvest(location) * exp(-discountRate * forestRotationPeriod(location)) -
-                            landCoverArea('timberForest', location) * forestEstablishmentCost     
-                        ] / (1 - exp(-discountRate * forestRotationPeriod(location)))
-                     ) +
-                     (woodExportPrice - vegClearingCost) * woodSold +
-                     sum(location, carbonFlux(location) * carbonPrice) -
-                     sum(import_crop, importPrices(import_crop) * importAmount(import_crop)) -
-                     sum(location, totalConversionCost(location)) -
-                     woodImportPrice * woodImported -
-                     sum(location, woodHarvestLuc(location) * vegClearingCost);
+ COST_CALC .. total_cost =E= sum((crop, location), cropArea(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) * domesticPriceMarkup +
+                       sum(location, totalConversionCost(location)) +
+                       sum(import_crop, importPrices(import_crop) * importAmount(import_crop)) +
+                       woodImported * woodImportPrice -
+                       sum(import_crop, exportPrices(import_crop) * exportAmount(import_crop)) -
+                       (woodExportPrice - vegClearingCost) * woodSold;
+                       
                   
 
 
@@ -344,8 +330,8 @@ $gdxin
 
 * display landCoverChange.L
 
-* SOLVE LAND_USE USING NLP MINIMIZING total_cost;
- SOLVE LAND_USE USING NLP MAXIMIZING npv;
+ SOLVE LAND_USE USING NLP MINIMIZING total_cost;
+* SOLVE LAND_USE USING NLP MAXIMIZING npv;
 
  display previousCropArea, irrigMaxRate, otherIntensity.L, fertI.L, irrigI.L, cropArea.L;
  display supply.l, demand, ruminantFeed.l, monogastricFeed.l, exportAmount.l;

src/ac/ed/lurg/forestry/WoodYieldItem.java

@@ -29,10 +29,6 @@ public class WoodYieldItem implements RasterItem {
 			LccKey key = entry.getKey();
 			Double[] yields = entry.getValue();
 			
-			if (!landUseTiles.containsKey(key.getFromLc())) {
-				int foo = 1;
-			}
-			
 			LandCoverTile tiles = landUseTiles.get(key.getFromLc());
 			
 			double totalArea = tiles.getTotalConvertibleArea(); // Assuming no harvest from protected areas
@@ -79,6 +75,7 @@ public class WoodYieldItem implements RasterItem {
 			currentRotationHarvest = 0.0;
 		} else {
 			currentRotationHarvest += timberTiles.getConvertibleArea(optimalRotation) * yieldAtRotation;
+			/*
 			timberTiles.resetAreaForAge(optimalRotation); // WARNING: side effect
 			
 			for (int age=optimalRotation+1; age<=LandCoverTile.getMaxAgeBin(); age++) { // TODO better way of dealing with stands over rotation age?
@@ -86,6 +83,7 @@ public class WoodYieldItem implements RasterItem {
 				currentRotationHarvest += areaHarvested * yields[age];
 				timberTiles.resetAreaForAge(age, areaHarvested); 
 			}
+			*/
 		
 		}
 	}