Merge branch 'develop' of https://github.com/jansteinhauser/magpie in…

…to develop

.codeCheck

@@ -0,0 +1,6 @@
+capitalExclusionList:
+  - land
+  - iso
+  - type
+  - age
+  - k

CHANGELOG.md

@@ -4,7 +4,6 @@
 All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
-
 ## [Unreleased]
 
 ### changed
@@ -26,12 +25,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **scripts** added automatic set writer for new bioenergy realization to `start_functions`
 - **scripts** added start scripts for the GENIE project
 - **scenario_config.csv** added preset for GENIE project
+- **config** added `.codeCheck` with additonal configuration when running `gms::codeCheck`
+- **30_crop** Improved representation of cropland requiring relocation in response to introducing semi-natural habitat at the 1 km level based on high-resolution satellite imagery.
 
 ### removed
--
+- **core**    removed no longer needed set `si` Suitability classes
 
 ### fixed
--
+- **52_carbon** i52_land_carbon_sink was not identical before 2020 for different RCPs. Fixed by setting to RCPBU until the year defined in sm_fix_cc.
 
 
 ## [4.7.1] - 2024-02-28

config/default.cfg

@@ -565,7 +565,7 @@ cfg$gms$c20_scp_type <- "sugar"                         # def = sugar
 # *                     to exports
 cfg$gms$trade <- "selfsuff_reduced"             # def = selfsuff_reduced
 
-# * option for `exo` realization only: 
+# * option for `exo` realization only:
 # * fix `v21_manna_from_heaven` to zero (0) or not (1)
 # Note: Without fixing to zero, v21_manna_from_heaven might be used unnecessarily in runs started with highres.R
 cfg$gms$s21_manna_from_heaven <- 0
@@ -695,15 +695,17 @@ cfg$gms$s30_rotation_scenario_target <- 2050       # def = 2050
 cfg$gms$c30_marginal_land <- "q33_marginal"   # def = "q33_marginal"
 
 # * Share of available cropland that is withheld for maintaining semi-natural vegetation (SNV)
-# * in cropland areas, including grassland, forest and other land. For example,
+# * in each square km in cropland areas, including grassland, forest and other land. For example,
 # * a share of 0.2 corresponds to 20 % of SNV in terms of the available cropland.
-# * Accepted values are between 0 and 1
+# * The amount of cropland relocation is estimated based on external high resolution
+# * land cover information from the Copernicus Global Land Service for the year 2019.
+# * Accepted sensible values are between 0 and 0.5
 # Note: s30_snv_shr applies to countries selected in policy_countries30
 # s30_snv_shr_noselect applies to all other countries.
 cfg$gms$s30_snv_shr <- 0                # def = 0
 cfg$gms$s30_snv_shr_noselect <- 0       # def = 0
 # * Year by which SNV policy ('s30_snv_shr') is fully implemented.
-# * Start year:
+# * Start year (should be close to 2019):
 cfg$gms$s30_snv_scenario_start <- 2020       # def = 2020
 # * Target year (year when full implementation is reached):
 cfg$gms$s30_snv_scenario_target <- 2030       # def = 2030

core/macros.gms

@@ -73,21 +73,27 @@ $macro m_fillmissingyears(input,sets) loop(t_all, \
           ct_all(t_all) = no;    \
        );
 
-* macro for linear interpolation
-$macro m_linear_interpol(input,start_year,target_year,start_value,target_value) \
+*** Time interpolation
+* macro for linear time interpolation
+$macro m_linear_time_interpol(input,start_year,target_year,start_value,target_value) \
   input(t_all)$(m_year(t_all) > start_year AND m_year(t_all) < target_year) = ((m_year(t_all)-start_year) / (target_year-start_year));  \
   input(t_all) = start_value + input(t_all) * (target_value-start_value); \
   input(t_all)$(m_year(t_all) <= start_year) = start_value; \
   input(t_all)$(m_year(t_all) >= target_year) = target_value;
 
