$title magpie_30
* 3 resolution, 2178 cells, 11 time slices until 2095
* Regional feed grain balances
* Water constraints: requirements for irrig. vs. discharge
* Cost structures to be linked to GDP?
* No climate change effects on yields
* Endogenous TC
* Status: 1 Feb 2008

$offupper
$offsymxref
$offsymlist
$offinclude
$offlisting


************
*Input data
************

sets
   t 10-year time periods
    /1995,2005,2015,2025,2035,2045,2055,2065,2075,2085,2095/
   ts(t) 20-year time steps
    /1995,2015,2035,2055,2075,2095/
   ts_100(t)
    /2005,2015,2025,2035,2045,2055,2065,2075,2085,2095/
   ts_50(t)
    /2005,2015,2025,2035,2045,2055/

   i economic regions /AFR,CPA,EUR,FSU,LAM,MEA,NAM,PAO,PAS,SAS/
   j number of LPJ cells /1*2178  /
   cell(i,j) number of LPJ cells per region i
       / AFR.1*262
         CPA.263*405
         EUR.406*580
         FSU.581*1042
         LAM.1043*1308
         MEA.1309*1438
         NAM.1439*1866
         PAO.1867*2002
         PAS.2003*2106
         SAS.2107*2178  /

   k activities
       / tece_food,tece_feed,maiz_food,maiz_feed,trce_food,trce_feed,rice_pro,
         soybean,rapeseed,groundnut,sunflower,oilcrop_o,
         puls_pro,potato,cassav_sp,sugr_cane,sugr_beet,veg_fr_pro,
         foddr_c3,foddr_c4,cottn_pro,bio_fuel,
         livst_rum,livst_non,milk_pro,
         irri_earl,irri_late,
         lndcon_cr,lndcon_pa,
         labor_inp,chemi_inp,capit_inp  /
   kpr(k) production activities
       / tece_food,tece_feed,maiz_food,maiz_feed,trce_food,trce_feed,rice_pro,
         soybean,rapeseed,groundnut,sunflower,oilcrop_o,
         puls_pro,potato,cassav_sp,sugr_cane,sugr_beet,veg_fr_pro,
         foddr_c3,foddr_c4,cottn_pro,bio_fuel,
         livst_rum,livst_non,milk_pro  /
   kcr(k) crop activities
       / tece_food,tece_feed,maiz_food,maiz_feed,trce_food,trce_feed,rice_pro,
         soybean,rapeseed,groundnut,sunflower,oilcrop_o,
         puls_pro,potato,cassav_sp,sugr_cane,sugr_beet,veg_fr_pro,
         foddr_c3,foddr_c4,cottn_pro,bio_fuel  /
   kpr_cr(kpr) /tece_food,tece_feed,maiz_food,maiz_feed,
         trce_food,trce_feed,rice_pro,
         soybean,rapeseed,groundnut,sunflower,oilcrop_o,
         puls_pro,potato,cassav_sp,sugr_cane,sugr_beet,veg_fr_pro,
         foddr_c3,foddr_c4,cottn_pro,bio_fuel /
   kce(k) cereals activities
       / tece_food,tece_feed,maiz_food,maiz_feed,trce_food,trce_feed /
   kpr_ce(kpr) / tece_food,tece_feed,maiz_food,maiz_feed,trce_food,trce_feed /
   koi(k) oilcrop activities / soybean,rapeseed,groundnut,sunflower,oilcrop_o /
   kpr_oi(kpr) / soybean,rapeseed,groundnut,sunflower,oilcrop_o /
   krt(k) root_tuber activities /potato,cassav_sp /
   kpr_rt(kpr) /potato,cassav_sp /
   ksu(k) sugar activities / sugr_cane, sugr_beet /
   kpr_su(kpr) / sugr_cane, sugr_beet /
   kfe(kpr) feed crop activities /tece_feed,maiz_feed,trce_feed/
   kfd(kpr) fodder crop activities /foddr_c3,foddr_c4/
   kli(kpr) livestock activities /livst_rum,livst_non,milk_pro  /
*  kpr_li(kpr) /livst_rum,livst_non,milk_pro  /
   kir(k) irrigation activities /irri_earl,irri_late  /
   klc(k) land conversion activities /lndcon_cr,lndcon_pa  /
   kin(k) input purchase activities /labor_inp,chemi_inp,capit_inp  /

   inp factors of production /labor,chemicals,capit_oth,crop_land,pasture  /
   var(inp) variable factors of production /labor,chemicals,capit_oth  /
   fix(inp) fixed factors of production /crop_land,pasture  /
   land land types /crop,past,non_ag  /

   dem demand types
   / cereals,rice_con,oilcrops,puls_con,roots_tub,sugar,veg_fr_con,
     meat_rumi,meat_nonr,milk_con,
     fiber, biof_ener  /
   food(dem) food demand types
   / cereals,rice_con,oilcrops,puls_con,roots_tub,sugar,veg_fr_con,
     meat_rumi,meat_nonr,milk_con  /
   animal(dem) animal-based food demand / meat_rumi,meat_nonr,milk_con /
   seas distinctive periods of plant growth /early,late /

   cft_all LPJ crop types (non-irrigated & irrigated)
   /CGC3,CGC4,TeCe,TrRi,TeCo,TrMi,TePu,TeSb,TrMa,TeSf,TeSo,TrPe,TeRa,
    CGC3_i,CGC4_i,TeCe_i,TrRi_i,TeCo_i,TrMi_i,TePu_i,TeSb_i,TrMa_i,TeSf_i,
    TeSo_i,TrPe_i,TeRa_i /
   cft(cft_all) LPJ crop types non-irrigated
   /CGC3,CGC4,TeCe,TrRi,TeCo,TrMi,TePu,TeSb,TrMa,TeSf,TeSo,TrPe,TeRa  /
   cft2(cft_all) LPJ crop types irrigated
   / CGC3_i,CGC4_i,TeCe_i,TrRi_i,TeCo_i,TrMi_i,TePu_i,TeSb_i,TrMa_i,TeSf_i,
     TeSo_i,TrPe_i,TeRa_i /
;

*Key parameters during model runs

scalars z count for timesteps                      / 2    /
        scal_lc_st land conversion start rate      / 0.001  /
        scal_lc_gr land conversion growth rate     / 0.001 /
        scal_tb trade balance minimum reduction
        scal_tb_st trade balance reduction start   / 0.95  /
        scal_tb_gr trade balance reduction growth  / 0.95  /
        wat_tc share of water-saving tech. chnge.  / 0.5  /
*        wat_tc: 0 = no water saving; 1 = full water saving
        disch_red reduction factor for runoff      / 0.8  /
        irr_eff irrigation efficiency              / 0.4  /
        num_cel total number of cells

;

*Read LPJ output files

parameters
 cel_siz(i,j) Size of cells in km3 (to be converted into mio ha) [LPJ]
 /
$ondelim
$include "./output_lpj/area.csv"
$offdelim
 /

 crop_shr_1995(i,j) Share of crop land in cells (convert % into decimal) [LPJ]
 /
$ondelim
$include "./output_lpj/agg_crop.csv"
$offdelim
 /

 irrig_shr(i,j) Share of irrigated land in cells (%) [LPJ]
 /
$ondelim
$include "./output_lpj/agg_irrig.csv"
$offdelim
 /

 past_shr(i,j) Share of pasture in cells (% - convert into decimal)[LPJ]
 /
$ondelim
$include "./output_lpj/agg_past.csv"
$offdelim
 /

 pop95(i,j) Gridded population 1995 [LPJ]
 /
$ondelim
$include "./output_lpj/pop95.csv"
$offdelim
 /

 precip_shr(i,j) Share of precipitation (early) [LPJ]
 /
$ondelim
$include "./input/precip_shr.csv"
$offdelim
 /

 discharge(t,i,j) discharge available for irrigation per cell (mio m3) (LPJ)
*Alternative: runoff_natveg.csv (runoff under nat. vegetation in mm p.a.)
/
$ondelim
$include "./output_lpj/discharge.csv"
$offdelim
/

;

*Additional irrigation water, received from somewhere if needed:
*(10 mio m3 per grid cell ~ 1 mm)
table irr_wat(i,j,seas) Available irrigation water in mio m3 (LPJ)
$ondelim
$include "./input/irr_wat.csv"
$offdelim       ;

table precip(i,j,t) Precipitation (mm p.a.) [multiplied by 10 to get m3_ha][LPJ]
$ondelim
$include "./output_lpj/rain.csv"
$offdelim ;