-* macro for sigmoid interpolation (S-shaped curve)
-$macro m_sigmoid_interpol(input,start_year,target_year,start_value,target_value) \
+* macro for sigmoid time interpolation (S-shaped curve)
+$macro m_sigmoid_time_interpol(input,start_year,target_year,start_value,target_value) \
   input(t_all)$(m_year(t_all) >= start_year AND m_year(t_all) <= target_year) = ((m_year(t_all)-start_year) / (target_year-start_year));  \
   input(t_all) = 1 / (1 + exp(-10*(input(t_all)-0.5))); \
   input(t_all) = start_value + input(t_all) * (target_value-start_value); \
     input(t_all)$(m_year(t_all) <= start_year) = start_value; \
     input(t_all)$(m_year(t_all) >= target_year) = target_value;
 
+*** Data interpolation
+* macro for linear cell data interpolation
+$macro m_linear_cell_data_interpol(output,x,input_x1,input_x2,input_y1,input_y2) \
+  output(j) = input_y1 + (input_y2 - input_y1)/(input_x2 - input_x1) * (x - input_x1);
+
 * macro for simple carbon stocks
 $macro m_carbon_stock(land,carbon_density,item) \
             (land(j2,item) * sum(ct,carbon_density(ct,j2,item,ag_pools)))$(sameas(stockType,"actual")) + \

core/sets.gms

@@ -257,9 +257,6 @@ sets
   forest_type forest type
          / plantations, natveg /
 
-   si Suitability classes
-        / si0, nsi0 /
-
 ***Forestry**
   ac Age classes  / ac0,ac5,ac10,ac15,ac20,ac25,ac30,ac35,ac40,ac45,ac50,
                     ac55,ac60,ac65,ac70,ac75,ac80,ac85,ac90,ac95,ac100,

literature.bib

@@ -1596,3 +1596,13 @@ @article{mishra_timbercities_2022
     keywords = {Forestry, Climate-change mitigation},
     pages = {4889},
 }