table yields_rf(t,i,j,cft) LPJ potential yields rainfed(t_DM_ha) per cell
$ondelim
$include "./output_lpj/yields_rf.csv"
$offdelim         ;

table yields_ir(t,i,j,cft) LPJ potential yields irrigated(t_DM_ha) per cell
$ondelim
$include "./output_lpj/yields_ir.csv"
$offdelim         ;

table airrig(t,i,j,cft) LPJ annual water demand for irrig. (mm p.a.) per cell
$ondelim
$include "./output_lpj/airrig.csv"
$offdelim       ;

*Read MAgPIE input files

table max_shr(i,kcr) Max shares for rotational constraints[Define?!]
$ondelim
$include "./input/max_shr_target_3.csv"
$offdelim       ;

table yld_corr(i,kpr) Yield calibration to FAO_95 (t_DM_ha) per region
$ondelim
$include "./input/yield_correct.csv"
$offdelim        ;

parameter tcc(i) inital technical change costs per region
 /
$ondelim
$include "./input/tcc_new.csv"
$offdelim
 /
;

table watreqfao(i,kpr) Average water requirements of crops per region
$ondelim
$include "./input/wat_req_fao.csv"
$offdelim          ;

table wat_shr(i,j,kpr) Share of water requirement early in the season (LPJ)
$ondelim
$include "./input/wat_shr.csv"
$offdelim   ;

table fac_req(i,inp,kpr) Factor requirements - constant in each region [GTAP]
$ondelim
$include "./input/fac_req.csv"
$offdelim       ;

table p_wat(i,seas) price for irrigation water in US$ per m3
$ondelim
$include "./input/p_wat.csv"
$offdelim       ;

table food_shr(i,dem) Food energy demand share from crop_animal products [FAO]
$ondelim
$include "./input/food_shr.csv"
$offdelim       ;

table self_food(i,dem) Self-sufficiency rate for food [FAO]
$ondelim
$include "./input/self_suff_food.csv"
$offdelim       ;

table feed_req(i,kli) livestock energy requirement (GJ per ton of output) [FAO]
$ondelim
$include "./input/feed_en_req.csv"
$offdelim       ;

table grain_shr(i,kli) share of energy requirement from feed grain [FAO]
$ondelim
$include "./input/grain_en_shr.csv"
$offdelim       ;

*parameter pop_mio(i) Population in regions in million [FAO] ;
*scalar count ;
**Calculate from gridded pop95!
*for (count = 1 to card(i),
*     pop_mio(i)$(ord(i)=count) = sum(j, pop95(i,j))/1000000 ;
*     ) ;

*Alternative:
table pop_mio(i,t) Population in million [FAO]
$ondelim
$include "./input/pop_mio.csv"
$offdelim       ;

*parameter tot_pop ;
*tot_pop = sum(i, pop_mio(i)) ;
*display pop_mio, tot_pop ;

parameters
* pop_growth(i) Population growth rate
* /
*$ondelim
*$include "./input/pop_growth_cpi_base.csv"
*$offdelim
* /

 gdp_pc(i) GDP per capita 1995 - US$ PPP (CPI baseline)
 /
$ondelim
$include "./input/gdp_pc_cpi_base.csv"
$offdelim
 /

 gdp_growth(i) GDP growth rate
 /
$ondelim
$include "./input/gdp_growth_cpi_base.csv"
$offdelim
 /

 kcal_corr(i) Food energy intake correction for total intake [FAO]
 /
$ondelim
$include "./input/kcal_corr.csv"
$offdelim
 /

 fibr_pc(i) Demand for fibre crops (e.g. cotton)(ton p.c. per year)[FAO]
 /
$ondelim
$include "./input/fibr_dem_pc.csv"
$offdelim
 /

 self_fibr(i) Self-sufficiency rate for fibre crops (e.g. cotton)[FAO]
 /
$ondelim
$include "./input/self_suff_fibr.csv"
$offdelim
 /

 ener_dem(i) Demand for energy (GJ per capita per year)[Define?!]
 /
$ondelim
$include "./input/ener_dem_zero.csv"
$offdelim
 /

 ener_growth(i) growth rate for total energy use
 /
$ondelim
$include "./input/ener_growth_cpi_base.csv"
$offdelim
 /

 ener_shr(i) Biofuel energy share in total energy demand (region)[Define?!]
 /
$ondelim
$include "./input/ener_shr_cpi_base.csv"
$offdelim
 /

 ener_shr_growth(i) Biofuel energy share growth rate (region)[Define?!]
 /
$ondelim
$include "./input/ener_shr_growth_zero.csv"
$offdelim
 /

 self_biof(i) Self-sufficiency rate for biofuels (e.g. maize)[Define?!]
 /
$ondelim
$include "./input/self_suff_biof.csv"
$offdelim
 /

 food_cont(kpr) Food energy content (GJ per ton)[FAO]
 /
$ondelim
$include "./input/food_en_cont.csv"
$offdelim
 /

 feed_cont(kpr) Feed energy content (GJ per ton)[FAO]
 /
$ondelim
$include "./input/feed_en_cont.csv"
$offdelim
 /

 cowmeat(i) Ruminant meat production from cows (GJ per ton of milk)[FAO]
 /
$ondelim
$include "./input/cow_meat_coeff.csv"
$offdelim
 /

 lcc(i) Land conversion costs (US$ per ha approx.) [GTAP]
 /
$ondelim
$include "./input/lnd_con_cst.csv"
$offdelim
 /

;

*derived parameters
parameters
 yld(i,j,kpr) Yield in ton per ha DM or ton livestock
 yld_lpj(i,j,cft) yield average (rainfed & irrigated)
 viwa(i,j,cft) average virtual water demand (only irrigated!!)
 nonag_shr(i,j) Share of non ag land (% - convert into decimal)
 crop_shr(t,i,j) Crop shares after optimization
 wat_req(i,j,kpr) Crop water requirement (m3_t_DM) per cell
 fei_pc Food energy intake (fei) (kcal_pc_day into GJ_pc_year)
 biofuel_pc(i) Biofuel consumption per capita (GJ per capita per year)
 food_deliv(i,j,kpr) Food energy delivery per activity unit
 feed_deliv(i,j,kpr) Calculate feed energy delivery per activity unit
 reg_siz(i) Size of each region (in mio. ha)
 reg_lndtyp(i,land) total regional land type areas (in mio. ha)
 watreqearl(i,j,kpr) water requirements by crops - early (m3_ha)
 watreqlate(i,j,kpr) water requirements by crops - late (m3_ha)
 land_const(i,j,land) Ag. land constraint by land type in each cell
 rot_const(i,j,kcr) Rotational constraints for each cell
 watcstearl(i,j) water constraints early (mio. m3) (for cropland + pasture!)
 watcstlate(i,j) water constraints late (mio. m3) (for cropland + pasture!)
 feedgrnbal(i,j,kli) Feed grain requirements (in GJ per ton)
 foddrbal(i,j,kli) Green fodder requirements (in GJ per ton)
 reg_dem(i,dem) regional demand (PJ or ton)
 reg_dem_net(i,dem) regional demand net of foreign trade (PJ or ton)
 glo_dem(dem) global demand (PJ or ton)
 pwatcell(i,j,seas) price of water in each cell
 maxshrcell(i,j,kcr) rotational share in each cell
 facreqcell(i,j,inp,kpr) factor requirements in each cell
 inpcost(i,j,kin) factor costs for variable inputs (0 or 1)
 lndconcost(i,j,klc) factor costs for land conversion
 scal_lc scaling factor for land conversion
 tcc_grow(i) scaled tcc according to GDP growth

*results - output
 pop_out(t,i) population in period t
 gdp_out(t,i) gdp per capital in period t
 kcal_out(t,i) kcal per capita per day in period t
 ener_reg(t,i) bio-energy demand per year in region i
 ener_glo(t) global bio-energy demand per year
 cr_lnd_new(t,i) original cropland plus land_con rel. to original (regions)
 pa_lnd_new(t,i) original pasture plus land_con rel. to original (regions)
 lu_shr_opt(t,i,j,kcr) land use shares of single crops in total area (cells)
 lu_shr_all(t,i,j) land use share of all crops in total area
 lu_opt_reg(t,i,kcr) average land use shares (regions)
 lu_tot_reg(t,i) total land use share of all crops
 cr_tot_ori(t,i) crop share in total area in previous period
 cr_tot_95(i) original crop share in total area in 1995
 cr_shr_cel(t,i,j,kcr) cropland shares of crops (cells)
 cr_shr_reg(t,i,kcr) average cropland share of crops (regions)
 cr_sum_reg(t,i) sum of crop areas (regions)
 yld_opt(t,i,kpr_cr) average crop yields (ton DM per ha)
 p_crlnd_rg(t,i) average shadow price of cropland in regions
 p_crlnd_wo(t) average world shadow price of cropland
 p_palnd_rg(t,i) average shadow price of pasture in regions
 p_palnd_wo(t) average world shadow price of pasture
 lnd_val_rg(t,i) total land value in regions
 lnd_val_wo(t) total land value - world
 p_watea_rg(t,i) average shadow price of water (early season) in regions
 p_watea_wo(t) average world shadow price for water (early season)
 p_watlt_rg(t,i) average shadow price of water (late season) in regions
 sp_wat_out(t,i,j,seas) shadow prices for water in cells
 wat_val_cl(t,i,j,seas) value of available discharge in cells
 wat_val_rg(t,i) total value of available discharge in regions
 wat_val_wo(t) total value of water - world
 meat_ru_rg(t,i) meat_rumi output regions
 meat_ru_wo(t) meat_rumi output world
 meat_nr_rg(t,i) meat_nonr output regions
 meat_nr_wo(t) meat_nonr output world
 milk_rg(t,i) milk output regions
 milk_wo(t) milk output world
 past_rg(t,i) pasture use in regions
 past_wo(t) pasture use world
 self_suff_reg(t,i) regional food supply_demand balance
 self_suff_glo(t,dem) global supply_demand balances
 trad_vol_glo(t) share of food exports in total production
 tot_cost(t) total cost of production in regions
 tc_rate(t,i) minimum technical change rates in regions per year
 tc_avrge_100(i) average tc rate over all periods
 tc_avrge_50(i) average tc rate until 2055
 yld_grow(t,i) yield update after TC in previous period
 sp_biof(t,dem) global shadow price for biofuels
;

crop_shr("1995",i,j) = crop_shr_1995(i,j) ;
yld_grow("1995",i) = 1 ;
*********
*model description
*********

variables
 x(i,j,k) activities in each cell
 yld_tc(i) yield increasing technical change (per region)
* yld_abs(i) yld_tc times 1000
 goal total cost of production
;

positive variable x ;
*positive variable yld_tc ;
*lower bound = 0 from second period onwards

equations
 cost objective function
* tcdef(i) definition of tc rate
 demand1(dem) global demand balance
 demand2(dem) global demand balance
 demand3(dem) global demand balance
 demand4(dem) global demand balance
 demand5(dem) global demand balance
 demand6(dem) global demand balance
 demand7(dem) global demand balance
 demand8(dem) global demand balance
 demand9(dem) global demand balance
 demand10(dem) global demand balance
 demand11(dem) global demand balance
 demand12(dem) global demand balance