+
+@dataset{buchhorn_copernicus_2020,
+  title = {Copernicus Global Land Service: {{Land}} Cover 100m: Collection 3: Epoch 2019: {{Globe}}},
+  author = {Buchhorn, Marcel and Smets, Bruno and Bertels, Luc and Roo, Bert De and Lesiv, Myroslava and Tsendbazar, Nandin-Erdene and Herold, Martin and Fritz, Steffen},
+  date = {2020-09},
+  publisher = {{Zenodo}},
+  doi = {10.5281/zenodo.3939050},
+  url = {https://doi.org/10.5281/zenodo.3939050},
+  version = {V3.0.1}
+}

main.gms

@@ -147,44 +147,44 @@ $title magpie
 *'  * Always try to access model outputs through the corresponding magpie package instead of accessing them directly with readGDX. It cannot be guaranteed that your script will work in the future if you do otherwise (as only the corresponding magpie package will be continuously adapted to changes in the GAMS code).
 
 *##################### R SECTION START (VERSION INFO) ##########################
-* 
+*
 * Used data set: rev4.99_h12_magpie.tgz
 * md5sum: NA
 * Repository: scp://cluster.pik-potsdam.de/p/projects/rd3mod/inputdata/output
-* 
+*
 * Used data set: rev4.99_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz
 * md5sum: NA
 * Repository: scp://cluster.pik-potsdam.de/p/projects/rd3mod/inputdata/output
-* 
+*
 * Used data set: rev4.99_h12_validation.tgz
 * md5sum: NA
 * Repository: scp://cluster.pik-potsdam.de/p/projects/rd3mod/inputdata/output
-* 
+*
 * Used data set: additional_data_rev4.47.tgz
 * md5sum: NA
 * Repository: scp://cluster.pik-potsdam.de/p/projects/landuse/data/input/archive
-* 
+*
 * Used data set: calibration_H12_per_ton_fao_may22_glo_08Aug23.tgz
 * md5sum: NA
 * Repository: https://rse.pik-potsdam.de/data/magpie/public
-* 
+*
 * Low resolution: c200
 * High resolution: 0.5
-* 
+*
 * Total number of cells: 200
-* 
+*
 * Number of cells per region:
 *   CAZ  CHA  EUR  IND  JPN  LAM  MEA  NEU  OAS  REF  SSA  USA
 *    14   23   10    7    4   26   21    9   16   23   32   15
-* 
+*
 * Regionscode: 62eff8f7
-* 
+*
 * Regions data revision: 4.99
-* 
+*
 * lpj2magpie settings:
 * * LPJmL data: MRI-ESM2-0:ssp370
 * * Revision: 4.99
-* 
+*
 * aggregation settings:
 * * Input resolution: 0.5
 * * Output resolution: c200
@@ -193,10 +193,10 @@ $title magpie
 *   CAZ  CHA  EUR  IND  JPN  LAM  MEA  NEU  OAS  REF  SSA  USA
 *    14   23   10    7    4   26   21    9   16   23   32   15
 * * Call: withCallingHandlers(expr, message = messageHandler, warning = warningHandler,     error = errorHandler)
-* 
-* 
+*
+*
 * Last modification (input data): Wed Feb  7 15:35:14 2024
-* 
+*
 *###################### R SECTION END (VERSION INFO) ###########################
 
 $offupper

modules/22_land_conservation/area_based_apr22/preloop.gms

@@ -26,7 +26,7 @@ p22_country_weight(i) = sum(i_to_iso(i,iso), p22_country_dummy(iso) * i22_land_i
 
 ** Trajectory for implementation of land conservation
 * sigmoidal interpolation between 2020 and target year
-m_sigmoid_interpol(p22_conservation_fader,s22_conservation_start,s22_conservation_target,0,1);
+m_sigmoid_time_interpol(p22_conservation_fader,s22_conservation_start,s22_conservation_target,0,1);
 
 ** Initialise additional conservation area
 p22_add_consv(t,j,consv22_all,land) = 0;

modules/30_crop/endo_apr21/declarations.gms

@@ -9,7 +9,10 @@ parameters
  p30_avl_cropland(t,j)              Total available land for crop cultivation (mio. ha)
  p30_country_snv_weight(i)          SNV policy country weight per region (1)
  p30_snv_shr(t,j)                   Share of semi-natural vegetation in cropland areas (1)
- p30_country_dummy(iso)             Dummy parameter indicating whether country is affected by selected SNV policy (1)
+ i30_snv_relocation_target(j)       Overall cropland area that requires relocation due to SNV policy (mio. ha)
+ p30_snv_relocation(t,j)            Cropland area that is actually relocated during time step (mio. ha)
+ p30_max_snv_relocation(t,j)        Maximum cropland relocation during time step (mio. ha)
+ p30_country_dummy(iso)             Dummy parameter indicating whether country is affected by selected cropland policy (1)
  i30_avl_cropland_iso(iso)          Available land area for cropland at ISO level (mio. ha)
  p30_snv_scenario_fader(t_all)      SNV scenario fader (1)
  p30_rotation_scenario_fader(t_all) Crop rotation scenario fader (1)
@@ -34,6 +37,7 @@ equations
  q30_bv_per(j,potnatveg)            Biodiversity value of perennial cropland (mio. ha)
  q30_land_snv(j)                    Land constraint for the SNV policy in cropland areas (mio. ha)
  q30_crop_reg(i)                    Total regional crop production area (mio. ha)
+ q30_land_snv_trans(j)              Land transition constraint for SNV policy in cropland areas (mio. ha)
 ;
 
 *#################### R SECTION START (OUTPUT DECLARATIONS) ####################
@@ -52,6 +56,7 @@ parameters
  oq30_bv_per(t,j,potnatveg,type)          Biodiversity value of perennial cropland (mio. ha)
  oq30_land_snv(t,j,type)                  Land constraint for the SNV policy in cropland areas (mio. ha)
  oq30_crop_reg(t,i,type)                  Total regional crop production area (mio. ha)
+ oq30_land_snv_trans(t,j,type)            Land transition constraint for SNV policy in cropland areas (mio. ha)
 ;
 *##################### R SECTION END (OUTPUT DECLARATIONS) #####################
 

modules/30_crop/endo_apr21/equations.gms

@@ -32,6 +32,14 @@
             sum(ct, p30_snv_shr(ct,j2)) * vm_land(j2,"crop")
           + sum((ct,land_snv,consv_type), pm_land_conservation(ct,j2,land_snv,consv_type));
 
+*' The semi-natural vegetation constraint for cropland areas has been suggested at the per square
+*' kilometer scale. The amount of cropland requiring relocation has therefore been derived from
+*' exogeneous high-resolution land cover information from the Copernicus Global Land Service
+*' (@buchhorn_copernicus_2020).
+
+ q30_land_snv_trans(j2) ..
+         sum(land_snv, vm_lu_transitions(j2,"crop",land_snv)) =g= sum(ct, p30_snv_relocation(ct,j2));
+
 *' As additional constraints minimum and maximum rotational constraints limit
 *' the placing of crops. On the one hand, these rotational constraints reflect
 *' crop rotations limiting the share a specific crop can cover of the total area

modules/30_crop/endo_apr21/input.gms

@@ -23,8 +23,10 @@ $setglobal c30_rotation_constraints  on
 scalars
 s30_snv_shr                     Share of available cropland that is witheld for other land cover types (1) / 0 /
 s30_snv_shr_noselect            Share of available cropland that is witheld for other land cover types (1) / 0 /
-s30_snv_scenario_start          SNV scenario start year     / 2020 /
+s30_snv_scenario_start          SNV scenario start year       / 2020 /
 s30_snv_scenario_target         SNV scenario target year      / 2030 /
+s30_snv_relocation_data_x1      First reference value in SNV target cropland data (1) / 0.2 /
+s30_snv_relocation_data_x2      Second reference value in SNV target cropland data (1) / 0.5 /
 s30_rotation_scenario_start     Rotation scenario start year      / 2020 /
 s30_rotation_scenario_target    Rotation scenario target year     / 2050 /
 s30_annual_max_growth Max annual cropland growth as share of previous cropland (-) / Inf /
@@ -114,3 +116,10 @@ $include "./modules/30_crop/endo_apr21/input/avl_cropland_iso.cs3"
 $offdelim
 ;
 
+********* SNV TARGET CROPLAND *******************************************
+
+table f30_snv_target_cropland(j,relocation_target30) Cropland in 2019 requiring relocation due to SNV policy (mio. ha)
+$ondelim
+$include "./modules/30_crop/endo_apr21/input/SNVTargetCropland.cs3"
+$offdelim
+;

modules/30_crop/endo_apr21/input/files

@@ -5,3 +5,5 @@ avl_cropland.cs3
 avl_cropland_0.5.mz
 f30_croparea_w_initialisation.cs3
 avl_cropland_iso.cs3
+SNVTargetCropland.cs3
+SNVTargetCropland_0.5.mz

modules/30_crop/endo_apr21/postsolve.gms

@@ -21,6 +21,7 @@
  oq30_bv_per(t,j,potnatveg,"marginal")          = q30_bv_per.m(j,potnatveg);
  oq30_land_snv(t,j,"marginal")                  = q30_land_snv.m(j);
  oq30_crop_reg(t,i,"marginal")                  = q30_crop_reg.m(i);
+ oq30_land_snv_trans(t,j,"marginal")            = q30_land_snv_trans.m(j);
  ov_fallow(t,j,"level")                         = vm_fallow.l(j);
  ov_area(t,j,kcr,w,"level")                     = vm_area.l(j,kcr,w);
  ov30_crop_area(t,i,"level")                    = v30_crop_area.l(i);
@@ -35,6 +36,7 @@
  oq30_bv_per(t,j,potnatveg,"level")             = q30_bv_per.l(j,potnatveg);
  oq30_land_snv(t,j,"level")                     = q30_land_snv.l(j);
  oq30_crop_reg(t,i,"level")                     = q30_crop_reg.l(i);
+ oq30_land_snv_trans(t,j,"level")               = q30_land_snv_trans.l(j);
  ov_fallow(t,j,"upper")                         = vm_fallow.up(j);
  ov_area(t,j,kcr,w,"upper")                     = vm_area.up(j,kcr,w);
  ov30_crop_area(t,i,"upper")                    = v30_crop_area.up(i);
@@ -49,6 +51,7 @@
  oq30_bv_per(t,j,potnatveg,"upper")             = q30_bv_per.up(j,potnatveg);
  oq30_land_snv(t,j,"upper")                     = q30_land_snv.up(j);
  oq30_crop_reg(t,i,"upper")                     = q30_crop_reg.up(i);
+ oq30_land_snv_trans(t,j,"upper")               = q30_land_snv_trans.up(j);
  ov_fallow(t,j,"lower")                         = vm_fallow.lo(j);
  ov_area(t,j,kcr,w,"lower")                     = vm_area.lo(j,kcr,w);
  ov30_crop_area(t,i,"lower")                    = v30_crop_area.lo(i);
@@ -63,4 +66,5 @@
  oq30_bv_per(t,j,potnatveg,"lower")             = q30_bv_per.lo(j,potnatveg);
  oq30_land_snv(t,j,"lower")                     = q30_land_snv.lo(j);
  oq30_crop_reg(t,i,"lower")                     = q30_crop_reg.lo(i);
+ oq30_land_snv_trans(t,j,"lower")               = q30_land_snv_trans.lo(j);
 *##################### R SECTION END (OUTPUT DEFINITIONS) ######################

modules/30_crop/endo_apr21/preloop.gms

@@ -7,8 +7,20 @@
 
 ** Trajectory for cropland scenarios
 * sigmoidal interpolation between start year and target year
-m_sigmoid_interpol(p30_snv_scenario_fader,s30_snv_scenario_start,s30_snv_scenario_target,0,1);
-m_sigmoid_interpol(p30_rotation_scenario_fader,s30_rotation_scenario_start,s30_rotation_scenario_target,0,1);
+m_sigmoid_time_interpol(p30_snv_scenario_fader,s30_snv_scenario_start,s30_snv_scenario_target,0,1);
+m_sigmoid_time_interpol(p30_rotation_scenario_fader,s30_rotation_scenario_start,s30_rotation_scenario_target,0,1);
+
+* linear interpolation to estimate the cropland that
+* requires relocation due to SNV policy
+if (s30_snv_shr = 0,
+i30_snv_relocation_target(j) = 0;
+
+elseif s30_snv_shr <= s30_snv_relocation_data_x1,
+m_linear_cell_data_interpol(i30_snv_relocation_target, s30_snv_shr,0,s30_snv_relocation_data_x1,0, f30_snv_target_cropland(j, "SNV20TargetCropland"));
+
+elseif s30_snv_shr > s30_snv_relocation_data_x1,
+m_linear_cell_data_interpol(i30_snv_relocation_target, s30_snv_shr,s30_snv_relocation_data_x1, s30_snv_relocation_data_x2,f30_snv_target_cropland(j, "SNV20TargetCropland"), f30_snv_target_cropland(j, "SNV50TargetCropland"));
+);
 
 *due to some rounding errors the input data currently may contain in some cases
 *very small, negative numbers. These numbers have to be set to 0 as area

modules/30_crop/endo_apr21/presolve.gms

@@ -33,6 +33,20 @@ p30_snv_shr(t,j) = p30_snv_scenario_fader(t) *
   (s30_snv_shr * sum(cell(i,j), p30_country_snv_weight(i))
   + s30_snv_shr_noselect * sum(cell(i,j), 1-p30_country_snv_weight(i)));
 
+*' Cropland relocation in response to SNV policy is based on exogeneous land
+*' cover information from the Copernicus Global Land Service (@buchhorn_copernicus_2020).
+*' The rate of the policy implementation is derived based
+*' on the difference of scenario fader values in consecutive time steps
+p30_snv_relocation(t,j) = (p30_snv_scenario_fader(t) - p30_snv_scenario_fader(t-1)) *
+  (i30_snv_relocation_target(j) * sum(cell(i,j), p30_country_snv_weight(i))
+  + s30_snv_shr_noselect * sum(cell(i,j), 1-p30_country_snv_weight(i)));
+*' The following lines take care of mismatches in the input
+*' data (derived from satellite imagery from the Copernicus
+*' Global Land Service (@buchhorn_copernicus_2020)) and in
+*' cases of cropland reduction
+p30_max_snv_relocation(t,j) = p30_snv_shr(t,j) * (p30_snv_scenario_fader(t) - p30_snv_scenario_fader(t-1)) * pcm_land(j,"crop");
+p30_snv_relocation(t,j)$(p30_snv_relocation(t, j) > p30_max_snv_relocation(t,j)) = p30_max_snv_relocation(t,j);
+
 *' Area potentially available for cropping
 p30_avl_cropland(t,j) = f30_avl_cropland(j,"%c30_marginal_land%") * (1 - p30_snv_shr(t,j));
 *' @stop

modules/30_crop/endo_apr21/sets.gms

@@ -53,4 +53,8 @@ sets
 
    policy_target30 Target year for cropland policy
    / none, by2030, by2050 /
+
+   relocation_target30 Cropland requiring relocation based on different SNV targets
+   / SNV20TargetCropland, SNV50TargetCropland /
+
 ;

modules/30_crop/penalty_apr22/declarations.gms

@@ -9,7 +9,10 @@ parameters
  p30_avl_cropland(t,j)                  Total available land for crop cultivation (mio. ha)
  p30_country_snv_weight(i)              SNV policy country weight per region (1)
  p30_snv_shr(t,j)                       Share of semi-natural vegetation in cropland areas (1)
- p30_country_dummy(iso)                 Dummy parameter indicating whether country is affected by selected SNV policy (1)
+ i30_snv_relocation_target(j)           Overall cropland area that requires relocation due to SNV policy (mio. ha)
+ p30_snv_relocation(t,j)                Cropland area that is actually relocated during time step (mio. ha)
+ p30_max_snv_relocation(t,j)            Maximum cropland relocation during time step (mio. ha)
+ p30_country_dummy(iso)                 Dummy parameter indicating whether country is affected by selected cropland policy (1)
  i30_avl_cropland_iso(iso)              Available land area for cropland at ISO level (mio. ha)
  i30_rotation_incentives(t_all,rota30)  Penalty for violating rotational constraints (USD05MER per ha)
  p30_snv_scenario_fader(t_all)          SNV scenario fader (1)
@@ -39,6 +42,7 @@ equations
  q30_bv_per(j,potnatveg)              Biodiversity value of perennial cropland (mio. ha)
  q30_land_snv(j)                    Land constraint for the SNV policy in cropland areas (mio. ha)
  q30_crop_reg(i)                      Total regional crop production area (mio. ha)
+ q30_land_snv_trans(j)              Land transition constraint for SNV policy in cropland areas (mio. ha)
 ;
 
 *#################### R SECTION START (OUTPUT DECLARATIONS) ####################
@@ -61,6 +65,7 @@ parameters
  oq30_bv_per(t,j,potnatveg,type)             Biodiversity value of perennial cropland (mio. ha)
  oq30_land_snv(t,j,type)                     Land constraint for the SNV policy in cropland areas (mio. ha)
  oq30_crop_reg(t,i,type)                     Total regional crop production area (mio. ha)
+ oq30_land_snv_trans(t,j,type)               Land transition constraint for SNV policy in cropland areas (mio. ha)
 ;
 *##################### R SECTION END (OUTPUT DECLARATIONS) #####################