* feedgrain global feed balance
 feedgrain(i) regional feed balances
 tradebal1(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal2(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal3(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal4(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal5(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal6(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal7(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal8(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal9(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal10(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal11(i,dem) regional trade balances i.e. minimum self-suff ratio
 tradebal12(i,dem) regional trade balances i.e. minimum self-suff ratio
 inputbal1(i,inp) regional (variable) input balances
 inputbal2(i,inp) regional (variable) input balances
 inputbal3(i,inp) regional (variable) input balances
* fodder(i,j) local green fodder balances
 fodder(i) regional green fodder balances
 cropconst(i,j) local cropland constraints
 pastconst(i,j) local pasture constraints
 lndcnvcst(i,j) local land conversion constraints
 rotation1(i,j,kcr) local rotational constraints
 rotation2(i,j,kcr) local rotational constraints
 rotation3(i,j,kcr) local rotational constraints
 rotation4(i,j,kcr) local rotational constraints
 rotation5(i,j,kcr) local rotational constraints
 rotation6(i,j,kcr) local rotational constraints
 rotation7(i,j,kcr) local rotational constraints
 rotation8(i,j,kcr) local rotational constraints
 rotation9(i,j,kcr) local rotational constraints
 rotation10(i,j,kcr) local rotational constraints
 rotation11(i,j,kcr) local rotational constraints
 rotation12(i,j,kcr) local rotational constraints
 rotation13(i,j,kcr) local rotational constraints
 rotation14(i,j,kcr) local rotational constraints
 rotation15(i,j,kcr) local rotational constraints
 rotation16(i,j,kcr) local rotational constraints
 rotation17(i,j,kcr) local rotational constraints
 rotation18(i,j,kcr) local rotational constraints
 waterconst1(i,j,seas) local seasonal water constraints
 waterconst2(i,j,seas) local seasonal water constraints
* irriconst1(i,j,seas) local seasonal irrigation constraints
* irriconst2(i,j,seas) local seasonal irrigation constraints
;

* objective function
cost .. goal =e= sum((cell,kin), x(cell,kin))
                 + sum((cell), x(cell,"irri_earl")*pwatcell(cell,"early")
                              + x(cell,"irri_late")*pwatcell(cell,"late") )
                 + sum(i, yld_tc(i)*tcc_grow(i));

*tcdef(i).. yld_abs(i) =e= 1000*yld_tc(i) ;

* global constraints
 demand1(dem) $(ord(dem)=1)
  .. sum((cell(i,j),kpr_ce), x(cell,kpr_ce)*food_deliv(cell,kpr_ce)
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand2(dem) $(ord(dem)=2)
  .. sum(cell(i,j), x(cell,"rice_pro")*food_deliv(cell,"rice_pro")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand3(dem) $(ord(dem)=3)
  .. sum((cell(i,j),kpr_oi), x(cell,kpr_oi)*food_deliv(cell,kpr_oi)
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand4(dem) $(ord(dem)=4)
  .. sum(cell(i,j), x(cell,"puls_pro")*food_deliv(cell,"puls_pro")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand5(dem) $(ord(dem)=5)
  .. sum((cell(i,j),kpr_rt), x(cell,kpr_rt)*food_deliv(cell,kpr_rt)
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand6(dem) $(ord(dem)=6)
  .. sum((cell(i,j),kpr_su), x(cell,kpr_su)*food_deliv(cell,kpr_su)
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand7(dem) $(ord(dem)=7)
  .. sum(cell(i,j), x(cell,"veg_fr_pro")*food_deliv(cell,"veg_fr_pro")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand8(dem) $(ord(dem)=8)
  .. sum(cell(i,j), x(cell,"livst_rum")*food_deliv(cell,"livst_rum")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand9(dem) $(ord(dem)=9)
  .. sum(cell(i,j), x(cell,"livst_non")*food_deliv(cell,"livst_non")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

 demand10(dem) $(ord(dem)=10)
  .. sum(cell(i,j), x(cell,"milk_pro")*food_deliv(cell,"milk_pro")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

*take only yield (for cotton)!
 demand11(dem) $(ord(dem)=11)
  .. sum(cell(i,j), x(cell,"cottn_pro")*yld(cell,"cottn_pro")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

* replace later from "food_deliv"!
 demand12(dem) $(ord(dem)=12)
  .. sum(cell(i,j), x(cell,"bio_fuel")*food_deliv(cell,"bio_fuel")
         *(1+yld_tc(i)))
     =g= glo_dem(dem) ;

* global feed balance
* feedgrain .. sum((cell,kfe), x(cell,kfe)*feed_deliv(cell,kfe))
*                  - sum((cell,kli), x(cell,kli)*feedgrnbal(cell,kli)) =g= 0 ;

*regional balances
 feedgrain(i) ..  sum((j,kfe), x(i,j,kfe)*feed_deliv(i,j,kfe)
                      *(1+yld_tc(i)))
                  - sum((j,kli), x(i,j,kli)*feedgrnbal(i,j,kli)) =g= 0 ;

* Determine different levels with scal_tb !
 tradebal1(i,dem) $(ord(dem)=1)
  .. sum((j,kpr_ce), x(i,j,kpr_ce)*food_deliv(i,j,kpr_ce)
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal2(i,dem) $(ord(dem)=2)
  .. sum(j, x(i,j,"rice_pro")* food_deliv(i,j,"rice_pro")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal3(i,dem) $(ord(dem)=3)
  .. sum((j,kpr_oi), x(i,j,kpr_oi)*food_deliv(i,j,kpr_oi)
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal4(i,dem) $(ord(dem)=4)
  .. sum(j, x(i,j,"puls_pro")*food_deliv(i,j,"puls_pro")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal5(i,dem) $(ord(dem)=5)
  .. sum((j,kpr_rt), x(i,j,kpr_rt)*food_deliv(i,j,kpr_rt)
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal6(i,dem) $(ord(dem)=6)
  .. sum((j,kpr_su), x(i,j,kpr_su)*food_deliv(i,j,kpr_su)
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal7(i,dem) $(ord(dem)=7)
  .. sum(j, x(i,j,"veg_fr_pro")*food_deliv(i,j,"veg_fr_pro")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal8(i,dem) $(ord(dem)=8)
  .. sum(j, x(i,j,"livst_rum")*food_deliv(i,j,"livst_rum")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal9(i,dem) $(ord(dem)=9)
  .. sum(j, x(i,j,"livst_non")*food_deliv(i,j,"livst_non")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal10(i,dem) $(ord(dem)=10)
  .. sum(j, x(i,j,"milk_pro")*food_deliv(i,j,"milk_pro")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal11(i,dem) $(ord(dem)=11)
  .. sum(j, x(i,j,"cottn_pro")*yld(i,j,"cottn_pro")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 tradebal12(i,dem) $(ord(dem)=12)
  .. sum(j, x(i,j,"bio_fuel")*food_deliv(i,j,"bio_fuel")
         *(1+yld_tc(i)))
     =g= reg_dem_net(i,dem) ;

 inputbal1(i,inp) $(ord(inp)=1)
  .. sum(j, x(i,j,"labor_inp")*inpcost(i,j,"labor_inp"))
     - sum((j,kpr), x(i,j,kpr)*facreqcell(i,j,inp,kpr)) =g= 0 ;

 inputbal2(i,inp) $(ord(inp)=2)
  .. sum(j, x(i,j,"chemi_inp")*inpcost(i,j,"chemi_inp"))
     - sum((j,kpr), x(i,j,kpr)*facreqcell(i,j,inp,kpr)) =g= 0 ;

* land conversion costs region-specific, later explicitly in fac_req !
* Scaled to GTAP land rents (total for regions)
 inputbal3(i,inp) $(ord(inp)=3)
  .. sum(j, x(i,j,"capit_inp")*inpcost(i,j,"capit_inp"))
     - sum((j,kpr), x(i,j,kpr)*facreqcell(i,j,inp,kpr))
     - sum((j,klc), x(i,j,klc)*lndconcost(i,j,klc)) =g= 0 ;

*Regional fodder balance (for simplicity w.g.t. local pasture use)
 fodder(i) .. sum((j,kfd) ,x(i,j,kfd)*feed_deliv(i,j,kfd)
                  *(1+yld_tc(i)))
                 - sum((j,kli), x(i,j,kli)*foddrbal(i,j,kli)) =g= 0 ;

*Local (cellular) balances

*Local fodder balance (alternative to regional balance)
* fodder(cell) .. sum(kfd, x(cell,kfd)*feed_deliv(cell,kfd))
*                 - sum(kli, x(cell,kli)*foddrbal(cell,kli)) =g= 0 ;

 cropconst(cell) .. sum(kpr, x(cell,kpr)*facreqcell(cell,"crop_land",kpr))
                   - x(cell,"lndcon_cr") =l= land_const(cell,"crop") ;
*Without lndcon_pa !
* pastconst(cell) .. sum(kpr, x(cell,kpr)*facreqcell(cell,"pasture",kpr))
*                     =l= land_const(cell,"past") ;

*Alternative: with lndcon_pa (important for future scenarios?)
 pastconst(cell(i,j)) .. sum(kpr, x(cell,kpr)*facreqcell(cell,"pasture",kpr)
                        /10)
                   - x(cell,"lndcon_pa") =l= land_const(cell,"past") ;

*Determine land conversion constraints with scal_lc !
*Idea: define absolute plus relative constraint!
 lndcnvcst(cell)
  .. sum(klc, x(cell,klc)) =l= land_const(cell,"non_ag")*scal_lc ;

* cereals
 rotation1(cell,kcr) $(ord(kcr)=1)
 .. x(cell,kcr) + x(cell,"tece_feed") - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation2(cell,kcr) $(ord(kcr)=3)
 .. x(cell,kcr) + x(cell,"maiz_feed") - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation3(cell,kcr) $(ord(kcr)=5)
 .. x(cell,kcr) + x(cell,"trce_feed") - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

* rice
 rotation4(cell,kcr) $(ord(kcr)=7)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr)  ;

* oilseeds
 rotation5(cell,kcr) $(ord(kcr)=8)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation6(cell,kcr) $(ord(kcr)=9)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation7(cell,kcr) $(ord(kcr)=10)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation8(cell,kcr) $(ord(kcr)=11)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation9(cell,kcr) $(ord(kcr)=12)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

* pulses
 rotation10(cell,kcr) $(ord(kcr)=13)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

*roots and tubers
 rotation11(cell,kcr) $(ord(kcr)=14)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation12(cell,kcr) $(ord(kcr)=15)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

* sugar crops
 rotation13(cell,kcr) $(ord(kcr)=16)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

 rotation14(cell,kcr) $(ord(kcr)=17)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

* fibre
 rotation15(cell,kcr) $(ord(kcr)=18)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
    =l= rot_const(cell,kcr) ;

*fodder crops combined
 rotation16(cell,kcr) $(ord(kcr)=19)
 .. x(cell,kcr) + x(cell,"foddr_c4")
    - x(cell,"lndcon_cr")*maxshrcell(cell,kcr) =l= rot_const(cell,kcr) ;

 rotation17(cell,kcr) $(ord(kcr)=21)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
 =l= rot_const(cell,kcr) ;

 rotation18(cell,kcr) $(ord(kcr)=22)
 .. x(cell,kcr) - x(cell,"lndcon_cr")*maxshrcell(cell,kcr)
 =l= rot_const(cell,kcr) ;

*$ontext
*irrigation water constraint - airrig vs. discharge
 waterconst1(cell(i,j),seas) $(ord(seas)=1)
  .. sum(kpr, x(cell,kpr)*watreqearl(cell,kpr)
         /(1+(yld_tc(i)*wat_tc))) =l= watcstearl(cell)  ;
 waterconst2(cell(i,j),seas) $(ord(seas)=2)
  .. sum(kpr, x(cell,kpr)*watreqlate(cell,kpr)
         /(1+(yld_tc(i)*wat_tc))) =l= watcstlate(cell)  ;
*$offtext

$ontext
*with potential additional irrigation
 waterconst1(cell,seas) $(ord(seas)=1)
  .. sum(kpr, x(cell,kpr)*watreqearl(cell,kpr)) - x(cell,"irri_earl")
     =l= watcstearl(cell)  ;
 waterconst2(cell,seas) $(ord(seas)=2)
  .. sum(kpr, x(cell,kpr)*watreqlate(cell,kpr)) - x(cell,"irri_late")
     =l= watcstlate(cell)  ;

 irriconst1(cell,seas) $(ord(seas)=1)
  .. x(cell,"irri_earl") =l= irr_wat(cell,"early") ;
 irriconst2(cell,seas) $(ord(seas)=2)
  .. x(cell,"irri_late") =l= irr_wat(cell,"late") ;
$offtext

model magpie_30 /all/ ;

option lp = cplex ;
option qcp = cplex ;
option nlp = conopt ;
option iterlim = 100000 ;
option reslim  = 100000 ;
option solprint = off;
option sysout = off;
option limcol = 0;
option limrow = 0;
option decimals = 3;
*option bratio = 1 ;

******************************
loop (t$(ord(t) le z),

tcc_grow(i) = tcc(i)*(1+gdp_growth(i))**((ord(t)-1)*10) ;

*Intermediate calculations
*Non-ag land shares
nonag_shr(cell) = 100 - crop_shr(t,cell) - (past_shr(cell) - 0.5);

*average yield and viwa irrig/non-irrig
*WITHOUT climate effects on yields
*$ontext
*Assumption: constant 1995 ratio irrig/non-irrig in each cell
yld_lpj(i,j,cft)$(cell(i,j) and crop_shr(t,i,j)>0)
       = ((irrig_shr(i,j)*yields_ir("1995",i,j,cft)
         + (crop_shr_1995(i,j) - irrig_shr(i,j))*yields_rf("1995",i,j,cft))
         /crop_shr_1995(i,j))
         *yld_grow(t,i)  ;

* viwa is here only additional water demand for irrigation
viwa(i,j,cft)$(cell(i,j) and crop_shr(t,i,j)>0 and yields_ir("1995",i,j,cft)>0)
            = ((airrig("1995",i,j,cft)*10/yields_ir("1995",i,j,cft))
              *irrig_shr(i,j)/crop_shr_1995(i,j)) ;

*water constraints early (mio. m3) (discharge in total cell, reduced by x %)
 watcstearl(i,j)$cell(i,j) = discharge("1995",i,j)*disch_red*precip_shr(i,j);
*water constraints late (mio. m3) (discharge in total cell, reduced by x %)
 watcstlate(i,j)$cell(i,j) = discharge("1995",i,j)
                             *disch_red*(1-precip_shr(i,j));
*$offtext

*ALTERNATIVE
*WITH climate effects on yields
$ontext
*Assumption: constant 1995 ratio irrig/non-irrig in each cell
yld_lpj(i,j,cft)$(cell(i,j) and crop_shr(t,i,j)>0)
       = ((irrig_shr(i,j)*yields_ir(t,i,j,cft)
         + (crop_shr_1995(i,j) - irrig_shr(i,j))*yields_rf(t,i,j,cft))
         /crop_shr_1995(i,j))
         *yld_grow(t,i)  ;

* viwa is here only additional water demand for irrigation
viwa(i,j,cft)$(cell(i,j) and crop_shr(t,i,j)>0 and yields_ir(t,i,j,cft)>0)
            = ((airrig(t,i,j,cft)*10/yields_ir(t,i,j,cft))
              *irrig_shr(i,j)/crop_shr_1995(i,j)) ;

*water constraints early (mio. m3) (discharge in total cell, reduced by x %)
 watcstearl(i,j)$cell(i,j) = discharge(t,i,j)*disch_red*precip_shr(i,j);
*water constraints late (mio. m3) (discharge in total cell, reduced by x %)
 watcstlate(i,j)$cell(i,j) = discharge(t,i,j)
                             *disch_red*(1-precip_shr(i,j));
$offtext

*Corrected yield levels
yld(i,j,"tece_food")$cell(i,j) = yld_corr(i,"tece_food")*yld_lpj(i,j,"TeCe") ;
yld(i,j,"tece_feed")$cell(i,j) = yld_corr(i,"tece_feed")*yld_lpj(i,j,"TeCe") ;
yld(i,j,"maiz_food")$cell(i,j) = yld_corr(i,"maiz_food")*yld_lpj(i,j,"TeCo") ;
yld(i,j,"maiz_feed")$cell(i,j) = yld_corr(i,"maiz_feed")*yld_lpj(i,j,"TeCo");
yld(i,j,"trce_food")$cell(i,j) = yld_corr(i,"trce_food")*yld_lpj(i,j,"TrMi");
yld(i,j,"trce_feed")$cell(i,j) = yld_corr(i,"trce_feed")*yld_lpj(i,j,"TrMi");
yld(i,j,"rice_pro")$cell(i,j) = yld_corr(i,"rice_pro")*yld_lpj(i,j,"TrRi");
yld(i,j,"soybean")$cell(i,j) = yld_corr(i,"soybean")*yld_lpj(i,j,"TeSo");
yld(i,j,"rapeseed")$cell(i,j) = yld_corr(i,"rapeseed")*yld_lpj(i,j,"TeRa");
yld(i,j,"groundnut")$cell(i,j) = yld_corr(i,"groundnut")*yld_lpj(i,j,"TrPe");
yld(i,j,"sunflower")$cell(i,j) = yld_corr(i,"sunflower")*yld_lpj(i,j,"TeSf");
* oilcrop_o -> TrPe !
yld(i,j,"oilcrop_o")$cell(i,j) = yld_corr(i,"oilcrop_o")*yld_lpj(i,j,"TrPe");
yld(i,j,"puls_pro")$cell(i,j) = yld_corr(i,"puls_pro")*yld_lpj(i,j,"TePu");
* potato -> TeSb !
yld(i,j,"potato")$cell(i,j) = yld_corr(i,"potato")*yld_lpj(i,j,"TeSb");
yld(i,j,"cassav_sp")$cell(i,j) = yld_corr(i,"cassav_sp")*yld_lpj(i,j,"TrMa");
* sugr_cane -> TrPe !
yld(i,j,"sugr_cane")$cell(i,j) = yld_corr(i,"sugr_cane")*yld_lpj(i,j,"TrPe");
yld(i,j,"sugr_beet")$cell(i,j) = yld_corr(i,"sugr_beet")*yld_lpj(i,j,"TeSb");
* veg_fr_pro -> TePu !
yld(i,j,"veg_fr_pro")$cell(i,j) = yld_corr(i,"veg_fr_pro")*yld_lpj(i,j,"TePu");
yld(i,j,"foddr_c3")$cell(i,j) = yld_corr(i,"foddr_c3")*yld_lpj(i,j,"CGC3");
yld(i,j,"foddr_c4")$cell(i,j) = yld_corr(i,"foddr_c4")*yld_lpj(i,j,"CGC4");
* cottn_pro -> TePu
yld(i,j,"cottn_pro")$cell(i,j) = yld_corr(i,"cottn_pro")*yld_lpj(i,j,"TePu");
yld(i,j,"bio_fuel")$cell(i,j) = yld_corr(i,"bio_fuel")*yld_lpj(i,j,"TeCo");

yld(i,j,"livst_rum")$cell(i,j) = yld_corr(i,"livst_rum")*yld_grow(t,i) ;
yld(i,j,"livst_non")$cell(i,j) = yld_corr(i,"livst_non")*yld_grow(t,i) ;
yld(i,j,"milk_pro")$cell(i,j) = yld_corr(i,"milk_pro")*yld_grow(t,i) ;

*Corrected water requirements
wat_req(i,j,"tece_food")$cell(i,j) = viwa(i,j,"TeCe") ;
wat_req(i,j,"tece_feed")$cell(i,j) = viwa(i,j,"TeCe") ;
wat_req(i,j,"maiz_food")$cell(i,j) = viwa(i,j,"TeCo") ;
wat_req(i,j,"maiz_feed")$cell(i,j) = viwa(i,j,"TeCo") ;
wat_req(i,j,"trce_food")$cell(i,j) = viwa(i,j,"TrMi") ;
wat_req(i,j,"trce_feed")$cell(i,j) = viwa(i,j,"TrMi") ;
wat_req(i,j,"rice_pro")$cell(i,j) = viwa(i,j,"TrRi") ;
wat_req(i,j,"soybean")$cell(i,j) = viwa(i,j,"TeSo") ;
wat_req(i,j,"rapeseed")$cell(i,j) = viwa(i,j,"TeRa") ;
wat_req(i,j,"groundnut")$cell(i,j) = viwa(i,j,"TrPe") ;
wat_req(i,j,"sunflower")$cell(i,j) = viwa(i,j,"TeSf") ;
* oilcrop_o -> TrPe
wat_req(i,j,"oilcrop_o")$cell(i,j) = viwa(i,j,"TrPe") ;
wat_req(i,j,"puls_pro")$cell(i,j) = viwa(i,j,"TePu") ;
* potato -> sugarbeet
wat_req(i,j,"potato")$cell(i,j) = viwa(i,j,"TeSb") ;
wat_req(i,j,"cassav_sp")$cell(i,j) = viwa(i,j,"TrMa") ;
* sugr_cane -> TrPe
wat_req(i,j,"sugr_cane")$cell(i,j) = viwa(i,j,"TrPe") ;
wat_req(i,j,"sugr_beet")$cell(i,j) = viwa(i,j,"TeSb") ;
* veg_fr_pro -> TePu
wat_req(i,j,"veg_fr_pro")$cell(i,j) = viwa(i,j,"TePu") ;
wat_req(i,j,"foddr_c3")$cell(i,j) = viwa(i,j,"CGC3") ;
wat_req(i,j,"foddr_c4")$cell(i,j) = viwa(i,j,"CGC4") ;
* cottn_pro -> TePu
wat_req(i,j,"cottn_pro")$cell(i,j) = viwa(i,j,"TePu") ;
wat_req(i,j,"bio_fuel")$cell(i,j) = viwa(i,j,"TeCo") ;
wat_req(i,j,"livst_rum")$cell(i,j) = watreqfao(i,"livst_rum") ;
wat_req(i,j,"livst_non")$cell(i,j) = watreqfao(i,"livst_non") ;
wat_req(i,j,"milk_pro")$cell(i,j) = watreqfao(i,"milk_pro") ;

*Food energy intake (fei) (kcal_pc_day into GJ_pc_year);
*Regression formula GDP p.c. (PPP) --> kcal
 fei_pc(i) = (kcal_corr(i)
   *802*(gdp_pc(i)*(1+gdp_growth(i))**((ord(t)-1)*10))**(0.142327))
   *4.184*365/1000000 ;

*Biofuel consumption per capita (GJ per capita per year);
 biofuel_pc(i) = ener_dem(i)*(1+ener_growth(i))**((ord(t)-1)*10)
                 *ener_shr(i)*(1+ener_shr_growth(i))**((ord(t)-1)*10) ;

*Food energy delivery per activity unit
 food_deliv(i,j,kpr)$cell(i,j) = yld(i,j,kpr)*food_cont(kpr) ;
*Calculate feed energy delivery per activity unit ;
 feed_deliv(i,j,kpr)$cell(i,j) = yld(i,j,kpr)*feed_cont(kpr) ;

* Feed energy requirements (for regional and cellular balances)
*Feed grain requirements (in GJ per ton)
*TC livestock: 50% through reduced feed demand!
 feedgrnbal(i,j,kli)$cell(i,j) = feed_req(i,kli)*grain_shr(i,kli) ;
*Green fodder requirements (in GJ per ton)
 foddrbal(i,j,kli)$cell(i,j) = feed_req(i,kli)*(1-grain_shr(i,kli)) ;

*Calculate final demand
*regional demand (PJ or ton)
 reg_dem(i,food) = pop_mio(i,t)*fei_pc(i)*food_shr(i,food) ;
 reg_dem(i,"fiber") = pop_mio(i,t)*fibr_pc(i) ;
 reg_dem(i,"biof_ener") = pop_mio(i,t)*biofuel_pc(i) ;

*trade balance reduction
 scal_tb = scal_tb_st  ;
 if ((ord(t)>1),
     scal_tb = scal_tb_gr ;
);

*regional demand net of trade (PJ or ton)
 reg_dem_net(i,food) = reg_dem(i,food)*self_food(i,food)*scal_tb ;
 reg_dem_net(i,"fiber") = reg_dem(i,"fiber")*self_fibr(i)*scal_tb ;
 reg_dem_net(i,"biof_ener") = reg_dem(i,"biof_ener")*self_biof(i)*scal_tb ;

*global demand (PJ or ton)
 glo_dem(dem) = sum(i, reg_dem(i,dem)) ;

*Total number of cells ;
 num_cel = card(j) ;

*Size of each region (in mio. ha) ;
 reg_siz(i) = sum(j, cel_siz(i,j))/10000 ;

*total regional land type areas (in mio. ha) ;
 reg_lndtyp(i,"crop") = sum(j, cel_siz(i,j)*crop_shr(t,i,j))/1000000 ;
 reg_lndtyp(i,"past") = sum(j, cel_siz(i,j)*past_shr(i,j))/1000000 ;
 reg_lndtyp(i,"non_ag") = sum(j, cel_siz(i,j)*nonag_shr(i,j))/1000000 ;

*Ag. land constraint by land type in each cell ;
*Clarify pasture constraint!
 land_const(i,j,"crop")$cell(i,j) = cel_siz(i,j)*crop_shr(t,i,j)/1000000 ;
 land_const(i,j,"past")$cell(i,j) = cel_siz(i,j)*past_shr(i,j)/1000000 ;
 land_const(i,j,"non_ag")$cell(i,j) = cel_siz(i,j)*nonag_shr(i,j)/1000000 ;

*Rotational constraints for each cell
 rot_const(i,j,kcr)$cell(i,j) = land_const(i,j,"crop")*max_shr(i,kcr) ;

*TC implementation via HI and/or LAI (defined by wat_tc)
*water requirements by crops - early (m3_ha)
 watreqearl(i,j,kpr)$cell(i,j)
   = yld(i,j,kpr)*wat_req(i,j,kpr)*wat_shr(i,j,kpr)*(1/irr_eff) ;
*water requirements by crops - late (m3_ha)
 watreqlate(i,j,kpr)$cell(i,j)
   = yld(i,j,kpr)*wat_req(i,j,kpr)*(1-wat_shr(i,j,kpr))*(1/irr_eff) ;

*price of water in each cell
 pwatcell(i,j,seas)$cell(i,j) = p_wat(i,seas);

*rotational share in each cell
 maxshrcell(i,j,kcr)$cell(i,j) = max_shr(i,kcr) ;

*factor requirements in each cell
*TC in livestock through reduced feed use incl. pasture (to be corr. later)
 facreqcell(i,j,inp,kpr)$cell(i,j) = fac_req(i,inp,kpr) ;
 facreqcell(i,j,"pasture",kpr)$cell(i,j) = fac_req(i,"pasture",kpr) ;

*factor costs for variable inputs (0 or 1)
 inpcost(i,j,kin)$cell(i,j) = 1 ;
*factor costs for land conversion
 lndconcost(i,j,klc)$cell(i,j) = lcc(i) ;

*scaling for maximum land conversion
 scal_lc = scal_lc_st ;
 if ((ord(t) > 1),
     scal_lc = scal_lc_gr ;
 );
* if ((ord(t)=3) or (ord(t)=6) or (ord(t)=9),
*     scal_lc = scal_lc_gr*20 ;
* );

*BOUNDS
yld_tc.up(i) = 1 ;

 if ((ord(t) > 1),
     yld_tc.lo(i) = 0 ;
 );

*yld_tc.fx(i)=0 ;

magpie_30.scaleopt = 1 ;

solve magpie_30 USING nlp MINIMIZING goal ;

*outputs
 pop_out(t,i)
   = pop_mio(i,t) ;
 gdp_out(t,i)
   = gdp_pc(i)*(1+gdp_growth(i))**((ord(t)-1)*10) ;
 kcal_out(t,i)
   = kcal_corr(i)
     *802*(gdp_pc(i)*(1+gdp_growth(i))**((ord(t)-1)*10))**(0.142327) ;
 ener_reg(t,i)
   = pop_mio(i,t)*biofuel_pc(i) ;

 ener_glo(t)=sum(i,ener_reg(t,i)) ;

 tot_cost(t) = goal.l ;

 tc_rate(t,i) = ((1+yld_tc.l(i))**(1/10))-1 ;

 meat_ru_rg(t,i) = sum(j, x.l(i,j,"livst_rum"))  ;
 meat_ru_wo(t) = sum(i, meat_ru_rg(t,i));
 meat_nr_rg(t,i) = sum(j, x.l(i,j,"livst_non"))  ;
 meat_nr_wo(t) = sum(i, meat_nr_rg(t,i));
 milk_rg(t,i) = sum(j, x.l(i,j,"milk_pro"))  ;
 milk_wo(t) = sum(i, milk_rg(t,i));
 past_rg(t,i) = sum((j,kpr), x.l(i,j,kpr)*facreqcell(i,j,"pasture",kpr))  ;
 past_wo(t) = sum(i, past_rg(t,i));

 cr_lnd_new(t,i)
   = sum(j, land_const(i,j,"crop")+x.l(i,j,"lndcon_cr"))/
     sum(j, land_const(i,j,"crop")) ;
 pa_lnd_new(t,i)
   = sum(j, land_const(i,j,"past")+x.l(i,j,"lndcon_pa"))/
     sum(j, land_const(i,j,"past")) ;

 lu_shr_opt(t,i,j,kcr)$cell(i,j) = x.l(i,j,kcr)/cel_siz(i,j)*10000 ;
 lu_shr_all(t,i,j)$cell(i,j) = sum(kcr, x.l(i,j,kcr))/cel_siz(i,j)*10000 ;
 lu_opt_reg(t,i,kcr) = sum(j, x.l(i,j,kcr))/reg_siz(i) ;
 lu_tot_reg(t,i) = sum((j,kcr), x.l(i,j,kcr))/reg_siz(i) ;

 cr_tot_ori(t,i) = reg_lndtyp(i,"crop")/reg_siz(i) ;
 cr_tot_95(i) = cr_tot_ori("1995",i);
 cr_shr_cel(t,i,j,kcr)
   $((cell(i,j)) and (land_const(i,j,"crop")+x.l(i,j,"lndcon_cr")>0))
   = x.l(i,j,kcr)/(land_const(i,j,"crop")+x.l(i,j,"lndcon_cr")) ;

 cr_sum_reg(t,i) = sum((j,kcr), x.l(i,j,kcr)) ;
 cr_shr_reg(t,i,kcr)
   = sum(j, x.l(i,j,kcr))/cr_sum_reg(t,i);

*old version (20 Jul 2006)
* cr_opt_reg(t,i,kcr)
*   = sum(j, x.l(i,j,kcr))/sum(j, land_const(i,j,"crop")+x.l(i,j,"lndcon_cr"));

 yld_opt(t,i,kpr_cr)
   $(sum(j, x.l(i,j,kpr_cr))>0)
   = sum(j, x.l(i,j,kpr_cr)*yld(i,j,kpr_cr))/sum(j, x.l(i,j,kpr_cr)) ;

*check later together with land conversion!
 p_crlnd_rg(t,i) =
   sum(j, -cropconst.m(i,j)*land_const(i,j,"crop"))/reg_lndtyp(i,"crop");
 p_crlnd_wo(t) = sum((i,j), -cropconst.m(i,j)*land_const(i,j,"crop"))
                 /sum(i, reg_lndtyp(i,"crop"));

 p_palnd_rg(t,i) =
   sum(j, -pastconst.m(i,j)*land_const(i,j,"past"))/reg_lndtyp(i,"past");
 p_palnd_wo(t) = sum((i,j), -pastconst.m(i,j)*land_const(i,j,"past"))
                 /sum(i, reg_lndtyp(i,"past"));

 lnd_val_rg(t,i) = sum(j, -cropconst.m(i,j)*land_const(i,j,"crop"))
                   + sum(j, -pastconst.m(i,j)*land_const(i,j,"past")) ;
 lnd_val_wo(t) = sum(i, lnd_val_rg(t,i)) ;

 sp_wat_out(t,i,j,"early") = -waterconst1.m(i,j,"early");
* sp_wat_out(t,i,j,"late")  = -waterconst2.m(i,j,"late");
 p_watea_rg(t,i) =
   sum(j, -waterconst1.m(i,j,"early")*watcstearl(i,j))/sum(j, watcstearl(i,j));
* p_watlt_rg(t,i) =
*   sum(j, waterconst2.m(i,j,"late")*watcstlate(i,j))/sum(j, watcstlate(i,j));
 p_watea_wo(t) = sum((i,j), -waterconst1.m(i,j,"early")*watcstearl(i,j))
                 /sum((i,j), watcstearl(i,j));

 wat_val_cl(t,i,j,"early")$(watcstearl(i,j)>0)
     = -waterconst1.m(i,j,"early") ;
 wat_val_rg(t,i) = sum(j, -waterconst1.m(i,j,"early")*watcstearl(i,j)) ;
 wat_val_wo(t) = sum(i, wat_val_rg(t,i)) ;

 self_suff_reg(t,i) = sum((j,kpr),
                       x.l(i,j,kpr)*food_deliv(i,j,kpr)*(1+yld_tc.l(i)))
                      /sum(food, reg_dem(i,food)) ;

 self_suff_glo(t,dem)$(ord(dem)=1) = demand1.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=2) = demand2.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=3) = demand3.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=4) = demand4.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=5) = demand5.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=6) = demand6.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=7) = demand7.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=8) = demand8.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=9) = demand9.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=10) = demand10.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=11) = demand11.l(dem)/glo_dem(dem);
 self_suff_glo(t,dem)$(ord(dem)=12)
               = demand12.l(dem)/glo_dem(dem)$(glo_dem(dem)>0);

 trad_vol_glo(t)
    = sum(i, ((self_suff_reg(t,i)-1)*sum(food, reg_dem(i,food)))
       $(self_suff_reg(t,i)>1))
      /sum((i,food), reg_dem(i,food)) ;

 sp_biof(t,dem)$(ord(dem)=12) = demand12.m(dem) ;

 crop_shr(t+1,i,j) = lu_shr_all(t,i,j)*100 ;
 yld_grow(t+1,i) = yld_grow(t,i)*(1+yld_tc.l(i))   ;

*display demand1.m, demand2.m, demand3.m, demand4.m, demand5.m,
*        demand6.m, demand7.m, demand8.m, demand9.m, demand10.m,
*        demand11.m, demand12.m ;

*display tradebal1.m, tradebal2.m, tradebal3.m, tradebal4.m, tradebal5.m,
*        tradebal6.m, tradebal7.m, tradebal8.m, tradebal9.m, tradebal10.m,
*        tradebal11.m ;

);
*loop close

tc_avrge_100(i) = (prod(t$ts_100(t), 1+tc_rate(t,i))**(1/10))-1 ;
tc_avrge_50(i) = (prod(t$ts_50(t), 1+tc_rate(t,i))**(1/6))-1 ;

display
tot_cost,
pop_out,
gdp_out,
kcal_out,
ener_reg,
ener_glo,
sp_biof,
glo_dem,
self_suff_glo,
self_suff_reg,
trad_vol_glo,
*meat_ru_rg,
*meat_ru_wo,
*meat_nr_rg,
*meat_nr_wo,
*milk_rg,
*milk_wo,
*lu_opt_reg,
p_watea_wo,
p_watea_rg,
*p_watlt_rg,
wat_val_rg,
wat_val_wo,
p_palnd_rg,
p_crlnd_rg,
lnd_val_rg,
lnd_val_wo,
*cr_shr_reg,
*yld_opt,
pa_lnd_new,
past_rg,
past_wo,
cr_lnd_new,
cr_tot_95,
lu_tot_reg,
tc_rate,
tc_avrge_100,
tc_avrge_50
;

file tech_change /"./output/tech_change.csv" / ;
tech_change.pc=5 ;
tech_change.pw=1000 ;
tech_change.nd=4;
put tech_change;
put "Tech Change Avrge (regions)", "2055", "2095" ;
loop((i),
     put / i.tl ;
     put tc_avrge_50(i) ;
     put tc_avrge_100(i) ;
);
putclose;

file lu_cell /"./output/lu_cell.csv" /  ;
lu_cell.pc=5 ;
lu_cell.pw=1000 ;
lu_cell.nd=4;
put lu_cell;
put "LU share (cells)", " ", " " ;
loop(kcr, put kcr.tl) ;
loop((t,cell(i,j)),
     put / t.tl, i.tl, j.tl ;
     loop(kcr, put lu_shr_opt(t,i,j,kcr)) ;
);
putclose;

file wat_val_cel /"./output/wat_val_cel.csv" /  ;
wat_val_cel.pc=5 ;
wat_val_cel.pw=1000 ;
wat_val_cel.nd=4;
put wat_val_cel;
put "Water value (cells)", " ", " " ;
loop(seas, put seas.tl) ;
loop((t,i,j)$cell(i,j),
     put / t.tl, i.tl, j.tl ;
     loop (seas, put wat_val_cl(t,i,j,seas)) ;
);
putclose;

*file sp_wat_cell /"./output/sp_wat_cell.csv" /  ;
*sp_wat_cell.pc=5 ;
*sp_wat_cell.pw=1000 ;
*sp_wat_cell.nd=4;
*put sp_wat_cell;
*put "Shadow price water (cells)", " ", " " ;
*loop(seas, put seas.tl) ;
*loop((t,i,j)$cell(i,j),
*     put / t.tl, i.tl, j.tl ;
*     loop (seas, put sp_wat_out(t,i,j,seas)) ;
*);
*putclose;

file lu_cr_reg /"./output/lu_reg.csv" /  ;
lu_cr_reg.pc=5 ;
lu_cr_reg.pw=1000 ;
lu_cr_reg.nd=4;
put lu_cr_reg;
put "LU share (regions)", " " ;
loop(kcr, put kcr.tl) ;
loop((t,i),
     put / t.tl, i.tl ;
     loop(kcr, put cr_shr_reg(t,i,kcr)) ;
);
putclose;

file lu_all_reg /"./output/lu_all_reg.csv" /  ;
lu_all_reg.pc=5 ;
lu_all_reg.pw=1000 ;
lu_all_reg.nd=4;
put lu_all_reg;
put "Land use share sum of crops (regions)", "1995_fao" ;
loop(t, put t.tl) ;
loop((i),
     put / i.tl, cr_tot_95(i) ;
     loop(t, put lu_tot_reg(t,i)) ;
);
putclose;

file yld_reg /"./output/yld_reg.csv" /  ;
yld_reg.pc=5 ;
yld_reg.pw=1000 ;
yld_reg.nd=4;
put yld_reg;
put "Average crop yields (regions)", " " ;
loop(kpr_cr, put kpr_cr.tl) ;
loop((t,i),
     put / t.tl, i.tl ;
     loop(kpr_cr, put yld_opt(t,i,kpr_cr)) ;
);
putclose;

file sp_crlnd_reg /"./output/sp_crlnd_reg.csv" /  ;
sp_crlnd_reg.pc=5 ;
sp_crlnd_reg.pw=1000 ;
sp_crlnd_reg.nd=4;
put sp_crlnd_reg;
put "Shadow price cropland (regions)" ;
loop(t, put t.tl) ;
put / "World" ;
loop(t, put p_crlnd_wo(t)) ;
loop((i),
     put / i.tl ;
     loop(t, put p_crlnd_rg(t,i)) ;
);
putclose;

file sp_palnd_reg /"./output/sp_palnd_reg.csv" /  ;
sp_palnd_reg.pc=5 ;
sp_palnd_reg.pw=1000 ;
sp_palnd_reg.nd=4;
put sp_palnd_reg;
put "Shadow price pasture (regions)" ;
loop(t, put t.tl) ;
put / "World" ;
loop(t, put p_palnd_wo(t)) ;
loop((i),
     put / i.tl ;
     loop(t, put p_palnd_rg(t,i)) ;
);
putclose;

file sp_wat_reg /"./output/sp_wat_reg.csv" /  ;
sp_wat_reg.pc=5 ;
sp_wat_reg.pw=1000 ;
sp_wat_reg.nd=4;
put sp_wat_reg;
put "Shadow price water (regions, early)" ;
loop(t, put t.tl) ;
put / "World" ;
loop(t, put p_watea_wo(t)) ;
loop((i),
     put / i.tl ;
     loop(t, put p_watea_rg(t,i)) ;
);
putclose;

file pop /"./output/pop.csv" /  ;
pop.pc=5 ;
pop.pw=1000 ;
pop.nd=4;
put pop;
put "Population (regions)" ;
loop(t, put t.tl) ;
loop((i),
     put / i.tl ;
     loop(t, put pop_out(t,i)) ;
);
putclose;

file gdp /"./output/gdp.csv" /  ;
gdp.pc=5 ;
gdp.pw=1000 ;
gdp.nd=4;
put gdp;
put "GDP per capita (regions)" ;
loop(t, put t.tl) ;
loop((i),
     put / i.tl ;
     loop(t, put gdp_out(t,i)) ;
);
putclose;

file kcal /"./output/kcal.csv" /  ;
kcal.pc=5 ;
kcal.pw=1000 ;
kcal.nd=4;
put kcal;
put "Food availability (kcal/capita/day, regions)" ;
loop(t, put t.tl) ;
loop((i),
     put / i.tl ;
     loop(t, put kcal_out(t,i)) ;
);
putclose;

file tot_values /"./output/tot_val.csv" /  ;
tot_values.pc=5 ;
tot_values.pw=1000 ;
tot_values.nd=4;
put tot_values ;
put " " ;
loop(t, put t.tl) ;
put / "Total cost of production (mio US$)" ;
loop(t, put tot_cost(t)) ;
put / ;
put / "Total value of land" ;
loop((i),
     put / i.tl ;
     loop(t, put lnd_val_rg(t,i)) ;
);
put / ;
put / "Total value of water" ;
loop((i),
     put / i.tl ;
     loop(t, put wat_val_rg(t,i)) ;
);
putclose;

***********************************************************
***********************************************************

This document contains some general information about the magpie-sourcecode. 
There are some restrictions concerning the code-structure that should help
to understand the code much faster. These restrictions will be explained 
in the following:

###########################
### 1.General structure ###
###########################
The Magpie-sourcecode consists of 1 main-file (magpie.gms in the main folder) 
and 6 sub-files (located in ./sourcecode).

### 1.1 magpie.gms ###
Main file. To run magpie you have to compile this file. All other 
sourcecode-files are included here. Furthermore one can configure the
behaviour of magpie and its outputs in magpie.lst.

### 1.2 sets.gms ###
Here one finds the definitions of all used sets and their relations
between each other.

### 1.3 load_input_files.gms ###
This file contains the include statements for all external data used in magpie.

### 1.4 definitions.gms ###
Definitions.gms contains all definitions of parameters, tables, variables and 
equations used in the following. Also one finds here some short explanations
of each object.

### 1.5 constraints.gms ###
Constraints.gms adheres all constraints used for the optimization process
separated in global, regional and cellular ones.

### 1.6 calculations.gms ###
Here one finds all calculations separated in three categories:
   
    (1.5.1) Preprocessing
            This part contains calculations that can be done before the 
            optimization process is started. Hence one finds here 
            conditioning and generation of input data.
            Everything in this section CAN INFLUENCE but IS NOT 
            INFLUENCED by the optimization process

    (1.5.2) Optimization process
            All calculations in this part are directly connected to the
            optimization process. Hence they use data provided by 
            the optimization routine and they produce data used
            by the optimization.
            So everything in this section CAN INFLUENCE and IS 
            INFLUENCED by the optimization process.

    (1.5.3) Postprocessing
            This part contains conditioning of the output data. 
            Hence everything in this section CANNOT INFLUENCE but IS 
            INFLUENCED by the optimization process.

### 1.7 export_data_to_files.gms ###
As the name already indicates one finds here the export-routine for 
all output-files produced by magpie.


#############################
### 2. Syntax declaration ###
#############################
For better readability tags are used for all parameters and variables. So every
parameter and variable has the structure:
 
            <tag>_<name>

Following tags are used at the moment:

            f: file   -   e.g. f_pop_mio
               every parameter with an file-tag contains unmodified data of an 
               external data source. These parameters are all defined in 
               load_input_files.gms

            i: independent of optimization process   -   e.g. i_glo_dem
               These parameters are independent of the optimization process.
               They are defined in definitions.gms and mostly calculated in
               the preprocessing-part of calculations.gms

            d: depending on the optimization process   -   e.g. d_watreq
               These parameters are influenced by the optimization. They are
               defined in definitions.gms and normally used in the 
               optimization-section of calculations.gms

            v: variables   -   e.g. v_goal
               Variables which values are obtained during the optimization
               process. They are defined in definitions.gms and mostly 
               used in constraints.gms

            s: scalars   -   e.g. s_max_timesteps
               These are some optional parameters that influence the behaviour
               of the optimization. They are defined and can be choosen 
               in magpie.gms. 


########################
### 3. Default Units ###
########################

Within Magpie all input data is calibrated to default units.
These units are:

	[1]        for shares
        [10^6]     for population
        [10^6 ha]  for areas
        [10^6]     dry matter ton] for yields
        [10^6 GJ]  for energies
        [10^6 US$] for values
        [10^6 m^3] e.g. for water

It is recommended to use these default units only.



#########################
### 4. PHP Extensions ###
#########################
Unfortunalety it is not possible to create sets in GAMS dynamically. Therefore 
these sets had to be calculated external. This is done with php. A php-script 
setup_resolution.php calculates these parts and adds them automatically into
the magpie-sourcecode. Affected files are magpie.gms, sets.gms and 
load_input_files.gms. These three files contain all a part which is clearly
marked as added by the php-script. Because these lines are modified 
automatically one MUST NOT MODIFY this code directly. Instead one has to use
setup_resolution.php. Further instructions can be found under
./php/tutorial.txt in the project folder.                            

q13_cost_tc(i2) ..
  v13_cost_tc(i2) =e= sum(ct, i13_land(i2) * i13_tc_factor(ct,i2)
                     * vm_tau(i2)**i13_tc_exponent(ct,i2)
                     * (1+pm_interest(i2))**15);

[...] 

*' @equations

*' ![Investment-yield ratio in relation to $\tau$-factor
*' [@dietrich_forecasting_2014]](tcc_regression.png){ width=60% }
*'
*' Relative technological change costs `v13_cost_tc` are calculated as a
*' heuristically derived power function of the land use intensity `vm_tau` for
*' the investment-yield-ratio (see figure above) multiplied by the current
*' regional crop areas `pc13_land` (taken from previous time step) and shifted
*' 15 years into the future using the region specific interest
*' rate `pm_interest`:

q13_cost_tc(i2) ..
  v13_cost_tc(i2) =e= sum(ct, pc13_land(i2) * i13_tc_factor(ct)
                     * sum(supreg(h2,i2),vm_tau(h2))**i13_tc_exponent(ct)
                     * (1+pm_interest(ct,i2))**15);

*' The shifting is performed because investments into technological change
*' require on average 15 years of research before a yield increase is achieved,
*' but the model has to see costs and benefits concurrently in order to take the
*' right investment decisions (see also @dietrich_forecasting_2014). Investment
*' costs are scaled in relation to crop area, since a wider areal coverage means
*' typically also higher variety in biophysical conditions, which would require
*' more research for the same overall intensity boost.
*'
*' In order to get the full investments required for the desired intensification
*' the relative technological change costs are multiplied with the given
*' intensification rate. These full costs are then distributed over an infinite
*' time horizon by multiplication with the interest rate `pm_interest(i)`
*' (annuity with infinite time horizon):

q13_tech_cost(i2) ..
 vm_tech_cost(i2) =e= sum(supreg(h2,i2), vm_tau(h2)/pcm_tau(h2)-1) * v13_cost_tc(i2)
                               * sum(ct,pm_interest(ct,i2)/(1+pm_interest(ct,i2)));

[...]