stable_sectoral_model_Feb_3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "The below is the python script that defines the dairy land footprint, runs LivSim, and calculates milk producitivity and\n",
    "GHG emissions from the base and mitigation scenarios described in the article. \n",
    "'''\n",
    "\n",
    "def model_in(data_wb,scenario_parameters,l,s,timestep,MC_sims,breed,cohort,subsectors,herd_size,herd_size_y_prev,herd_size_base,check_feasibility):\n",
    "\n",
    "\n",
    "    # Load data sheets from excel\n",
    "    ghg_sheet= data_wb.sheet_by_index(1)\n",
    "    uncertainty_sheet= data_wb.sheet_by_index(2)\n",
    "       \n",
    "    if (l == 'MRT'):\n",
    "        sheet = data_wb.sheet_by_index(3)\n",
    "        diet_data=(cwd+str('\\\\diet_MRT.xlsx'))\n",
    "        \n",
    "    elif (l == 'MRH'):\n",
    "        sheet = data_wb.sheet_by_index(4) \n",
    "        diet_data=(cwd+str('\\\\diet_MRH.xlsx'))\n",
    "        \n",
    "    yield_sheet=data_wb.sheet_by_index(5)\n",
    "    TLU_dat=data_wb.sheet_by_index(6)\n",
    "    feed_util_effcy=data_wb.sheet_by_index(8)\n",
    "    diet_wb = xlrd.open_workbook(diet_data)\n",
    "\n",
    "    \n",
    "# DEFINE DATA PARAMETERS/VARIABLES USED IN SIMULATIONS #\n",
    "\n",
    "    # milk yield\n",
    "    milk_yield_reg = {} \n",
    "    milk_yield_lact={}\n",
    "    \n",
    "    # feed intake dictionaries\n",
    "    grass_intake={}\n",
    "    stover_intake={}\n",
    "    maize_bran_intake={}\n",
    "    sfsm_intake={}\n",
    "    napier_intake={}\n",
    "    napier_hy_intake={}\n",
    "    pasture_intake={}\n",
    "    \n",
    "    # dietary and excretion variables\n",
    "    faecal_Nitrogen = {}\n",
    "    urinary_Nitrogen = {}\n",
    "    fraction_roughage_eaten = {}\n",
    "    fraction_conc_eaten={}\n",
    "    gross_energy = {}\n",
    "    crude_protein_intake={}\n",
    "    DMI={}\n",
    "    DMD={}\n",
    "    acid_detergent_fibre={}\n",
    "    neutral_detergent_fibre={}\n",
    "    Diet_N={}\n",
    "    ADF={}\n",
    "    NDF={}\n",
    "    ME={}\n",
    "    Ym={}\n",
    "    VS={}\n",
    "    \n",
    "    TLU_equiv={}\n",
    "    utn_efy={}\n",
    "    \n",
    "    # GHG emissions results\n",
    "    enteric_CH4={}\n",
    "    manure_CH4={}\n",
    "    \n",
    "    soil_N2O_cropland={}\n",
    "    soil_N2O_direct_cropland={}\n",
    "    soil_N2O_indirect_cropland={}\n",
    "    soil_N2O_indirect_vol_cropland={}\n",
    "    soil_N2O_indirect_leach_cropland={}\n",
    "    \n",
    "    soil_N2O_grassland={}\n",
    "    soil_N2O_direct_grassland={}\n",
    "    soil_N2O_indirect_grassland={}\n",
    "    soil_N2O_indirect_vol_grassland={}\n",
    "    soil_N2O_indirect_leach_grassland={}\n",
    "\n",
    "    manure_N2O={}\n",
    "    manure_N2O_direct={}\n",
    "    manure_N2O_vol={}\n",
    "    manure_N2O_leach={}\n",
    "\n",
    "    soil_N2O_feed={}\n",
    "    soil_N2O_feed_direct={}\n",
    "    soil_N2O_feed_adn={}\n",
    "    soil_N2O_feed_lch={}\n",
    "    \n",
    "    energy_use_CO2={}\n",
    "    \n",
    "    # as intensities  \n",
    "    enteric_CH4_intensity={}\n",
    "    manure_CH4_intensity={}\n",
    "    soil_N2O_intensity={}\n",
    "    manure_N2O_intensity={}\n",
    "    direct_emission_intensity={}\n",
    "    \n",
    "\n",
    "    # parameters and results that are random (for uncertainty analysis)\n",
    "    Ym_rd={}\n",
    "    B0_rd={}\n",
    "    MCF_rd={}\n",
    "    EF_3_FAECAL_rd={}\n",
    "    EF_3_URINARY_rd={}\n",
    "     \n",
    "    enteric_CH4_rd={}\n",
    "    manure_CH4_rd={}\n",
    "    manure_N2O_rd={}\n",
    "    manure_N2O_direct_rd={}\n",
    "    manure_N2O_leach_rd={}\n",
    "    manure_N2O_vol_rd={}\n",
    "    energy_use_CO2_rd={}\n",
    "    \n",
    "    soil_N2O_feed_direct_rd={}\n",
    "    soil_N2O_feed_adn_rd={}\n",
    "    soil_N2O_feed_lch_rd={}\n",
    "    soil_N2O_feed_rd={}\n",
    "    soil_N2O_rd={}\n",
    "    \n",
    "    grassland_area_expansion_emissions_rd={}\n",
    "    cropland_area_expansion_emissions_rd={}\n",
    "\n",
    "    herd_size_rd={}\n",
    "    direct_emission_intensity_rd={}\n",
    "    direct_emission_absolute_rd={}\n",
    "    \n",
    "    crop_primary_area_rd={}\n",
    "    feed_primary_area_rd={}\n",
    "    \n",
    "    grassland_area_expansion_emissions_rd_ss={}\n",
    "    cropland_area_expansion_emissions_rd_ss={}\n",
    "    \n",
    "    for r in range(0,MC_sims):\n",
    "        grassland_area_expansion_emissions_rd[r]=0\n",
    "        cropland_area_expansion_emissions_rd[r]=0\n",
    "        for ss in subsectors:\n",
    "            grassland_area_expansion_emissions_rd_ss[r,ss]=0\n",
    "            cropland_area_expansion_emissions_rd_ss[r,ss]=0\n",
    "    \n",
    "    RESULTS={} # dictionary to store results of the given iteration\n",
    "    \n",
    "    # land footprint dictionaries\n",
    "    feed_area_base={}\n",
    "    feed_primary_area_base={}\n",
    "    \n",
    "    feed_area={}\n",
    "    feed_primary_area={}\n",
    "    \n",
    "    total_grasslands_area={}\n",
    "    \n",
    "    crop_primary_area={}\n",
    "    crop_primary_area_base={}\n",
    "    \n",
    "    pasture_area={}\n",
    "    pasture_area_base={}\n",
    "    \n",
    "    grass_area={}\n",
    "    grass_area_base={}\n",
    "    \n",
    "    stover_area={}\n",
    "    stover_area_base={}\n",
    "    \n",
    "    napier_area={}\n",
    "    napier_area_base={}\n",
    "    \n",
    "    napier_hy_area={}\n",
    "    napier_hy_area_base={}\n",
    "    \n",
    "    maize_as_bran_area={}\n",
    "    maize_as_bran_area_base={}\n",
    "    \n",
    "    sunflower_as_meal_area={}\n",
    "    sunflower_as_meal_area_base={}\n",
    "    \n",
    "    # monthly feed parameter dictionaries\n",
    "    mo_me_short={}\n",
    "    mo_mp_short={}\n",
    "    \n",
    "    mo_dp_cp ={}\n",
    "    mo_dp_me ={}\n",
    "    \n",
    "    mo_fi_nhy = {}\n",
    "    mo_fi_np={}\n",
    "    \n",
    "    mo_fi_ns = {}\n",
    "    mo_fi_sto={}\n",
    "    mo_fi_mb = {}\n",
    "    mo_fi_sc={}\n",
    "    mo_fi_grs = {}\n",
    "    mo_fi_pt={}\n",
    "    \n",
    "    \n",
    "    count_calendar_month={} # dictionary to keep track of calendar months of year to account for seasonality\n",
    "    \n",
    "    # Populate Parameters\n",
    "    \n",
    "    # Emission factors\n",
    "    if (l == 'MRT'):\n",
    "        col_ghg=2\n",
    "    elif (l=='MRH'):\n",
    "        col_ghg=3\n",
    "            \n",
    "    # manure EFs\n",
    "    B0=ghg_sheet.cell_value(2,col_ghg)\n",
    "    MCF_LOC=ghg_sheet.cell_value(3,col_ghg)\n",
    "    MCF_IMP=ghg_sheet.cell_value(4,col_ghg)\n",
    "    EF_1=ghg_sheet.cell_value(5,col_ghg)\n",
    "    EF_3_FAECAL=ghg_sheet.cell_value(6,col_ghg)\n",
    "    EF_3_URINARY=ghg_sheet.cell_value(7,col_ghg)\n",
    "    EF_3_PRP=ghg_sheet.cell_value(8,col_ghg)\n",
    "    FRAC_MANURE_PRP_LOC=ghg_sheet.cell_value(9,col_ghg)\n",
    "    FRAC_MANURE_MMT_LOC=ghg_sheet.cell_value(10,col_ghg)\n",
    "    FRAC_MANURE_PRP_IMP=ghg_sheet.cell_value(11,col_ghg)\n",
    "    FRAC_MANURE_MMT_IMP=ghg_sheet.cell_value(12,col_ghg)\n",
    "    EF_4_LOC=ghg_sheet.cell_value(13,col_ghg)\n",
    "    EF_4_IMP=ghg_sheet.cell_value(14,col_ghg)\n",
    "    EF_5_LOC=ghg_sheet.cell_value(15,col_ghg)\n",
    "    EF_5_IMP=ghg_sheet.cell_value(16,col_ghg)\n",
    "    FRAC_vol_LOC=ghg_sheet.cell_value(17,col_ghg)\n",
    "    FRAC_vol_IMP=ghg_sheet.cell_value(18,col_ghg)\n",
    "    FRAC_leach_LOC=ghg_sheet.cell_value(19,col_ghg)\n",
    "    FRAC_leach_IMP=ghg_sheet.cell_value(20,col_ghg)\n",
    "\n",
    "    # soil related EFs\n",
    "    EF_4_adn=ghg_sheet.cell_value(21,col_ghg)\n",
    "    EF_5_lch=ghg_sheet.cell_value(22,col_ghg)\n",
    "    FRAC_gasf=ghg_sheet.cell_value(23,col_ghg)\n",
    "    FRAC_gasm=ghg_sheet.cell_value(24,col_ghg)\n",
    "    FRAC_loss_MMS=ghg_sheet.cell_value(25,col_ghg)\n",
    "   \n",
    "    # land C stock densities\n",
    "    mean_C_density=ghg_sheet.cell_value(26,col_ghg)\n",
    "    Maximal_Forest_C_Density=ghg_sheet.cell_value(27,col_ghg)\n",
    "    MEAN_C_DENSITY_GRASSLAND=ghg_sheet.cell_value(28,col_ghg)\n",
    "    luc_emission_coefficient_grass_to_crop= 11 # as described in article\n",
    "   \n",
    "    luc_emission_coefficient_total_area=mean_C_density-MEAN_C_DENSITY_GRASSLAND\n",
    "    \n",
    "    # period over which LUC emissions are amortized\n",
    "    amort_period=ghg_sheet.cell_value(29,col_ghg)\n",
    "    \n",
    "    ## allocation factors\n",
    "    allocation_factor_grass=ghg_sheet.cell_value(30,col_ghg)\n",
    "    allocation_factor_napier=ghg_sheet.cell_value(31,col_ghg)\n",
    "    allocation_factor_maize_bran=ghg_sheet.cell_value(32,col_ghg)\n",
    "    allocation_factor_sunflower_seed_meal=ghg_sheet.cell_value(33,col_ghg)\n",
    "    allocation_factor_stover=ghg_sheet.cell_value(34,col_ghg)\n",
    "    \n",
    "    # Energy related emissions \n",
    "    feed_embodied_CO2={\"maize_bran\": ghg_sheet.cell_value(35,col_ghg),\"sunflower_seed_meal\":ghg_sheet.cell_value(35,col_ghg)} # kg co2eq per Mg \n",
    "    fertilizer_embodied_CO2=ghg_sheet.cell_value(36,col_ghg)\n",
    "        \n",
    "        \n",
    "    #  livestock weights in tlu\n",
    "    TLU_equiv['local','cow']=TLU_dat.cell_value(1, 2)\n",
    "    TLU_equiv['improved','cow']=TLU_dat.cell_value(6, 2)\n",
    "    \n",
    "    TLU_equiv['local','heifer']=TLU_dat.cell_value(2, 2)\n",
    "    TLU_equiv['improved','heifer']=TLU_dat.cell_value(7, 2)\n",
    "    \n",
    "    TLU_equiv['local','bull']=TLU_dat.cell_value(4, 2)\n",
    "    TLU_equiv['improved','bull']=TLU_dat.cell_value(9, 2)\n",
    "    \n",
    "    TLU_equiv['local','fe_calf']=TLU_dat.cell_value(3, 2)\n",
    "    TLU_equiv['improved','fe_calf']=TLU_dat.cell_value(9, 2)\n",
    "    \n",
    "    TLU_equiv['local','ml_calf']=TLU_dat.cell_value(3, 2)\n",
    "    TLU_equiv['improved','ml_calf']=TLU_dat.cell_value(9, 2)\n",
    "    \n",
    "    TLU_equiv['local','juv_male']=TLU_dat.cell_value(5, 2)\n",
    "    TLU_equiv['improved','juv_male']=TLU_dat.cell_value(10, 2)\n",
    "    \n",
    "    # utilization efficiencies of biomass sources\n",
    "    utn_efy[('grass')]=feed_util_effcy.cell_value(1,1)\n",
    "    utn_efy[('pasture')]=feed_util_effcy.cell_value(2,1)\n",
    "    utn_efy[('napier')]=feed_util_effcy.cell_value(3,1)\n",
    "    utn_efy[('napier_hy')]=feed_util_effcy.cell_value(4,1)\n",
    "    utn_efy[('maize_bran')]=feed_util_effcy.cell_value(5,1)\n",
    "    utn_efy[('sunflower_seed_meal')]=feed_util_effcy.cell_value(6,1)\n",
    "    utn_efy[('stover')]=feed_util_effcy.cell_value(7,1)\n",
    "\n",
    "    # DEFINE AND INITIALIZE RESULTS DATAFRAME VARIABLES\n",
    "    RESULTS['v1_2_Milk_yield_total']=0\n",
    "    RESULTS['v1_2_Milk_yield_base_year_kg_yr']=0\n",
    "\n",
    "    RESULTS['v6_1_Enteric_CH4_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_2_Manure_CH4_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_3_Manure_N2O_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_4_Soil_N2O_cropland_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_4_Soil_N2O_grassland_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_4_2_Feed_CO2_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_5_Direct_emissions_aggregate_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']=0\n",
    "    RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']=0\n",
    "   \n",
    "    RESULTS['v9_2_1_Grasslands_area_ha_per_lps'] =0\n",
    "    RESULTS['v9_2_2_Feed_area_ha_per_lps'] =0\n",
    "    RESULTS['v9_2_3_Crop_primary_area_ha_per_lps']=0\n",
    "\n",
    "    RESULTS['v3_2_1_Feed_intake_grass_%_dm']=0\n",
    "    RESULTS['v3_2_2_Feed_intake_stover_%_dm']=0\n",
    "    RESULTS['v3_2_3_Feed_intake_napier_%_dm']=0\n",
    "    RESULTS['v3_2_4_Feed_intake_maize_bran_%_dm']=0\n",
    "    RESULTS['v3_2_5_Feed_intake_sfsm_%_dm']=0   \n",
    "\n",
    "    RESULTS['v3_2_6_Feed_intake_pasture_%_dm']=0\n",
    "    RESULTS['v3_2_6_Feed_intake_grass_kg_dm']=0\n",
    "    RESULTS['v3_2_7_Feed_intake_stover_kg_dm']=0\n",
    "    RESULTS['v3_2_8_Feed_intake_napier_kg_dm']=0\n",
    "    RESULTS['v3_2_8_Feed_intake_napier_hy_kg_dm']=0\n",
    "    RESULTS['v3_2_9_Feed_intake_maize_bran_kg_dm']=0\n",
    "    RESULTS['v3_2_10_Feed_intake_sfsm_kg_dm']=0\n",
    "    RESULTS['v3_2_11_Feed_intake_pasture_kg_dm']=0\n",
    "    \n",
    "    RESULTS['v5_6_Total_feed_primary_area_expansion']=0\n",
    "                         \n",
    "    RESULTS['v9_15_Grass_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_16_Stover_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_17_Napier_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_17_Napier_hy_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_18_Sunflower_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_19_Maize_land_use_ha_per_TLU']=0\n",
    "    RESULTS['v9_20_Pasture_land_use_ha_per_TLU']=0\n",
    "    \n",
    "    RESULTS['v8_4_Herd_size_total_1000_head_per_lps']=0\n",
    "    RESULTS['v8_6_Herd_size_total_TLU_per_lps']=0\n",
    "    \n",
    "    RESULTS['v6_0_Units_limited']=0\n",
    "\n",
    "    RESULTS['v9_4_1_Grass_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_2_Stover_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_3_Napier_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_4_Napier_hy_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_5_Sunflower_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_6_Maize_area_ha_per_lps']=0\n",
    "    RESULTS['v9_4_7_Pasture_area_ha_per_lps']=0\n",
    "         \n",
    "# Initialize simulations\n",
    "\n",
    "    for ss in subsectors:      \n",
    "        for b in breed:\n",
    "            if (herd_size[(ss,b,'herd_pop')] ==0):\n",
    "                continue\n",
    "                \n",
    "            for c in cohort:\n",
    "                \n",
    "                # initialize variables to evaluate seasonal (monthly) nutrient availability per animal (crude protein and metabolisable energy)               \n",
    "                for i in range(0,12):\n",
    "                    mo_me_short[ss,b,c,i]  = 0\n",
    "                    mo_mp_short[ss,b,c,i]  = 0\n",
    "                    mo_dp_cp[ss,b,c,i]  = 0\n",
    "                    mo_dp_me[ss,b,c,i]  = 0\n",
    "                    \n",
    "                    mo_fi_pt[ss,b,c,i]=0\n",
    "                    mo_fi_nhy[ss,b,c,i]=0\n",
    "                    mo_fi_np[ss,b,c,i]=0\n",
    "                    mo_fi_sto[ss,b,c,i]=0\n",
    "                    mo_fi_mb[ss,b,c,i]=0\n",
    "                    mo_fi_sc[ss,b,c,i]=0\n",
    "                    mo_fi_grs[ss,b,c,i]=0\n",
    "                    mo_fi_ns[ss,b,c,i]=0\n",
    "                    \n",
    "\n",
    "                # specify animal parameters for current run\n",
    "            \n",
    "                if (c == 'cow'):\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=0\n",
    "                        start_age=28\n",
    "                        max_simu_range=12*13-20\n",
    "                        if (scenario_parameters[(ss,'artificial_insemination',s)]==1):\n",
    "                            d = Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                            d.sex='f'\n",
    "                            d.parameters['potential_annual_calving_rate']['value']=1.0\n",
    "\n",
    "                        else:\n",
    "                            d = Cow(age=start_age,breed='Zebu -- TZ south high') \n",
    "                            d.parameters['potential_annual_calving_rate']['value']=0.45\n",
    "                            d.sex='f'\n",
    "\n",
    "                    else:\n",
    "                        bc_count=5\n",
    "                        start_age=20\n",
    "                        max_simu_range=12*13-20\n",
    "                        if (scenario_parameters[(ss,'artificial_insemination',s)]==1):\n",
    "                            d = Cow(breed='Improved Tanz southern highlands')\n",
    "                            d.sex='f'\n",
    "                            d.parameters['potential_annual_calving_rate']['value']=1.0\n",
    "                        else:\n",
    "                            d = Cow(breed='Improved Tanz southern highlands') \n",
    "                            d.sex='f'\n",
    "                            \n",
    "                elif (c == 'bull'):\n",
    "                    start_age=36\n",
    "                    max_simu_range=24\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=3\n",
    "                        d=Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                        d.sex='m'\n",
    "                    else:\n",
    "                        bc_count=8\n",
    "                        d=Cow(age=start_age,breed='Improved Tanz southern highlands')\n",
    "                        d.sex='m'\n",
    "                \n",
    "                elif (c == 'ml_calf'):\n",
    "                    start_age=1\n",
    "                    max_simu_range=12\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=2\n",
    "                        d=Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                        d.sex='m'\n",
    "                    else:\n",
    "                        bc_count=7\n",
    "                        d=Cow(age=start_age,breed='Improved Tanz southern highlands')\n",
    "                        d.sex='m'\n",
    "\n",
    "                elif (c == 'fe_calf'):\n",
    "                    start_age=1\n",
    "                    max_simu_range=12\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=2\n",
    "                        d=Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                        d.sex='f'\n",
    "                    else:\n",
    "                        bc_count=7\n",
    "                        d=Cow(age=start_age,breed='Improved Tanz southern highlands')\n",
    "                        d.sex='f'\n",
    "\n",
    "                elif (c == 'heifer'):\n",
    "                    start_age=12\n",
    "                    max_simu_range=36\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=1\n",
    "                        d=Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                        d.sex='f'\n",
    "                    else:\n",
    "                        bc_count=6\n",
    "                        d=Cow(age=start_age,breed='Improved Tanz southern highlands')\n",
    "                        d.sex='f'\n",
    "                \n",
    "                elif (c == 'juv_male'):\n",
    "                    start_age=12\n",
    "                    max_simu_range=24\n",
    "                    if (b == 'local'):\n",
    "                        bc_count=4\n",
    "                        d=Cow(age=start_age,breed='Zebu -- TZ south high')\n",
    "                        d.sex='m'\n",
    "                    else:\n",
    "                        bc_count=9\n",
    "                        d=Cow(age=start_age,breed='Improved Tanz southern highlands')\n",
    "                        d.sex='m'\n",
    "                      \n",
    "                # general livsim parameters     \n",
    "                d._main_variables = None \n",
    "                d._DETERMINISTIC= True\n",
    "                d.is_bull_present = 1\n",
    "                d._MAX_ITER = 20\n",
    "                \n",
    "\n",
    "                ############################### FEEDING SCRIPT ########################\n",
    "                \n",
    "                # define feeds used in simulations from library\n",
    "                feeds=['grass','stover','napier','napier_hy','maize_bran','sunflower_seed_meal','pasture']\n",
    "\n",
    "                GRASS=Feed(parameters='Native grass') \n",
    "                PASTURE=Feed(parameters='Pasture')\n",
    "                NAPIER=Feed(parameters='napier grass') \n",
    "                NAPIER_HY=Feed(parameters='napier-legume') \n",
    "                NAPIER_SILAGE=Feed(parameters='napier silage') \n",
    "                SUNFLOWER_SEED_MEAL=Feed(parameters='sunflower seed meal')\n",
    "                MAIZE_BRAN=Feed(parameters='maize bran')\n",
    "                \n",
    "                \n",
    "                # Specify the diets for cows in line with the scenario\n",
    "                feed_cow_base=diet_wb.sheet_by_index(0)\n",
    "                \n",
    "                if (b=='local'):\n",
    "                    feed_cow_scen=diet_wb.sheet_by_index(int(scenario_parameters[(ss,'diet_id_cow_local',s)]))\n",
    "   \n",
    "                elif (b=='improved'):\n",
    "                    feed_cow_scen=diet_wb.sheet_by_index(int(scenario_parameters[(ss,'diet_id_cow_improved',s)]))\n",
    "\n",
    "                # Load stover nutrient properties \n",
    "                if ((c=='cow') and  (scenario_parameters[(ss,'urea_treatment',s)]==1)):\n",
    "                    STOVER=Feed(parameters='maize stover urea treated')\n",
    "                    m_stover=\"maize stover urea treated\"\n",
    "                    \n",
    "                else:\n",
    "                    STOVER=Feed(parameters='maize stover')\n",
    "                    m_stover=\"maize stover\"\n",
    "\n",
    "                    \n",
    "                # Define feed storage object\n",
    "                feed_storage = FeedStorage()\n",
    "\n",
    "                # Cow specific feed intake variables\n",
    "                # Feed intake per day (kg dm/head/day) at specific stage of cow's production cycle (lactating, gestating, other)\n",
    "                feed_intake_cow_early_lac={}\n",
    "                feed_intake_cow_late_lac={}\n",
    "                feed_intake_cow_gest={}\n",
    "                feed_intake_cow_other={}\n",
    "                feed_intake_cow_early_lac_base={}\n",
    "                feed_intake_cow_late_lac_base={}\n",
    "                feed_intake_cow_gest_base={}\n",
    "                feed_intake_cow_other_base={}\n",
    "                \n",
    "                # Feed intake per year from specific production stages of cow (kg DM/head/yr)\n",
    "                cons_per_year_early_lac_base={}\n",
    "                cons_per_year_late_lac_base={}\n",
    "                cons_per_year_lac_base={}\n",
    "                cons_per_year_oth_base={}\n",
    "                cons_per_year_early_lac={}\n",
    "                cons_per_year_late_lac={}\n",
    "                cons_per_year_lac={}\n",
    "                cons_per_year_oth={}\n",
    "                \n",
    "                \n",
    "                # All feed intake variables\n",
    "                feed_intake={} # feed intake on an annual basis per breed/cohort (kg dm/yr/head)\n",
    "                feed_intake_base={} # reference feed intake on an annual basis per breed/cohort (kg dm/yr/head)\n",
    "\n",
    "                # Actual feed intake for this iteration\n",
    "                # Feed for cows\n",
    "                if (c =='cow'):\n",
    "                    \n",
    "                    lactation_duration=30*d.parameters['lactation_duration']['value']\n",
    "                    early_lactation_duration=30*5\n",
    "                    late_lactation_duration=lactation_duration-early_lactation_duration\n",
    "                    calving_interval=365/d.parameters['potential_annual_calving_rate']['value']\n",
    "                    FRAC_time_lactation=lactation_duration/calving_interval\n",
    "                    FRAC_time_early_lactation=early_lactation_duration/calving_interval\n",
    "                    FRAC_time_late_lactation=late_lactation_duration/calving_interval\n",
    "                    days_not_in_lac=(1-FRAC_time_lactation)*365\n",
    "\n",
    "                    # Load diets from excel\n",
    "                    if (b=='local'):\n",
    "                        breed_set=0\n",
    "                    elif (b=='improved'):\n",
    "                        breed_set=1\n",
    "                        \n",
    "                    for f in feeds:\n",
    "                            row=(subsectors.index(ss)*3)+breed_set*(9)\n",
    "                            feed_intake_cow_early_lac_base[ss,b,c,f]=feed_cow_base.cell(row+1,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_late_lac_base[ss,b,c,f]=feed_cow_base.cell(row+2,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_gest_base[ss,b,c,f]=feed_cow_base.cell(row+2,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_other_base[ss,b,c,f]=feed_cow_base.cell(row+3,3+feeds.index(f)).value/1000\n",
    "    \n",
    "                    \n",
    "                    cons_per_year_early_lac_base[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_early_lac_base[ss,b,c,'sunflower_seed_meal']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_late_lac_base[ss,b,c,'sunflower_seed_meal']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_other_base[ss,b,c,'sunflower_seed_meal']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'maize_bran']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'maize_bran']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'maize_bran']=feed_intake_cow_other_base[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'napier']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'napier']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'napier']=feed_intake_cow_other_base[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'napier_hy']=feed_intake_cow_early_lac_base[ss,b,c,'napier_hy']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'napier_hy']=feed_intake_cow_late_lac_base[ss,b,c,'napier_hy']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'napier_hy']=feed_intake_cow_other_base[ss,b,c,'napier_hy']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'grass']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'grass']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'grass']=feed_intake_cow_other_base[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'stover']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'stover']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'stover']=feed_intake_cow_other_base[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_early_lac_base[ss,b,c,'pasture']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_late_lac_base[ss,b,c,'pasture']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_oth_base[ss,b,c,'pasture']=feed_intake_cow_other_base[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    # the annual feed intake of a feed category is equal to the sum over all production stages\n",
    "                    feed_intake_base[ss,b,c,'napier']=cons_per_year_oth_base[ss,b,c,'napier']+cons_per_year_early_lac_base[ss,b,c,'napier']+cons_per_year_late_lac_base[ss,b,c,'napier']\n",
    "                    feed_intake_base[ss,b,c,'napier_hy']=cons_per_year_oth_base[ss,b,c,'napier_hy']+cons_per_year_early_lac_base[ss,b,c,'napier_hy']+cons_per_year_late_lac_base[ss,b,c,'napier_hy']\n",
    "                    feed_intake_base[ss,b,c,'maize_bran']=cons_per_year_oth_base[ss,b,c,'maize_bran']+cons_per_year_early_lac_base[ss,b,c,'maize_bran']+cons_per_year_late_lac_base[ss,b,c,'maize_bran']\n",
    "                    feed_intake_base[ss,b,c,'sunflower_seed_meal']=cons_per_year_oth_base[ss,b,c,'sunflower_seed_meal']+cons_per_year_early_lac_base[ss,b,c,'sunflower_seed_meal']+cons_per_year_late_lac_base[ss,b,c,'sunflower_seed_meal']\n",
    "                    feed_intake_base[ss,b,c,'grass']=cons_per_year_oth_base[ss,b,c,'grass']+cons_per_year_early_lac_base[ss,b,c,'grass']+cons_per_year_late_lac_base[ss,b,c,'grass']\n",
    "                    feed_intake_base[ss,b,c,'stover']=cons_per_year_oth_base[ss,b,c,'stover']+cons_per_year_early_lac_base[ss,b,c,'stover']+cons_per_year_late_lac_base[ss,b,c,'stover']\n",
    "                    feed_intake_base[ss,b,c,'pasture']=cons_per_year_oth_base[ss,b,c,'pasture']+cons_per_year_early_lac_base[ss,b,c,'pasture']+cons_per_year_late_lac_base[ss,b,c,'pasture']\n",
    "\n",
    "\n",
    "                    # if the current scenario is a feeding scenario, diet is different from base\n",
    "                    if (scenario_parameters[(ss,'diet_scenario',s)]==1):    \n",
    "                        \n",
    "                        # load diets from excel\n",
    "                        for f in feeds:\n",
    "                            row=(subsectors.index(ss)*3)+breed_set*9\n",
    "                            \n",
    "                            feed_intake_cow_early_lac[ss,b,c,f]=feed_cow_scen.cell(row+1,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_late_lac[ss,b,c,f]=feed_cow_scen.cell(row+2,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_gest[ss,b,c,f]=feed_cow_scen.cell(row+2,3+feeds.index(f)).value/1000\n",
    "                            feed_intake_cow_other[ss,b,c,f]=feed_cow_scen.cell(row+3,3+feeds.index(f)).value/1000\n",
    "\n",
    "                        \n",
    "                        \n",
    "                    else: ## if base scenario \n",
    "\n",
    "                        feed_intake_cow_early_lac[ss,b,c,'napier']=feed_intake_cow_early_lac_base[ss,b,c,'napier']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'napier']=feed_intake_cow_late_lac_base[ss,b,c,'napier']\n",
    "                        feed_intake_cow_gest[ss,b,c,'napier']=feed_intake_cow_gest_base[ss,b,c,'napier']\n",
    "                        feed_intake_cow_other[ss,b,c,'napier']=feed_intake_cow_other_base[ss,b,c,'napier']\n",
    "                            \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'napier_hy']=feed_intake_cow_early_lac_base[ss,b,c,'napier_hy']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'napier_hy']=feed_intake_cow_late_lac_base[ss,b,c,'napier_hy']\n",
    "                        feed_intake_cow_gest[ss,b,c,'napier_hy']=feed_intake_cow_gest_base[ss,b,c,'napier_hy']\n",
    "                        feed_intake_cow_other[ss,b,c,'napier_hy']=feed_intake_cow_other_base[ss,b,c,'napier_hy']\n",
    "  \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'maize_bran']=feed_intake_cow_early_lac_base[ss,b,c,'maize_bran']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'maize_bran']=feed_intake_cow_late_lac_base[ss,b,c,'maize_bran']\n",
    "                        feed_intake_cow_gest[ss,b,c,'maize_bran']=feed_intake_cow_gest_base[ss,b,c,'maize_bran']\n",
    "                        feed_intake_cow_other[ss,b,c,'maize_bran']=feed_intake_cow_other_base[ss,b,c,'maize_bran']\n",
    "                            \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_early_lac_base[ss,b,c,'sunflower_seed_meal']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_late_lac_base[ss,b,c,'sunflower_seed_meal']\n",
    "                        feed_intake_cow_gest[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_gest_base[ss,b,c,'sunflower_seed_meal']\n",
    "                        feed_intake_cow_other[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_other_base[ss,b,c,'sunflower_seed_meal']\n",
    "                        \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'grass']=feed_intake_cow_early_lac_base[ss,b,c,'grass']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'grass']=feed_intake_cow_late_lac_base[ss,b,c,'grass']\n",
    "                        feed_intake_cow_gest[ss,b,c,'grass']=feed_intake_cow_gest_base[ss,b,c,'grass']\n",
    "                        feed_intake_cow_other[ss,b,c,'grass']=feed_intake_cow_other_base[ss,b,c,'grass']\n",
    "                            \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'stover']=feed_intake_cow_early_lac_base[ss,b,c,'stover']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'stover']=feed_intake_cow_late_lac_base[ss,b,c,'stover']\n",
    "                        feed_intake_cow_gest[ss,b,c,'stover']=feed_intake_cow_gest_base[ss,b,c,'stover']\n",
    "                        feed_intake_cow_other[ss,b,c,'stover']=feed_intake_cow_other_base[ss,b,c,'stover']\n",
    "                        \n",
    "                        feed_intake_cow_early_lac[ss,b,c,'pasture']=feed_intake_cow_early_lac_base[ss,b,c,'pasture']\n",
    "                        feed_intake_cow_late_lac[ss,b,c,'pasture']=feed_intake_cow_late_lac_base[ss,b,c,'pasture']\n",
    "                        feed_intake_cow_gest[ss,b,c,'pasture']=feed_intake_cow_gest_base[ss,b,c,'pasture']\n",
    "                        feed_intake_cow_other[ss,b,c,'pasture']=feed_intake_cow_other_base[ss,b,c,'pasture']\n",
    "                           \n",
    "                    ## set annual consumption based on fraction of time in lactation\n",
    "                        \n",
    "                    cons_per_year_late_lac[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_late_lac[ss,b,c,'sunflower_seed_meal']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_early_lac[ss,b,c,'sunflower_seed_meal']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'sunflower_seed_meal']=feed_intake_cow_other[ss,b,c,'sunflower_seed_meal']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_late_lac[ss,b,c,'maize_bran']=feed_intake_cow_late_lac[ss,b,c,'maize_bran']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'maize_bran']=feed_intake_cow_early_lac[ss,b,c,'maize_bran']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'maize_bran']=feed_intake_cow_other[ss,b,c,'maize_bran']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_late_lac[ss,b,c,'napier']=feed_intake_cow_late_lac[ss,b,c,'napier']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'napier']=feed_intake_cow_early_lac[ss,b,c,'napier']*lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'napier']=feed_intake_cow_other[ss,b,c,'napier']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_late_lac[ss,b,c,'napier_hy']=feed_intake_cow_late_lac[ss,b,c,'napier_hy']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'napier_hy']=feed_intake_cow_early_lac[ss,b,c,'napier_hy']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'napier_hy']=feed_intake_cow_other[ss,b,c,'napier_hy']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_late_lac[ss,b,c,'grass']=feed_intake_cow_late_lac[ss,b,c,'grass']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'grass']=feed_intake_cow_early_lac[ss,b,c,'grass']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'grass']=feed_intake_cow_other[ss,b,c,'grass']*days_not_in_lac\n",
    "\n",
    "                    cons_per_year_late_lac[ss,b,c,'stover']=feed_intake_cow_late_lac[ss,b,c,'stover']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'stover']=feed_intake_cow_early_lac[ss,b,c,'stover']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'stover']=feed_intake_cow_other[ss,b,c,'stover']*days_not_in_lac\n",
    "                    \n",
    "                    cons_per_year_late_lac[ss,b,c,'pasture']=feed_intake_cow_late_lac[ss,b,c,'pasture']*late_lactation_duration*FRAC_time_late_lactation\n",
    "                    cons_per_year_early_lac[ss,b,c,'pasture']=feed_intake_cow_early_lac[ss,b,c,'pasture']*early_lactation_duration*FRAC_time_early_lactation\n",
    "                    cons_per_year_oth[ss,b,c,'pasture']=feed_intake_cow_other[ss,b,c,'pasture']*days_not_in_lac\n",
    "                    \n",
    "\n",
    "                    ## sum intake over production stages to get the intake on an annual basis\n",
    "                    feed_intake[ss,b,c,'napier']=cons_per_year_oth[ss,b,c,'napier']+cons_per_year_early_lac[ss,b,c,'napier']+cons_per_year_late_lac[ss,b,c,'napier']\n",
    "                    feed_intake[ss,b,c,'napier_hy']=cons_per_year_oth[ss,b,c,'napier_hy']+cons_per_year_early_lac[ss,b,c,'napier_hy']+cons_per_year_late_lac[ss,b,c,'napier_hy']\n",
    "                    feed_intake[ss,b,c,'maize_bran']=cons_per_year_oth[ss,b,c,'maize_bran']+cons_per_year_early_lac[ss,b,c,'maize_bran']+cons_per_year_late_lac[ss,b,c,'maize_bran']\n",
    "                    feed_intake[ss,b,c,'sunflower_seed_meal']=cons_per_year_oth[ss,b,c,'sunflower_seed_meal']+cons_per_year_early_lac[ss,b,c,'sunflower_seed_meal']+cons_per_year_late_lac[ss,b,c,'sunflower_seed_meal']\n",
    "                    feed_intake[ss,b,c,'grass']=cons_per_year_oth[ss,b,c,'grass']+cons_per_year_early_lac[ss,b,c,'grass']+cons_per_year_late_lac[ss,b,c,'grass']\n",
    "                    feed_intake[ss,b,c,'stover']=cons_per_year_oth[ss,b,c,'stover']+cons_per_year_early_lac[ss,b,c,'stover']+cons_per_year_late_lac[ss,b,c,'stover']\n",
    "                    feed_intake[ss,b,c,'pasture']=cons_per_year_oth[ss,b,c,'pasture']+cons_per_year_early_lac[ss,b,c,'pasture']+cons_per_year_late_lac[ss,b,c,'pasture']\n",
    "  \n",
    "\n",
    "                else: ## if not a cow, no stage specific feeding, and therefore annual intake is daily intake * 365\n",
    "                    \n",
    "                    feed_intake[ss,b,c,'grass']=365*sheet.cell_value(bc_count+1, 2)/1000\n",
    "                    feed_intake[ss,b,c,'stover']=365*sheet.cell_value(bc_count+1, 3)/1000\n",
    "                    feed_intake[ss,b,c,'napier']=365*sheet.cell_value(bc_count+1, 4)/1000 \n",
    "                    feed_intake[ss,b,c,'napier_hy']=365*sheet.cell_value(bc_count+1, 5)/1000 \n",
    "                    feed_intake[ss,b,c,'maize_bran']=365*sheet.cell_value(bc_count+1, 6)/1000 \n",
    "                    feed_intake[ss,b,c,'sunflower_seed_meal']=365*sheet.cell_value(bc_count+1, 7)/1000 \n",
    "                    feed_intake[ss,b,c,'pasture']=365*sheet.cell_value(bc_count+1, 8)/1000\n",
    "            \n",
    "\n",
    "                if (c != 'cow'):\n",
    "                    feed_intake_base[ss,b,c,'grass']=365*sheet.cell_value(bc_count+1, 2)/1000\n",
    "                    feed_intake_base[ss,b,c,'stover']=365*sheet.cell_value(bc_count+1, 3)/1000\n",
    "                    feed_intake_base[ss,b,c,'napier']=365*sheet.cell_value(bc_count+1, 4)/1000\n",
    "                    feed_intake_base[ss,b,c,'napier_hy']=365*sheet.cell_value(bc_count+1, 5)/1000 \n",
    "                    feed_intake_base[ss,b,c,'maize_bran']=365*sheet.cell_value(bc_count+1, 6)/1000 \n",
    "                    feed_intake_base[ss,b,c,'sunflower_seed_meal']=365*sheet.cell_value(bc_count+1, 7)/1000 \n",
    "                    feed_intake_base[ss,b,c,'pasture']=365*sheet.cell_value(bc_count+1, 8)/1000 \n",
    "                    \n",
    "                \n",
    "                '''\n",
    "                LAND FOOTPRINT indicators\n",
    "                \n",
    "                This section takes the feed intake for the current animal as an annual value and calculates area of land #\n",
    "                needed to produce this feed. This is done for both the reference (base) scenario and the current iteration,#\n",
    "                which is required in order to calculate LUC emissions between scenarios \n",
    "                \n",
    "                It is important to understand that in this framework, the concentrate feeds are fed according to the \n",
    "                production cycles of the cow (lactating, gestating, etc.) whereas the other feeds are specifed as annual intakes,\n",
    "                and then the intake per month depends on the seasonality parameters (e.g the variation in yields of feeds, and rationing parameters, etc.)\n",
    "                '''\n",
    "\n",
    "                # reference yield (Mg DM/ha/yr) (corrected for efficiency of harvest/grazing)\n",
    "                grass_annual_yield=(utn_efy[('grass')])*yield_sheet.cell_value(1, 2)\n",
    "                pasture_annual_yield=(utn_efy[('pasture')])*yield_sheet.cell_value(2, 2)\n",
    "                napier_annual_yield=(utn_efy[('napier')])*yield_sheet.cell_value(3, 2)\n",
    "                napier_hy_annual_yield=(utn_efy[('napier_hy')])*yield_sheet.cell_value(4, 2)\n",
    "                stover_annual_yield=(utn_efy[('stover')])*yield_sheet.cell_value(5, 2)\n",
    "                maize_annual_yield=(utn_efy[('maize_bran')])*yield_sheet.cell_value(6, 2)*(scenario_parameters[(ss,'yield_scale_factor_maize',s)])\n",
    "                sunflower_annual_yield=(utn_efy[('sunflower_seed_meal')])*yield_sheet.cell_value(7, 2)*(scenario_parameters[(ss,'yield_scale_factor_sunflower',s)])\n",
    "                \n",
    "                # reference area (ha)\n",
    "                grass_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'grass']/grass_annual_yield\n",
    "                pasture_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'pasture']/pasture_annual_yield\n",
    "                stover_area_base[ss,b,c]=(1/1000)*(feed_intake_base[ss,b,c,'stover'])/stover_annual_yield\n",
    "                sunflower_as_meal_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'sunflower_seed_meal']/(scenario_parameters[(ss,'yield_scale_factor_sunflower','BAU')]*sunflower_annual_yield)\n",
    "                maize_as_bran_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'maize_bran']/(scenario_parameters[(ss,'yield_scale_factor_maize','BAU')]*maize_annual_yield)\n",
    "                napier_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'napier']/napier_annual_yield       \n",
    "                napier_hy_area_base[ss,b,c]=(1/1000)*feed_intake_base[ss,b,c,'napier_hy']/napier_annual_yield       \n",
    "                \n",
    "                # summarize area for total feed, primary feeds, and primary crops\n",
    "                feed_area_base[(ss,b,c)]=pasture_area_base[ss,b,c]+grass_area_base[ss,b,c]+stover_area_base[ss,b,c]+sunflower_as_meal_area_base[ss,b,c]+maize_as_bran_area_base[ss,b,c]+napier_area_base[ss,b,c]+napier_hy_area_base[ss,b,c]\n",
    "                feed_primary_area_base[(ss,b,c)]= pasture_area_base[ss,b,c]+grass_area_base[ss,b,c]+sunflower_as_meal_area_base[ss,b,c]+maize_as_bran_area_base[ss,b,c]+napier_area_base[ss,b,c]+napier_hy_area_base[ss,b,c]\n",
    "                crop_primary_area_base[(ss,b,c)]=sunflower_as_meal_area_base[ss,b,c]+maize_as_bran_area_base[ss,b,c]\n",
    "\n",
    "                # actual area for this simulation (ha)\n",
    "                pasture_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'pasture']/(pasture_annual_yield)\n",
    "                grass_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'grass']/(grass_annual_yield)\n",
    "                stover_area[ss,b,c]=(1/1000)*(feed_intake[ss,b,c,'stover'])/(stover_annual_yield)                \n",
    "                sunflower_as_meal_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'sunflower_seed_meal']/(sunflower_annual_yield*scenario_parameters[(ss,'yield_scale_factor_sunflower',s)])\n",
    "                maize_as_bran_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'maize_bran']/(maize_annual_yield*scenario_parameters[(ss,'yield_scale_factor_maize',s)])\n",
    "                napier_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'napier']/(napier_annual_yield)            \n",
    "                napier_hy_area[ss,b,c]=(1/1000)*feed_intake[ss,b,c,'napier_hy']/(napier_hy_annual_yield) \n",
    "                \n",
    "                # summarize area for total feed, primary feeds, and primary crops\n",
    "                feed_area[(ss,b,c)]=grass_area[ss,b,c]+stover_area[ss,b,c]+sunflower_as_meal_area[ss,b,c]+maize_as_bran_area[ss,b,c]+napier_area[ss,b,c]+pasture_area[ss,b,c]\n",
    "                feed_primary_area[(ss,b,c)]=grass_area[ss,b,c]+sunflower_as_meal_area[ss,b,c]+maize_as_bran_area[ss,b,c]+napier_area[ss,b,c]+pasture_area[ss,b,c]\n",
    "                crop_primary_area[(ss,b,c)]=sunflower_as_meal_area[ss,b,c]+maize_as_bran_area[ss,b,c]\n",
    "                total_grasslands_area[(ss,b,c)]=grass_area[ss,b,c]+napier_area[ss,b,c]+pasture_area[ss,b,c]+napier_hy_area[ss,b,c]\n",
    "                \n",
    "                \n",
    "                # Define feed available for feeding based on yield and area (defined as 'production') (kg DM/animal/d for every month)\n",
    "                pasture_production=[]## production available for feeding (kg DM/head/day)\n",
    "                napier_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                napier_hy_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                napier_silage_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                grass_production=[]  ## production available for feeding (kg DM/head/day)\n",
    "                stover_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                maize_bran_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                sunflower_seed_meal_production=[] ## production available for feeding (kg DM/head/day)\n",
    "                \n",
    "                \n",
    "                for i in range(0,12):\n",
    "                \n",
    "                    napier_production.append(0)\n",
    "                    napier_hy_production.append(0)\n",
    "                    napier_silage_production.append(0)\n",
    "                    grass_production.append(0)\n",
    "                    stover_production.append(0)\n",
    "                    pasture_production.append(0)\n",
    "                \n",
    "                \n",
    "                # For maize bran and sunflower cake there is no seasonal variation and therefore \n",
    "                # monthly consumption is simply annual production divided by 12\n",
    "                \n",
    "                sunflower_yld=sunflower_annual_yield*1000/365\n",
    "                maize_bran_yld=maize_annual_yield*1000/365\n",
    "\n",
    "                maize_bran_production=maize_as_bran_area[ss,b,c]*maize_bran_yld\n",
    "                sunflower_seed_meal_production=sunflower_as_meal_area[ss,b,c]*sunflower_yld\n",
    "\n",
    "                # Define feed available in each month for stovers, pasture, and fodder crops, based on seasonality parameters\n",
    "                # as defined in SM of article\n",
    "                \n",
    "                # Specify months pertaining to different seasons \n",
    "                # The first 6 months of the year are set as the dry season, the last 6 are the rainy season\n",
    "                \n",
    "            \n",
    "                am=[0,1,2,3,4,5,6,7,8,9,10,11]  ## all months\n",
    "                ld=[0,1,2,3,4,5]  # long dry season\n",
    "                lr=[6,7,8,9,10,11] # long rainy season\n",
    "                sr1=[0,1,2,3,4,5]  # months in which stover is rationed (after harvest) (see SM to article)\n",
    "                \n",
    "                \n",
    "                 # Load relative yield data from excel; these parameters are described in the SM of article\n",
    "                wsgn = ryld_data.cell_value(1,1) \n",
    "                dsdn = ryld_data.cell_value(2,1) \n",
    "                wsg_p= ryld_data.cell_value(1,2) \n",
    "                dsd_p= ryld_data.cell_value(2,2) \n",
    "            \n",
    "                # define the fraction of fodder from wet season that is ensiled for dry season feeding (only for cows)   \n",
    "                ensilage_fodder_frac={}\n",
    "                \n",
    "                if (c == 'cow'):\n",
    "                    ensilage_fodder_frac[(ss,b)]=scenario_parameters[(ss,'ensilage_fodder_frac',s)]\n",
    "              \n",
    "                elif (c != 'cow'):\n",
    "                    ensilage_fodder_frac[(ss,b)]=0\n",
    "                   \n",
    "\n",
    "                # For every month of the year specify feed available based on land footprint and seasonality parameters\n",
    "                for i in lr:\n",
    "                    napier_production[i]=wsgn*(1/12)*(1-ensilage_fodder_frac[(ss,b)])*napier_area[ss,b,c]*1000*napier_annual_yield/(30)\n",
    "                    grass_production[i]=wsg_p*(1/30)*grass_area[ss,b,c]*1000*grass_annual_yield/12\n",
    "                    pasture_production[i]=wsg_p*(1/30)*pasture_area[ss,b,c]*1000*pasture_annual_yield/12\n",
    "                    napier_silage_production[i]=0\n",
    "                    stover_production[i]=(.2)*(1/6)*stover_area[ss,b,c]*1000*stover_annual_yield/30\n",
    "                \n",
    "                for i in ld:\n",
    "                    napier_production[i]=dsdn*(1/12)*napier_area[ss,b,c]*1000*napier_annual_yield/(30)\n",
    "                    napier_hy_production[i]=dsdn*(1/12)*napier_hy_area[ss,b,c]*1000*napier_hy_annual_yield/(30)\n",
    "                    grass_production[i]=(1-ensilage_fodder_frac[(ss,b)])*dsd_p*(1/30)*grass_area[ss,b,c]*1000*grass_annual_yield/12\n",
    "                    pasture_production[i]=(1-ensilage_fodder_frac[(ss,b)])*dsd_p*(1/30)*pasture_area[ss,b,c]*1000*pasture_annual_yield/12\n",
    "                    napier_silage_production[i]=(1/6)*(ensilage_fodder_frac[(ss,b)])*(wsgn*(1/2)*napier_area[ss,b,c]*1000*napier_annual_yield/(30))\n",
    "                    \n",
    "                # Stover is rationed equally over 6 months of the year (immediately following crop harvest)    \n",
    "                for i in sr1:\n",
    "                    stover_production[i]=(1/6)*stover_area[ss,b,c]*1000*stover_annual_yield/30\n",
    "                       \n",
    "                ## Specify actual diets for the current animal based on the above simulated seasonal feed availability\n",
    "                \n",
    "                if (c =='cow'):\n",
    "    \n",
    "                    # Specify feed intakes for cows\n",
    "                    # For concentrates feeds are specified based on production cycle\n",
    "                    # For non-concentrates, feeds are specified by month based on seasonality of these feeds as defined above\n",
    "                    \n",
    "                    feed_storage.add_external_roughage(grass_production,'Native grass', 'default')\n",
    "                    feed_storage.add_external_roughage(pasture_production ,'Pasture', 'default')\n",
    "                    feed_storage.add_external_roughage(stover_production, m_stover, 'default')\n",
    "                    feed_storage.add_external_roughage(napier_production, 'napier grass', 'default')\n",
    "                    feed_storage.add_external_roughage(napier_silage_production, 'napier silage', 'default')\n",
    "                    feed_storage.add_external_roughage(napier_hy_production, 'napier-legume', 'default')\n",
    "                \n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_early_lac[ss,b,c,'maize_bran'],'maize bran', 'early lactation')\n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_early_lac[ss,b,c,'sunflower_seed_meal'],'sunflower seed meal', 'early lactation')\n",
    "\n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_late_lac[ss,b,c,'maize_bran'],'maize bran', 'lactating')\n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_late_lac[ss,b,c,'sunflower_seed_meal'],'sunflower seed meal', 'lactating')\n",
    "\n",
    "                   \n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_other[ss,b,c,'maize_bran'],'maize bran', 'default')\n",
    "                    feed_storage.add_external_concentrate(feed_intake_cow_other[ss,b,c,'sunflower_seed_meal'],'sunflower seed meal', 'default')\n",
    "                \n",
    "                    \n",
    "                else: \n",
    "                    feed_storage.add_external_roughage(pasture_production,'Pasture', 'default')\n",
    "                    feed_storage.add_external_concentrate(maize_bran_production,'maize bran', 'default')\n",
    "                    feed_storage.add_external_concentrate(sunflower_seed_meal_production,'sunflower seed meal', 'default')\n",
    "                    feed_storage.add_external_roughage(napier_production, 'napier grass', 'default')\n",
    "                    feed_storage.add_external_roughage(napier_hy_production, 'napier-legume', 'default')\n",
    "                    feed_storage.add_external_roughage(napier_silage_production, 'napier silage', 'default')\n",
    "                    feed_storage.add_external_roughage(grass_production, 'Native grass', 'default')\n",
    "                    feed_storage.add_external_roughage(stover_production, 'maize stover', 'default')\n",
    "\n",
    "                d.feed_supply=feed_storage\n",
    "\n",
    "                #  RUN LIVSIM SIMULATIONS for a given breed and cohort #\n",
    "\n",
    "                output_df = d.run(max_simu_range)\n",
    "                \n",
    "                final_month = start_age + len(output_df)\n",
    "\n",
    "                #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
    "\n",
    "                '''\n",
    "                Results and summation/scaling \n",
    "                \n",
    "                Calculate the needed output parameters as averages from the entire simulation period.\n",
    "                This is performed by summing all output variables over simulation period, then dividing by time span\n",
    "                of simulation period to get annual averages\n",
    "                '''\n",
    "                \n",
    "                start_lact = 0 # has lactation started for cow , boolean\n",
    "                start_count = 0 # has the counting period for simulation started, boolean\n",
    "                months_lactating=0 # number of months cow has been lactating\n",
    "                simu_period=0 # months of simulation period\n",
    "                cmc=0 # calendar month count. Used for keeping track of the season/month when calculating outputs from livsim\n",
    "                \n",
    "                # initialize calendar month variables as zero\n",
    "                for i in range(0,12):\n",
    "                    count_calendar_month[i,ss,b,c]=0\n",
    "                \n",
    "                \n",
    "\n",
    "                # tallying variables\n",
    "                total_milk_yield = 0 ## kg/cow\n",
    "                faecal_Nitrogen_tally=0  # g N\n",
    "                urinary_Nitrogen_tally=0 # g N\n",
    "                dry_matter_intake_tally=0 # kg DM\n",
    "                body_weight_kg_tally=0 # body weight in kg\n",
    "                diet_Nitrogen_tally=0 # kg\n",
    "                frac_conc_tally=0 # fraction\n",
    "                frac_roughage_tally=0  # fraction\n",
    "                neutral_detergent_fibre_tally=0 # kg \n",
    "                acid_detergent_fibre_tally=0 # kg\n",
    "                metabolisable_energy_tally=0 # MJ\n",
    "                feed_intake_grass_tally=0 #kg\n",
    "                feed_intake_stover_tally=0 #kg\n",
    "                feed_intake_napier_tally=0 #kg\n",
    "                feed_intake_napier_hy_tally=0 #kg\n",
    "                feed_intake_pasture_tally=0 #kg\n",
    "                feed_intake_maize_bran_tally=0  #kg\n",
    "                feed_intake_sfsm_tally=0  #kg\n",
    "                crude_protein_tally=0  #kg\n",
    "\n",
    "                # summation variables\n",
    "                lifetime_milk = 0  ## kg/cow\n",
    "                average_body_weight_kg=0 # kg\n",
    "\n",
    "                total_diet_Nitrogen=0 # kg\n",
    "\n",
    "\n",
    "                for m in range(len(output_df)):\n",
    "\n",
    "                    # if cow start counting after lactation begins\n",
    "                    if ((((d.sex == 'f') and output_df.get_value(m,'Cow_milk_yield') > 0 ) and (start_lact == 0)) or ((d.sex == 'f') and (start_age == m))):\n",
    "\n",
    "                        if (c == 'cow'):\n",
    "                            if ((output_df.get_value(m,'Cow_age')-start_age)<12):\n",
    "                                cmc=(output_df.get_value(m,'Cow_age') - start_age)\n",
    "                                cmc=0\n",
    "                                start_lact = 1\n",
    "                            elif ((output_df.get_value(m,'Cow_age')-start_age)>12):\n",
    "                                cmc=int(12*(((output_df.get_value(m,'Cow_age') - start_age)/12)-1))\n",
    "                                cmc=0\n",
    "                                start_lact = 1\n",
    "                        else:\n",
    "                            cmc=0\n",
    "                        start_count = m\n",
    "\n",
    "                    # if male start counting immediately\n",
    "                    elif ((d.sex == 'm') and (start_count == 0)):\n",
    "                        start_count = m\n",
    "\n",
    "\n",
    "                    if (start_count > 0 ):\n",
    "                        simu_period += 1\n",
    "\n",
    "                        if (output_df.get_value(m,'Cow_dead') == 'True'):\n",
    "                            print('Animal ', b, ' ', c,' is dead')\n",
    "                            break\n",
    "\n",
    "\n",
    "                        # sum results variables over simulation period\n",
    "                        if (c == 'cow'):\n",
    "\n",
    "                            lifetime_milk += output_df.get_value(m,'Cow_milk_yield')\n",
    "\n",
    "                            if (output_df.get_value(m,'Cow_is_lactating')):\n",
    "\n",
    "                                months_lactating += 1\n",
    "\n",
    "                                if (output_df.get_value(m,'Cow_months_after_calving')<=5):\n",
    "                                    feed_intake_maize_bran =  feed_offer_cow[ss,b,c,'e_lac','maize_bran']*MAIZE_BRAN.__getattribute__('DM')*(1/1000)*30                \n",
    "\n",
    "\n",
    "                                    feed_intake_sfsm=feed_offer_cow[ss,b,c,'e_lac','sunflower_seed_meal']*SUNFLOWER_SEED_MEAL.__getattribute__('DM')*(1/1000)*30         \n",
    "\n",
    "                                elif (output_df.get_value(m,'Cow_months_after_calving')>5):\n",
    "                                    feed_intake_maize_bran=feed_offer_cow[ss,b,c,'l_lac','maize_bran']*MAIZE_BRAN.__getattribute__('DM')*(1/1000)*30                \n",
    "                                    feed_intake_sfsm=feed_offer_cow[ss,b,c,'l_lac','sunflower_seed_meal']*SUNFLOWER_SEED_MEAL.__getattribute__('DM')*(1/1000)*30         \n",
    "\n",
    "\n",
    "                                feed_intake_grass=grass_production[cmc]*GRASS.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_pasture=pasture_production[cmc]*PASTURE.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_stover=stover_production[cmc]*STOVER.__getattribute__('DM')*(1/1000)*30                   \n",
    "                                feed_intake_napier=napier_production[cmc]*NAPIER.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_napier_hy=napier_hy_production[cmc]*NAPIER_HY.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_napier_silage=napier_silage_production[cmc]*NAPIER_SILAGE.__getattribute__('DM')*(1/1000)*30                  \n",
    "\n",
    "\n",
    "                            else:\n",
    "\n",
    "                                feed_intake_sfsm=feed_offer_cow[ss,b,c,'other','sunflower_seed_meal']*SUNFLOWER_SEED_MEAL.__getattribute__('DM')*(1/1000)*30            \n",
    "                                feed_intake_maize_bran=feed_offer_cow[ss,b,c,'other','maize_bran']*MAIZE_BRAN.__getattribute__('DM')*(1/1000)*30           \n",
    "\n",
    "\n",
    "                                feed_intake_grass=grass_production[cmc]*GRASS.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_pasture=pasture_production[cmc]*PASTURE.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_stover=stover_production[cmc]*STOVER.__getattribute__('DM')*(1/1000)*30 \n",
    "                                feed_intake_napier=napier_production[cmc]*NAPIER.__getattribute__('DM')*(1/1000)*30 \n",
    "                                feed_intake_napier_hy=napier_hy_production[cmc]*NAPIER_HY.__getattribute__('DM')*(1/1000)*30                  \n",
    "                                feed_intake_napier_silage=napier_silage_production[cmc]*NAPIER_SILAGE.__getattribute__('DM')*(1/1000)*30                  \n",
    "\n",
    "\n",
    "                        else:\n",
    "\n",
    "                            feed_intake_maize_bran=maize_bran_production*MAIZE_BRAN.__getattribute__('DM')*(1/1000)*30               \n",
    "                            feed_intake_sfsm=sunflower_seed_meal_production*SUNFLOWER_SEED_MEAL.__getattribute__('DM')*(1/1000)*30             \n",
    "                            feed_intake_grass=grass_production[cmc]*GRASS.__getattribute__('DM')*(1/1000)*30                \n",
    "                            feed_intake_stover=stover_production[cmc]*STOVER.__getattribute__('DM')*(1/1000)*30                  \n",
    "                            feed_intake_napier=napier_production[cmc]*NAPIER.__getattribute__('DM')*(1/1000)*30 \n",
    "                            feed_intake_napier_hy=napier_hy_production[cmc]*NAPIER_HY.__getattribute__('DM')*(1/1000)*30\n",
    "                            feed_intake_pasture=pasture_production[cmc]*PASTURE.__getattribute__('DM')*(1/1000)*30                  \n",
    "\n",
    "\n",
    "                        total_rough_supply=(feed_intake_pasture+\n",
    "                                            feed_intake_napier+\n",
    "                                            feed_intake_napier_hy+\n",
    "                                            feed_intake_grass+\n",
    "                                            feed_intake_stover+\n",
    "                                            feed_intake_napier_silage)\n",
    "\n",
    "                        total_conc_supply=feed_intake_maize_bran+feed_intake_sfsm\n",
    "\n",
    "\n",
    "                        if (total_conc_supply > 0):\n",
    "                            fraction_concentrate_supply_eaten=output_df.get_value(m,'FI_concentrate_intake')/total_conc_supply\n",
    "                        else:\n",
    "                            print('Fraction conc eaten is zero')\n",
    "                            fraction_concentrate_supply_eaten = 1\n",
    "\n",
    "\n",
    "\n",
    "                        faecal_Nitrogen_tally +=output_df.get_value(m,'GrPr_faecal_N')\n",
    "                        urinary_Nitrogen_tally+=output_df.get_value(m,'GrPr_urinary_N')\n",
    "\n",
    "                        fraction_roughage_supply_eaten=output_df.get_value(m,'FI_roughage_intake')/total_rough_supply\n",
    "                        frac_conc_tally += fraction_concentrate_supply_eaten\n",
    "                        frac_roughage_tally += fraction_roughage_supply_eaten\n",
    "\n",
    "                        dry_matter_intake_tally += 30*output_df.get_value(m,'FI_feed_intake') \n",
    "                        body_weight_kg_tally += output_df.get_value(m,'Cow_body_weight') \n",
    "                        neutral_detergent_fibre_tally +=30*fraction_concentrate_supply_eaten*(sunflower_seed_meal_production*SUNFLOWER_SEED_MEAL.parameters['NDF']['value']+maize_bran_production*MAIZE_BRAN.parameters['NDF']['value'])+ 30*fraction_roughage_supply_eaten*(napier_silage_production[cmc]*NAPIER_SILAGE.parameters['NDF']['value']  +       pasture_production[cmc]*PASTURE.parameters['NDF']['value'][cmc]+grass_production[cmc]*GRASS.parameters['NDF']['value'][cmc]+stover_production[cmc]*STOVER.parameters['NDF']['value'][cmc]+napier_production[cmc]*NAPIER.parameters['NDF']['value'][cmc]+napier_hy_production[cmc]*NAPIER_HY.parameters['NDF']['value'][cmc])\n",
    "                        acid_detergent_fibre_tally += 30*fraction_concentrate_supply_eaten*(sunflower_seed_meal_production*SUNFLOWER_SEED_MEAL.parameters['ADF']['value']+maize_bran_production*MAIZE_BRAN.parameters['ADF']['value'])+30*fraction_roughage_supply_eaten*(napier_silage_production[cmc]*NAPIER_SILAGE.parameters['ADF']['value']       + pasture_production[cmc]*PASTURE.parameters['ADF']['value'][cmc]+grass_production[cmc]*GRASS.parameters['ADF']['value'][cmc]+stover_production[cmc]*STOVER.parameters['ADF']['value'][cmc]+napier_production[cmc]*NAPIER.parameters['ADF']['value'][cmc]+napier_hy_production[cmc]*NAPIER_HY.parameters['ADF']['value'][cmc]) \n",
    "                        crude_protein_tally += 30*fraction_concentrate_supply_eaten*(sunflower_seed_meal_production*SUNFLOWER_SEED_MEAL.parameters['CP']['value']+maize_bran_production*MAIZE_BRAN.parameters['CP']['value'])+30*fraction_roughage_supply_eaten*( napier_silage_production[cmc]*NAPIER_SILAGE.parameters['CP']['value'] +        pasture_production[cmc]*PASTURE.parameters['CP']['value'][cmc]+grass_production[cmc]*GRASS.parameters['CP']['value'][cmc]+stover_production[cmc]*STOVER.parameters['CP']['value'][cmc]+napier_production[cmc]*NAPIER.parameters['CP']['value'][cmc]+napier_hy_production[cmc]*NAPIER_HY.parameters['CP']['value'][cmc]) \n",
    "                        metabolisable_energy_tally += 30*(fraction_concentrate_supply_eaten*(sunflower_seed_meal_production*SUNFLOWER_SEED_MEAL.parameters['ME']['value']+maize_bran_production*MAIZE_BRAN.parameters['ME']['value'])+fraction_roughage_supply_eaten*( napier_silage_production[cmc]*NAPIER_SILAGE.parameters['ME']['value'][cmc]   +        pasture_production[cmc]*PASTURE.parameters['ME']['value'][cmc]+grass_production[cmc]*GRASS.parameters['ME']['value'][cmc]+stover_production[cmc]*STOVER.parameters['ME']['value'][cmc]+napier_production[cmc]*NAPIER.parameters['ME']['value'][cmc]+napier_hy_production[cmc]*NAPIER_HY.parameters['ME']['value'][cmc]))\n",
    "\n",
    "\n",
    "                        # Cow_MP_shortage\n",
    "                        mo_mp_short[ss,b,c,cmc] += output_df.get_value(m,'Cow_MP_shortage')\n",
    "                        mo_me_short[ss,b,c,cmc] += output_df.get_value(m,'Cow_ME_shortage')\n",
    "\n",
    "                        mo_fi_mb[ss,b,c,cmc] += fraction_concentrate_supply_eaten*maize_bran_production\n",
    "                        mo_fi_sc[ss,b,c,cmc] += fraction_concentrate_supply_eaten*sunflower_seed_meal_production\n",
    "                        mo_fi_sto[ss,b,c,cmc] += fraction_roughage_supply_eaten*stover_production[cmc]\n",
    "\n",
    "                        mo_fi_pt[ss,b,c,cmc] += fraction_roughage_supply_eaten*pasture_production[cmc]\n",
    "                        mo_fi_np[ss,b,c,cmc] += fraction_roughage_supply_eaten*napier_production[cmc]\n",
    "                        mo_fi_nhy[ss,b,c,cmc] += fraction_roughage_supply_eaten*stover_production[cmc]\n",
    "                        mo_fi_ns[ss,b,c,cmc] += fraction_roughage_supply_eaten*napier_silage_production[cmc]\n",
    "                        mo_fi_grs[ss,b,c,cmc] += fraction_roughage_supply_eaten*grass_production[cmc]\n",
    "\n",
    "\n",
    "                        count_calendar_month[cmc,ss,b,c] += 1\n",
    "\n",
    "\n",
    "                        feed_intake_grass_tally += fraction_roughage_supply_eaten*feed_intake_grass*GRASS.__getattribute__('DM')/1000              \n",
    "                        feed_intake_stover_tally += fraction_roughage_supply_eaten*feed_intake_stover*STOVER.__getattribute__('DM')/1000                 \n",
    "                        feed_intake_napier_tally += fraction_roughage_supply_eaten*(feed_intake_napier+feed_intake_napier_silage)*NAPIER.__getattribute__('DM')/1000                 \n",
    "                        feed_intake_napier_hy_tally += fraction_roughage_supply_eaten*(feed_intake_napier_hy)*NAPIER.__getattribute__('DM')/1000                 \n",
    "                        feed_intake_maize_bran_tally += fraction_concentrate_supply_eaten*feed_intake_maize_bran*MAIZE_BRAN.__getattribute__('DM')/1000         \n",
    "                        feed_intake_sfsm_tally += fraction_concentrate_supply_eaten*feed_intake_sfsm*SUNFLOWER_SEED_MEAL.__getattribute__('DM')/1000                 \n",
    "                        feed_intake_pasture_tally += fraction_roughage_supply_eaten*feed_intake_pasture*PASTURE.__getattribute__('DM')/1000\n",
    "\n",
    "                        if (m == (len(output_df)-1)):\n",
    "\n",
    "                            final_bodyweight_kg[ss,b,c]=output_df.get_value(m,'Cow_body_weight')\n",
    "                            final_age_months[ss,b,c]=output_df.get_value(m,'Cow_age')\n",
    "\n",
    "\n",
    "\n",
    "                    if (cmc < 11):\n",
    "                        cmc += 1\n",
    "                    else: \n",
    "                        cmc=0\n",
    "\n",
    "\n",
    "                # Average dietary and excretion parameters over the entire simulation period of the animal\n",
    "\n",
    "                faecal_Nitrogen[ss,b,c] = 12*faecal_Nitrogen_tally/simu_period\n",
    "                urinary_Nitrogen[ss,b,c] = 12*urinary_Nitrogen_tally/simu_period\n",
    "                NDF[ss,b,c] = 12*neutral_detergent_fibre_tally/simu_period\n",
    "                ADF[ss,b,c] = 12*acid_detergent_fibre_tally/simu_period\n",
    "                ME[ss,b,c] = 12*metabolisable_energy_tally/simu_period\n",
    "                DMI[ss,b,c] = 12*dry_matter_intake_tally/simu_period\n",
    "                Body_weight[ss,b,c] = body_weight_kg_tally/simu_period\n",
    "                gross_energy[ss,b,c]=18.4*DMI[ss,b,c]\n",
    "                crude_protein_intake[ss,b,c]=(1/1000)*12*crude_protein_tally/simu_period\n",
    "\n",
    "                ## calculated variables\n",
    "                DMD[ss,b,c] = 83.58 - 0.824*100*(1/1000)*ADF[ss,b,c] /DMI[ss,b,c]  + 2.626*100*(1/1000)*(crude_protein_intake[ss,b,c])/DMI[ss,b,c]\n",
    "                VS[ss,b,c]=((18.4*DMI[ss,b,c] *DMD[ss,b,c]/100)+0.04*gross_energy[ss,b,c])*(0.92/18.45)\n",
    "                Ym[ss,b,c]= 3.1 - 0.243*DMI[ss,b,c]/365 + 0.0059*NDF[ss,b,c]/DMI[ss,b,c]  + 0.0057*10*DMD[ss,b,c]/100\n",
    "\n",
    "                pasture_intake[ss,b,c]=12*feed_intake_pasture_tally/simu_period\n",
    "                grass_intake[ss,b,c]=12*feed_intake_grass_tally/simu_period\n",
    "                stover_intake[ss,b,c]=12*feed_intake_stover_tally/simu_period\n",
    "                napier_intake[ss,b,c]=12*feed_intake_napier_tally/simu_period\n",
    "                napier_hy_intake[ss,b,c]=12*feed_intake_napier_hy_tally/simu_period\n",
    "                maize_bran_intake[ss,b,c]=12*feed_intake_maize_bran_tally/simu_period\n",
    "                sfsm_intake[ss,b,c]=12*feed_intake_sfsm_tally/simu_period\n",
    "\n",
    "                if (c == 'cow'):\n",
    "                    # convert milk to FPCM\n",
    "                    milk_yield_reg[ss,b]=12*lifetime_milk*(0.337+.116*(1/10)*d.parameters['milk_butterfat_content']['value']+0.06*(1/10)*d.parameters['milk_crude_protein_content']['value'])/simu_period\n",
    "                    milk_yield_lact[ss,b]=d.parameters['lactation_duration']['value']*lifetime_milk* (0.337+.116*(1/10)*d.parameters['milk_butterfat_content']['value']+0.06*(1/10)*d.parameters['milk_crude_protein_content']['value'])/months_lactating\n",
    "                   \n",
    "                    RESULTS['v3_0_simu_period_months_'+str(ss)+'_'+str(b)]=simu_period\n",
    "                    RESULTS['v3_0_lactation_months_'+str(ss)+'_'+str(b)]=months_lactating\n",
    "                    \n",
    "                fraction_roughage_eaten[ss,b,c]=frac_roughage_tally/simu_period\n",
    "                fraction_conc_eaten[ss,b,c]=frac_conc_tally/simu_period\n",
    "\n",
    "                faecal_Nitrogen[ss,b,c]=total_faecal_N \n",
    "                urinary_Nitrogen[ss,b,c]=total_urinary_N               \n",
    "                Ym[ss,b,c]=Ym_per_cow_per_year\n",
    "                VS[ss,b,c]=volatile_solids_per_cow_per_year\n",
    "                DMI[ss,b,c]=total_DMI\n",
    "                DMD[ss,b,c]=total_diet_digestibility\n",
    "                ME[ss,b,c]=total_ME /(simu_period)\n",
    "                NDF[ss,b,c]=12*neutral_detergent_fibre_tally/(simu_period)\n",
    "                ADF[ss,b,c]=12*acid_detergent_fibre_tally/(simu_period)\n",
    "                crude_protein_intake[ss,b,c]=total_crude_protein/1000\n",
    "\n",
    "\n",
    "                # agricultural soils N2O based on IPCC framework\n",
    "                N_MMT_Frac={} # Fraction of Nitrogen excretion that is managed vs. excreted on pasture\n",
    "                \n",
    "                # for EFs specific to a breed, specify current iteration's EF\n",
    "                if (b=='improved'):\n",
    "                    MCF=MCF_IMP\n",
    "                    FRAC_MANURE_PRP=FRAC_MANURE_PRP_IMP\n",
    "                    FRAC_MANURE_MMT=FRAC_MANURE_MMT_IMP\n",
    "                    EF_4=EF_4_IMP\n",
    "                    EF_5=EF_5_IMP\n",
    "                    FRAC_vol=FRAC_vol_IMP\n",
    "                    FRAC_leach=FRAC_leach_IMP\n",
    "                     \n",
    "                elif (b=='local'):\n",
    "                    MCF=MCF_LOC\n",
    "                    FRAC_MANURE_PRP=FRAC_MANURE_PRP_LOC\n",
    "                    FRAC_MANURE_MMT=FRAC_MANURE_MMT_LOC\n",
    "                    EF_4=EF_4_LOC\n",
    "                    EF_5=EF_5_LOC\n",
    "                    FRAC_vol=FRAC_vol_LOC\n",
    "                    FRAC_leach=FRAC_leach_LOC\n",
    "                    \n",
    "                \n",
    "                N_manure_to_Pasture=(1/1000)*(faecal_Nitrogen[ss,b,c]+urinary_Nitrogen[ss,b,c])*FRAC_MANURE_PRP\n",
    "                N_manure_to_Storage=(1/1000)*(faecal_Nitrogen[ss,b,c]*(1-EF_3_FAECAL)+urinary_Nitrogen[ss,b,c]*(1-EF_3_URINARY))*FRAC_MANURE_MMT\n",
    "                \n",
    "                N_manure_for_Pasture=N_manure_to_Pasture\n",
    "                N_manure_for_Soil=N_manure_to_Storage*(1-FRAC_loss_MMS)\n",
    "                \n",
    "                # Allocate total n from manure management thats applied to soils to different land categories\n",
    "                # This is a management parameter; how N from management is allocated to different land uses\n",
    "                N_MMT_Frac['grass']=N_manure_for_Soil*grass_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                N_MMT_Frac['napier']=N_manure_for_Soil*napier_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                N_MMT_Frac['napier_hy']=N_manure_for_Soil*napier_hy_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                N_MMT_Frac['maize_bran']=N_manure_for_Soil*maize_as_bran_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                N_MMT_Frac['sunflower_seed_meal']=N_manure_for_Soil*sunflower_as_meal_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                N_MMT_Frac['stover']=N_manure_for_Soil*stover_area[ss,b,c]/feed_area[(ss,b,c)]\n",
    "                \n",
    "                # Fertilizer, manure, and cr residue inputs into land categories (kg n/ha/yr)\n",
    "                Fsn={}\n",
    "                Fon={}\n",
    "                Fcr={}\n",
    "                Fprp={}\n",
    "                \n",
    "                '''\n",
    "                This section accounts for all N inputs into different land categories based on IPCC (2006) chapter on Soils\n",
    "                This includes synthetic fertilizer (sn)\n",
    "                Organic fertilizer (fon) (manure)\n",
    "                Crop residues (fcr)\n",
    "                Manure deposited from grazing cattle\n",
    "                '''\n",
    "                \n",
    "                # Synthetic fertilizer n (kg n/ha/yr) (sources for these values are desribed in article)\n",
    "                Fsn['grass']=0\n",
    "                Fsn['pasture']=0\n",
    "                Fsn['napier']=0\n",
    "                Fsn['napier_hy']=0\n",
    "                Fsn['maize_bran']=20+scenario_parameters[(ss,'nitrous_oxide_factor_maize',s)]\n",
    "                Fsn['sunflower_seed_meal']=20+scenario_parameters[(ss,'nitrous_oxide_factor_sunflower',s)]\n",
    "                Fsn['stover']=10\n",
    "                \n",
    "                # Organic fertilizer n (kg n/ha/yr)\n",
    "                Fon['grass']=N_MMT_Frac['grass']/grass_area[ss,b,c]\n",
    "                Fon['pasture']=0\n",
    "                Fon['napier']=(N_MMT_Frac['napier']+N_MMT_Frac['maize_bran'])/napier_area[ss,b,c]\n",
    "                Fon['napier_hy']=(N_MMT_Frac['napier_hy'])/napier_hy_area[ss,b,c] \n",
    "                Fon['maize_bran']=0 \n",
    "                Fon['sunflower_seed_meal']=0\n",
    "                Fon['stover']=(N_MMT_Frac['stover']+N_MMT_Frac['sunflower_seed_meal'])/stover_area[ss,b,c]\n",
    "                \n",
    "                # Residue n (kg n/ha/yr) (these are calculated outside of model based on IPCC 2006 equations)\n",
    "                Fcr['grass']=2.2\n",
    "                Fcr['pasture']=4.5\n",
    "                Fcr['napier']=16.9\n",
    "                Fcr['napier_hy']=16.9\n",
    "                Fcr['maize_bran']=15.2\n",
    "                Fcr['sunflower_seed_meal']=10.6\n",
    "                Fcr['stover']=15.2\n",
    "                \n",
    "                # Tally total synthetic N applied to get emissions from manufacturing below\n",
    "                Synth_N_tot_kg_yr=(stover_area[ss,b,c]*Fsn['stover']+\n",
    "                                   sunflower_as_meal_area[ss,b,c]*Fsn['sunflower_seed_meal']+\n",
    "                                   Fsn['maize_bran']*maize_as_bran_area[ss,b,c])\n",
    "                \n",
    "                if (grass_area[ss,b,c]>0):\n",
    "                    \n",
    "                    Fprp['grass']=N_manure_for_Pasture/grass_area[ss,b,c]\n",
    "                    \n",
    "                elif (grass_area[ss,b,c]==0):\n",
    "                    \n",
    "                    Fprp['grass']=0\n",
    "    \n",
    "                # Direct N2O from all sources per land category (kg n2o-n/yr)\n",
    "                soil_N2O_feed_direct['pasture']=pasture_area[ss,b,c]*(EF_1*(Fsn['pasture']+Fcr['pasture']+Fon['pasture']))\n",
    "                soil_N2O_feed_direct['grass']=grass_area[ss,b,c]*(EF_1*(Fsn['grass']+Fcr['grass']+Fon['grass'])+EF_3_PRP*Fprp['grass'])\n",
    "                soil_N2O_feed_direct['napier']=napier_area[ss,b,c]*EF_1*(Fsn['napier']+ Fon['napier']+Fcr['napier'])\n",
    "                soil_N2O_feed_direct['napier_hy']=napier_area[ss,b,c]*EF_1*(Fsn['napier_hy']+ Fon['napier_hy']+Fcr['napier_hy'])\n",
    "                soil_N2O_feed_direct['stover']=stover_area[ss,b,c]*allocation_factor_stover*EF_1*(Fsn['stover']+ Fon['stover']+Fcr['stover'])\n",
    "                soil_N2O_feed_direct['maize_bran']=maize_as_bran_area[ss,b,c]*allocation_factor_maize_bran*EF_1*(Fsn['maize_bran']+Fcr['maize_bran'])\n",
    "                soil_N2O_feed_direct['sunflower_seed_meal']=sunflower_as_meal_area[ss,b,c]*allocation_factor_sunflower_seed_meal*EF_1*(Fsn['sunflower_seed_meal']+Fcr['sunflower_seed_meal'])\n",
    "                \n",
    "                # indirect (volatilization) (kg n2o-n/yr)\n",
    "                soil_N2O_feed_adn['pasture']=pasture_area[ss,b,c]*EF_4_adn*(Fsn['pasture']*FRAC_gasf+(Fon['pasture'])*FRAC_gasm)#\n",
    "                soil_N2O_feed_adn['grass']=grass_area[ss,b,c]*EF_4_adn*(Fsn['grass']*FRAC_gasf+(Fon['grass']+Fprp['grass'])*FRAC_gasm)\n",
    "                soil_N2O_feed_adn['napier']=napier_area[ss,b,c]*EF_4_adn*(Fsn['napier']*FRAC_gasf+Fon['napier']*FRAC_gasm)\n",
    "                soil_N2O_feed_adn['napier_hy']=napier_area[ss,b,c]*EF_4_adn*(Fsn['napier_hy']*FRAC_gasf+Fon['napier_hy']*FRAC_gasm)\n",
    "                soil_N2O_feed_adn['stover']=stover_area[ss,b,c]*allocation_factor_stover*EF_4_adn*(Fsn['stover']*FRAC_gasf+Fon['stover']*FRAC_gasm)\n",
    "                soil_N2O_feed_adn['maize_bran']=maize_as_bran_area[ss,b,c]*allocation_factor_maize_bran*EF_4_adn*(Fsn['maize_bran']*FRAC_gasf+Fon['maize_bran']*FRAC_gasm)\n",
    "                soil_N2O_feed_adn['sunflower_seed_meal']=sunflower_as_meal_area[ss,b,c]*allocation_factor_sunflower_seed_meal*EF_4_adn*(Fsn['sunflower_seed_meal']*FRAC_gasf+Fon['sunflower_seed_meal']*FRAC_gasm)\n",
    "              \n",
    "                # indirect (leaching) (kg n2o-n/yr)\n",
    "                soil_N2O_feed_lch['pasture']=pasture_area[ss,b,c]*EF_5_lch*(Fsn['pasture']+Fcr['pasture']+Fon['pasture'])\n",
    "                soil_N2O_feed_lch['grass']=grass_area[ss,b,c]*EF_5_lch*(Fsn['grass']+Fcr['grass']+Fon['grass']+Fprp['grass'])\n",
    "                soil_N2O_feed_lch['napier']=napier_area[ss,b,c]*EF_5_lch*(Fsn['napier']+ Fon['napier']+Fcr['napier'])\n",
    "                soil_N2O_feed_lch['napier_hy']=napier_area[ss,b,c]*EF_5_lch*(Fsn['napier_hy']+ Fon['napier_hy']+Fcr['napier_hy'])\n",
    "                soil_N2O_feed_lch['stover']=stover_area[ss,b,c]*allocation_factor_stover*EF_5_lch*(Fsn['stover']+ Fon['stover']+Fcr['stover'])\n",
    "                soil_N2O_feed_lch['maize_bran']=maize_as_bran_area[ss,b,c]*allocation_factor_maize_bran*EF_5_lch*(Fsn['maize_bran']+Fcr['maize_bran'])\n",
    "                soil_N2O_feed_lch['sunflower_seed_meal']=sunflower_as_meal_area[ss,b,c]*allocation_factor_sunflower_seed_meal*EF_5_lch*(Fsn['sunflower_seed_meal']+Fcr['sunflower_seed_meal'])\n",
    "               \n",
    "                # total (sum of direct and indirect)\n",
    "                soil_N2O_feed['pasture']=soil_N2O_feed_direct['pasture']+soil_N2O_feed_adn['pasture']+soil_N2O_feed_lch['pasture']\n",
    "                soil_N2O_feed['grass']=soil_N2O_feed_direct['grass']+soil_N2O_feed_adn['grass']+soil_N2O_feed_lch['grass']\n",
    "                soil_N2O_feed['napier']=soil_N2O_feed_direct['napier']+soil_N2O_feed_adn['napier']+soil_N2O_feed_lch['napier']\n",
    "                soil_N2O_feed['napier_hy']=soil_N2O_feed_direct['napier_hy']+soil_N2O_feed_adn['napier_hy']+soil_N2O_feed_lch['napier_hy']\n",
    "                soil_N2O_feed['stover']=soil_N2O_feed_direct['stover']+soil_N2O_feed_adn['stover']+soil_N2O_feed_lch['stover']\n",
    "                soil_N2O_feed['maize_bran']=soil_N2O_feed_direct['maize_bran']+soil_N2O_feed_adn['maize_bran']+soil_N2O_feed_lch['maize_bran']\n",
    "                soil_N2O_feed['sunflower_seed_meal']=soil_N2O_feed_direct['sunflower_seed_meal']+soil_N2O_feed_adn['sunflower_seed_meal']+soil_N2O_feed_lch['sunflower_seed_meal']\n",
    "                \n",
    "                \n",
    "                # calculate GHG emissions from all sources\n",
    "                gross_energy[ss,b,c]=18.4*DMI[ss,b,c]\n",
    "                \n",
    "                enteric_CH4[ss,b,c]=28*gross_energy[ss,b,c]*Ym[ss,b,c]/(100*55.65)\n",
    "                \n",
    "                manure_CH4[ss,b,c]=(28*VS[ss,b,c]*B0*MCF*0.67)\n",
    "                manure_N2O_direct[ss,b,c]=265*(44/28)*FRAC_MANURE_MMT*((1/1000)*faecal_Nitrogen[ss,b,c]*EF_3_FAECAL+(1/1000)*urinary_Nitrogen[ss,b,c]*EF_3_URINARY)\n",
    "                manure_N2O_vol[ss,b,c]=265*(44/28)*((1/1000)*FRAC_MANURE_MMT*(faecal_Nitrogen[ss,b,c]+urinary_Nitrogen[ss,b,c]))*FRAC_vol*EF_4\n",
    "                manure_N2O_leach[ss,b,c]=265*(44/28)*((1/1000)*FRAC_MANURE_MMT*(faecal_Nitrogen[ss,b,c]+urinary_Nitrogen[ss,b,c]))*FRAC_leach*EF_5\n",
    "                manure_N2O[ss,b,c]=manure_N2O_direct[ss,b,c]+manure_N2O_vol[ss,b,c]+manure_N2O_leach[ss,b,c]\n",
    "               \n",
    "                soil_N2O_direct_grassland[ss,b,c]=265*(44/28)*(soil_N2O_feed_direct['grass']+soil_N2O_feed_direct['napier_hy']+soil_N2O_feed_direct['napier']+soil_N2O_feed_direct['pasture'])\n",
    "                soil_N2O_indirect_leach_grassland[ss,b,c]=265*(44/28)*(soil_N2O_feed_lch['grass']+soil_N2O_feed_lch['napier_hy']+soil_N2O_feed_lch['napier']+soil_N2O_feed_lch['pasture'])\n",
    "                soil_N2O_indirect_vol_grassland[ss,b,c]=265*(44/28)*(soil_N2O_feed_adn['grass']+soil_N2O_feed_adn['napier_hy']+soil_N2O_feed_adn['napier']+soil_N2O_feed_adn['pasture'])\n",
    "                \n",
    "                soil_N2O_direct_cropland[ss,b,c]=265*(44/28)*(soil_N2O_feed_direct['maize_bran']+soil_N2O_feed_direct['sunflower_seed_meal']+soil_N2O_feed_direct['stover'])\n",
    "                soil_N2O_indirect_leach_cropland[ss,b,c]=265*(44/28)*(soil_N2O_feed_lch['maize_bran']+soil_N2O_feed_lch['sunflower_seed_meal']+soil_N2O_feed_lch['stover'])\n",
    "                soil_N2O_indirect_vol_cropland[ss,b,c]=265*(44/28)*(soil_N2O_feed_adn['maize_bran']+soil_N2O_feed_adn['sunflower_seed_meal']+soil_N2O_feed_adn['stover'])\n",
    "                \n",
    "                soil_N2O_indirect_cropland[ss,b,c]=soil_N2O_indirect_vol_cropland[ss,b,c]+soil_N2O_indirect_leach_cropland[ss,b,c]\n",
    "                soil_N2O_cropland[ss,b,c]=soil_N2O_direct_cropland[ss,b,c]+soil_N2O_indirect_cropland[ss,b,c]\n",
    "                \n",
    "                soil_N2O_indirect_grassland[ss,b,c]=soil_N2O_indirect_vol_grassland[ss,b,c]+soil_N2O_indirect_leach_grassland[ss,b,c]\n",
    "                soil_N2O_grassland[ss,b,c]=soil_N2O_direct_grassland[ss,b,c]+soil_N2O_indirect_grassland[ss,b,c]\n",
    "                \n",
    "                energy_use_CO2[ss,b,c]=feed_embodied_CO2['maize_bran']*maize_bran_production*365/1000 + feed_embodied_CO2['sunflower_seed_meal']*sunflower_seed_meal_production*365/1000\n",
    "                energy_use_CO2[ss,b,c]= energy_use_CO2[ss,b,c]+Synth_N_tot_kg_yr*fertilizer_embodied_CO2\n",
    "            \n",
    "                # Allocate to milk\n",
    "                # Allocation factor is 0.9\n",
    "                alloc_factor=0.9\n",
    "                enteric_CH4[ss,b,c]=alloc_factor*enteric_CH4[ss,b,c]\n",
    "                manure_CH4[ss,b,c]=alloc_factor*manure_CH4[ss,b,c]\n",
    "                manure_N2O_direct[ss,b,c]=alloc_factor*manure_N2O_direct[ss,b,c]\n",
    "                manure_N2O_vol[ss,b,c]=alloc_factor*manure_N2O_vol[ss,b,c]\n",
    "                manure_N2O_leach[ss,b,c]=alloc_factor*manure_N2O_leach[ss,b,c]\n",
    "                manure_N2O[ss,b,c]=alloc_factor*manure_N2O[ss,b,c]\n",
    "                soil_N2O_grassland[ss,b,c]=alloc_factor*soil_N2O_grassland[ss,b,c]\n",
    "                soil_N2O_cropland[ss,b,c]=alloc_factor*soil_N2O_cropland[ss,b,c]\n",
    "                energy_use_CO2[ss,b,c]=alloc_factor*energy_use_CO2[ss,b,c]\n",
    "\n",
    "                ''' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
    "                Monte Carlo simulations \n",
    "                \n",
    "                ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n",
    "                '''\n",
    "                \n",
    "                \n",
    "                uncertainty_sheet= data_wb.sheet_by_index(2)\n",
    "                rvars=0\n",
    "                random_var={}\n",
    "                num_vars=uncertainty_sheet.nrows\n",
    "\n",
    "                for r in range(0,MC_sims):\n",
    "                    for i in range(0,num_vars+1):\n",
    "                        random_var[r,i]=np.random.normal(100,uncertainty_sheet.cell_value(rvars+1, 1))/100\n",
    "\n",
    "  \n",
    "                    Ym_rd[ss,b,c]=Ym[ss,b,c]*random_var[r,0]\n",
    "                    B0_rd[ss,b,c]=B0*random_var[r,1]\n",
    "                    MCF_rd[ss,b,c]=MCF*random_var[r,2]\n",
    "                    EF_3_FAECAL_rd[ss,b,c]=EF_3_FAECAL*random_var[r,3]\n",
    "                    EF_3_URINARY_rd[ss,b,c]=EF_3_URINARY*random_var[r,4]\n",
    "                    EF_4_rd=EF_4*random_var[r,5]\n",
    "                    EF_5_rd=EF_5*random_var[r,6]\n",
    "                    FRAC_vol_rd= FRAC_vol*random_var[r,7]\n",
    "                    FRAC_leach_rd=FRAC_leach*random_var[r,8]\n",
    "                    EF_1_rd=EF_1*random_var[r,9]\n",
    "                    EF_3_PRP_rd=EF_3_PRP*random_var[r,10]\n",
    "                    EF_4_adn_rd=EF_4_adn*random_var[r,11]\n",
    "                    EF_5_lch_rd=EF_5_lch*random_var[r,12]\n",
    "                    \n",
    "\n",
    "                    # Pasture\n",
    "                    if (s == 'Base'):\n",
    "                        grass_area_rd=grass_area[ss,b,c]*random_var[r,14]\n",
    "                        pasture_area_rd=pasture_area[ss,b,c]*random_var[r,14]\n",
    "                        napier_area_rd=napier_area[ss,b,c]*random_var[r,15]\n",
    "                        napier_hy_area_rd=napier_area[ss,b,c]*random_var[r,15]\n",
    "                        maize_as_bran_area_rd=maize_as_bran_area[ss,b,c]*random_var[r,16]\n",
    "                        sunflower_as_meal_area_rd=sunflower_as_meal_area[ss,b,c]*random_var[r,17]\n",
    "                        stover_area_rd=stover_area[ss,b,c]*random_var[r,18]\n",
    "                        herd_size_rd[(ss,b,'herd_pop')]=herd_size[(ss,b,'herd_pop')]*random_var[r,19]\n",
    "                        \n",
    "                    else:\n",
    "                        grass_area_rd=grass_area[ss,b,c]\n",
    "                        pasture_area_rd=pasture_area[ss,b,c]\n",
    "                        napier_area_rd=napier_area[ss,b,c]\n",
    "                        napier_hy_area_rd=napier_area[ss,b,c]\n",
    "                        maize_as_bran_area_rd=maize_as_bran_area[ss,b,c]\n",
    "                        sunflower_as_meal_area_rd=sunflower_as_meal_area[ss,b,c]\n",
    "                        stover_area_rd=stover_area[ss,b,c]\n",
    "                        herd_size_rd[(ss,b,'herd_pop')]=herd_size[(ss,b,'herd_pop')]\n",
    "                        \n",
    "                    crop_primary_area_rd=crop_primary_area[ss,b,c]*random_var[r,23]\n",
    "                    feed_primary_area_rd=feed_primary_area[ss,b,c]*random_var[r,24]\n",
    "                        \n",
    "                    \n",
    "                # direct n2o from all sources per land category (kg n2o-n/yr)\n",
    "                    soil_N2O_feed_direct_rd['pasture']=pasture_area_rd*(EF_1_rd*(Fsn['pasture']+Fcr['pasture']))\n",
    "                    soil_N2O_feed_direct_rd['grass']=grass_area_rd*(EF_1_rd*(Fsn['grass']+Fcr['grass']+Fon['grass'])+EF_3_PRP_rd*Fprp['grass'])\n",
    "                    soil_N2O_feed_direct_rd['napier']=napier_area_rd*EF_1_rd*(Fsn['napier']+ Fon['napier']+Fcr['napier'])\n",
    "                    soil_N2O_feed_direct_rd['napier_hy']=napier_hy_area_rd*EF_1_rd*(Fsn['napier_hy']+ Fon['napier_hy']+Fcr['napier_hy'])\n",
    "                    soil_N2O_feed_direct_rd['stover']=stover_area_rd*allocation_factor_stover*EF_1_rd*(Fsn['stover']+ Fon['stover']+Fcr['stover'])\n",
    "                    soil_N2O_feed_direct_rd['maize_bran']=maize_as_bran_area_rd*allocation_factor_maize_bran*EF_1_rd*(Fsn['maize_bran']+Fcr['maize_bran'])\n",
    "                    soil_N2O_feed_direct_rd['sunflower_seed_meal']=sunflower_as_meal_area_rd*allocation_factor_sunflower_seed_meal*EF_1_rd*(Fsn['sunflower_seed_meal']+Fcr['sunflower_seed_meal'])\n",
    "                \n",
    "                # indirect (volatilization) (kg n2o-n/yr)\n",
    "                    soil_N2O_feed_adn_rd['pasture']=pasture_area_rd*EF_4_adn_rd*(Fsn['pasture']*FRAC_gasf+(Fon['pasture'])*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['grass']=grass_area_rd*EF_4_adn_rd*(Fsn['grass']*FRAC_gasf+(Fon['grass']+Fprp['grass'])*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['napier']=napier_area_rd*EF_4_adn_rd*(Fsn['napier']*FRAC_gasf+Fon['napier']*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['napier_hy']=napier_hy_area_rd*EF_4_adn_rd*(Fsn['napier_hy']*FRAC_gasf+Fon['napier_hy']*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['stover']=stover_area_rd*allocation_factor_stover*EF_4_adn_rd*(Fsn['stover']*FRAC_gasf+Fon['stover']*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['maize_bran']=maize_as_bran_area_rd*allocation_factor_maize_bran*EF_4_adn_rd*(Fsn['maize_bran']*FRAC_gasf+Fon['maize_bran']*FRAC_gasm)\n",
    "                    soil_N2O_feed_adn_rd['sunflower_seed_meal']=sunflower_as_meal_area_rd*allocation_factor_sunflower_seed_meal*EF_4_adn_rd*(Fsn['sunflower_seed_meal']*FRAC_gasf+Fon['sunflower_seed_meal']*FRAC_gasm)\n",
    "              \n",
    "                # indirect (leaching) (kg n2o-n/yr)\n",
    "                    soil_N2O_feed_lch_rd['pasture']=pasture_area_rd*EF_5_lch_rd*(Fsn['pasture']+Fcr['pasture']+Fon['pasture'])\n",
    "                    soil_N2O_feed_lch_rd['grass']=grass_area_rd*EF_5_lch_rd*(Fsn['grass']+Fcr['grass']+Fon['grass']+Fprp['grass'])\n",
    "                    soil_N2O_feed_lch_rd['napier']=napier_area_rd*EF_5_lch_rd*(Fsn['napier']+ Fon['napier']+Fcr['napier'])\n",
    "                    soil_N2O_feed_lch_rd['napier_hy']=napier_hy_area_rd*EF_5_lch_rd*(Fsn['napier_hy']+ Fon['napier_hy']+Fcr['napier_hy'])\n",
    "                    soil_N2O_feed_lch_rd['stover']=stover_area_rd*allocation_factor_stover*EF_5_lch_rd*(Fsn['stover']+ Fon['stover']+Fcr['stover'])\n",
    "                    soil_N2O_feed_lch_rd['maize_bran']=maize_as_bran_area_rd*allocation_factor_maize_bran*EF_5_lch_rd*(Fsn['maize_bran']+Fcr['maize_bran'])\n",
    "                    soil_N2O_feed_lch_rd['sunflower_seed_meal']=sunflower_as_meal_area_rd*allocation_factor_sunflower_seed_meal*EF_5_lch_rd*(Fsn['sunflower_seed_meal']+Fcr['sunflower_seed_meal'])\n",
    "               \n",
    "\n",
    "                    soil_N2O_feed_grass_rd=(soil_N2O_feed_direct_rd['grass']+soil_N2O_feed_lch_rd['grass']+soil_N2O_feed_adn_rd['grass'])\n",
    "                    soil_N2O_feed_napier_rd=(soil_N2O_feed_direct_rd['napier']+soil_N2O_feed_lch_rd['napier']+soil_N2O_feed_adn_rd['napier'])\n",
    "                    soil_N2O_feed_napier_hy_rd=(soil_N2O_feed_direct_rd['napier_hy']+soil_N2O_feed_lch_rd['napier_hy']+soil_N2O_feed_adn_rd['napier_hy'])\n",
    "                    soil_N2O_feed_maize_bran_rd=(soil_N2O_feed_direct_rd['maize_bran']+soil_N2O_feed_lch_rd['maize_bran']+soil_N2O_feed_adn_rd['maize_bran'])\n",
    "                    soil_N2O_feed_sunflower_seed_meal_rd=(soil_N2O_feed_direct_rd['sunflower_seed_meal']+soil_N2O_feed_lch_rd['sunflower_seed_meal']+soil_N2O_feed_adn_rd['sunflower_seed_meal'])\n",
    "                    soil_N2O_feed_stover_rd=(soil_N2O_feed_direct_rd['stover']+soil_N2O_feed_lch_rd['stover']+soil_N2O_feed_adn_rd['stover'])\n",
    "                    soil_N2O_feed_pasture_rd=(soil_N2O_feed_direct_rd['pasture']+soil_N2O_feed_lch_rd['pasture']+soil_N2O_feed_adn_rd['pasture'])\n",
    "\n",
    "                    enteric_CH4_rd[r,ss,b,c]=alloc_factor*28*gross_energy[ss,b,c]*Ym_rd[ss,b,c]*(1/(100*55.65))\n",
    "                    manure_CH4_rd[r,ss,b,c]=alloc_factor*(28*VS[ss,b,c]*B0_rd[ss,b,c]*MCF_rd[ss,b,c]*0.67)\n",
    "                    manure_N2O_direct_rd[r,ss,b,c]=alloc_factor*265*(44/28)*FRAC_MANURE_MMT*((1/1000)*faecal_Nitrogen[ss,b,c]*EF_3_FAECAL_rd[ss,b,c]+(1/1000)*urinary_Nitrogen[ss,b,c]*EF_3_URINARY_rd[ss,b,c])\n",
    "                    manure_N2O_leach_rd[r,ss,b,c]=alloc_factor*265*(44/28)*FRAC_MANURE_MMT*((1/1000)*(faecal_Nitrogen[ss,b,c]+urinary_Nitrogen[ss,b,c]))*FRAC_leach_rd*EF_5_rd\n",
    "                    manure_N2O_vol_rd[r,ss,b,c]=alloc_factor*265*(44/28)*((1/1000)*FRAC_MANURE_MMT*(faecal_Nitrogen[ss,b,c]+urinary_Nitrogen[ss,b,c]))*FRAC_vol_rd*EF_4_rd\n",
    "                    manure_N2O_rd[r,ss,b,c]=alloc_factor*manure_N2O_direct_rd[r,ss,b,c]+manure_N2O_leach_rd[r,ss,b,c]+manure_N2O_vol_rd[r,ss,b,c]\n",
    "                    soil_N2O_rd[r,ss,b,c]=alloc_factor*265*(44/28)*(soil_N2O_feed_grass_rd*grass_area_rd+soil_N2O_feed_napier_rd*napier_area_rd+soil_N2O_feed_maize_bran_rd*maize_as_bran_area_rd+soil_N2O_feed_sunflower_seed_meal_rd*sunflower_as_meal_area_rd +soil_N2O_feed_stover_rd*stover_area_rd+soil_N2O_feed_pasture_rd*pasture_area_rd)\n",
    "                    energy_use_CO2_rd[r,ss,b,c]=alloc_factor*energy_use_CO2[ss,b,c]*random_var[r,25]\n",
    "\n",
    "                    grassland_area_expansion_emissions_rd_ss[r,ss]+=alloc_factor*random_var[r,20]*(44/12)*(herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*feed_primary_area_rd-herd_size_y_prev[(ss,b)]*herd_size[(ss,b,c)]*feed_primary_area_rd)*(luc_emission_coefficient_total_area)*(timestep/amort_period)/timestep \n",
    "                    cropland_area_expansion_emissions_rd_ss[r,ss]+=alloc_factor*random_var[r,21]*(44/12)*(herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*(crop_primary_area_rd)-herd_size_y_prev[(ss,b)]*herd_size[(ss,b,c)]*(crop_primary_area_rd))*(luc_emission_coefficient_grass_to_crop)*(timestep/amort_period)/timestep \n",
    "\n",
    "\n",
    "    '''\n",
    "    Aggregation and scaling                  \n",
    "                    \n",
    "                    \n",
    "    This section sums numeric results of subsectors, breeds, and cohorts to report final results for a given \n",
    "    production system. \n",
    "                \n",
    "    '''\n",
    "\n",
    "        \n",
    "    for ss in subsectors:\n",
    "        for b in breed:\n",
    "            if (herd_size[(ss,b,'herd_pop')] ==0):\n",
    "                \n",
    "                continue # If the current subsector has no animals of a given breed, skip to next iteration of loop\n",
    "                \n",
    "            '''\n",
    "            Results are grouped numerically from v1 to v9, with each grouping for a specific type of result variable\n",
    "            This is necessary in for the pandas dataframe storing results to keep the variables in order\n",
    "            \n",
    "            v1 -- main production/productivity results\n",
    "            v2 -- nutrient intakes and N excretion\n",
    "            v3 -- feed intakes, simulation variables \n",
    "            v4 -- animal specific emissions\n",
    "            v5 -- feed emission properties\n",
    "            v6 -- aggregated emissions \n",
    "            v7 -- aggregated emissions intensity\n",
    "            v8 -- dairy herd populations\n",
    "            v9 -- land footprint\n",
    "            '''\n",
    "            \n",
    "            \n",
    "            for c in cohort:\n",
    "                if (c == 'cow'):\n",
    "                    \n",
    "                    RESULTS['v1_1_Milk_yield_ss'+str(ss)+'_'+str(b)]=milk_yield_reg[ss,b]\n",
    "                    RESULTS['v1_1_Milk_yield_lactation_kg_per_lactation__ss'+str(ss)+'_'+str(b)]=milk_yield_lact[ss,b]\n",
    "                    RESULTS['v1_2_Milk_yield_total']+= herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,'cow')]*milk_yield_reg[ss,b]\n",
    "                    RESULTS['v1_2_Milk_yield_base_year_kg_yr']+=herd_size_y_prev[(ss,b)]*herd_size[(ss,b,'cow')]*milk_yield_reg[ss,b]\n",
    "                \n",
    "                # results per breed and cohort and subsector\n",
    "                RESULTS['v2_1_Dry_matter_intake_kg_per_head_per_yr_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=DMI[ss,b,c]\n",
    "                RESULTS['v2_2_Dry_matter_intake_kg_per_head_per_day_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=DMI[ss,b,c]/365\n",
    "                RESULTS['v2_3_Fraction_Roughage_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=fraction_roughage_eaten[ss,b,c]\n",
    "                RESULTS['v2_3_Fraction_Concentrate_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=fraction_conc_eaten[ss,b,c]\n",
    "                RESULTS['v2_3_Metabolisable_Energy_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=ME[ss,b,c]\n",
    "                RESULTS['v2_4_Dry_matter_digestibility_%_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=DMD[ss,b,c]\n",
    "                RESULTS['v2_5_1_Gross_energy_MJ_per_head_per_yr_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=gross_energy[ss,b,c]\n",
    "                RESULTS['v2_5_2_Gross_energy_MJ_per_head_per_day_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=gross_energy[ss,b,c]/365\n",
    "                RESULTS['v2_6_faecal_nitrogen_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=faecal_Nitrogen[ss,b,c]\n",
    "                RESULTS['v2_6_urinary_nitrogen_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=urinary_Nitrogen[ss,b,c] \n",
    "                RESULTS['v2_7_Ym_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=Ym[ss,b,c]\n",
    "                RESULTS['v2_8_VS_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=VS[ss,b,c]\n",
    "                RESULTS['v2_9_Acid_detergent_fibre_g_per_head_per_year_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=ADF[ss,b,c]\n",
    "                RESULTS['v2_9_Acid_detergent_fibre_%_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*ADF[ss,b,c]*(1/1000)/DMI[ss,b,c]\n",
    "                RESULTS['v2_10_Neutral_detergent_fibre_g_per_head_per_year_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=NDF[ss,b,c]\n",
    "                RESULTS['v2_10_Neutral_detergent_fibre_%_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*NDF[ss,b,c]*(1/1000)/DMI[ss,b,c]\n",
    "                RESULTS['v2_11_Crude_protein_kg_per_head_per_year'+str(b)+'_'+str(c)+'_'+str(ss)+'']=crude_protein_intake[ss,b,c]\n",
    "                RESULTS['v2_11_Crude_protein_%_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*crude_protein_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                \n",
    "                # feed intake results (% dm for specific ss, b, c)\n",
    "                RESULTS['v3_1_1_Feed_intake_grass_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*grass_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                RESULTS['v3_1_2_Feed_intake_stover_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*stover_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                RESULTS['v3_1_3_Feed_intake_napier_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*napier_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                RESULTS['v3_1_4_Feed_intake_maize_bran_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*maize_bran_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                RESULTS['v3_1_5_Feed_intake_sfsm_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*sfsm_intake[ss,b,c]/DMI[ss,b,c]\n",
    "                RESULTS['v3_1_6_Feed_intake_pasture_%_dm'+str(b)+'_'+str(c)+'_'+str(ss)+'']=100*pasture_intake[ss,b,c]/DMI[ss,b,c]\n",
    "\n",
    "                # feed intake results (kg dm for entire herd)\n",
    "                RESULTS['v3_2_6_Feed_intake_grass_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*grass_intake[ss,b,c]\n",
    "                RESULTS['v3_2_7_Feed_intake_stover_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*stover_intake[ss,b,c]\n",
    "                RESULTS['v3_2_8_Feed_intake_napier_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*napier_intake[ss,b,c]\n",
    "                RESULTS['v3_2_8_Feed_intake_napier_hy_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*napier_hy_intake[ss,b,c]\n",
    "                RESULTS['v3_2_9_Feed_intake_maize_bran_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*maize_bran_intake[ss,b,c]\n",
    "                RESULTS['v3_2_10_Feed_intake_sfsm_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*sfsm_intake[ss,b,c]\n",
    "                RESULTS['v3_2_11_Feed_intake_pasture_kg_dm']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*pasture_intake[ss,b,c]\n",
    "\n",
    "                \n",
    "                # seasonality results\n",
    "                if (c == 'cow'):\n",
    "                \n",
    "                    \n",
    "                    for i in range(0,cmc+1):\n",
    "                        # Met energy intake in month i\n",
    "                        RESULTS['v3_2_12_M_'+str(i)+'_metenergy_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=mo_dp_me[ss,b,c,i]/(count_calendar_month[i,ss,b,c])\n",
    "                        \n",
    "                        # Crude protein intake in month i\n",
    "                        RESULTS['v3_2_12_M_'+str(i)+'_cp_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=mo_dp_cp[ss,b,c,i]/(count_calendar_month[i,ss,b,c])\n",
    "            \n",
    "                        # Shortage of MP in month i \n",
    "                        RESULTS['v3_2_12_M'+str(i)+'_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']= mo_mp_short[ss,b,c,i]/(count_calendar_month[i,ss,b,c])\n",
    "                   \n",
    "                         # Shortage of ME in month i \n",
    "                        RESULTS['v3_2_12_M'+str(i)+'_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']= mo_me_short[ss,b,c,i]/(count_calendar_month[i,ss,b,c])\n",
    "                   \n",
    "                        \n",
    "                   \n",
    "             # average protein or energy deficit during dry season\n",
    "                    RESULTS['v3_3_mean_dryse_prot_deficit']=(\n",
    "                        RESULTS['v3_2_12_M1_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']\n",
    "                        +RESULTS['v3_2_12_M2_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']\n",
    "                        +RESULTS['v3_2_12_M3_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']\n",
    "                        +RESULTS['v3_2_12_M4_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']\n",
    "                        +RESULTS['v3_2_12_M5_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']\n",
    "                        +RESULTS['v3_2_12_M6_MPS_g_d_'+str(b)+'_'+str(c)+'_'+str(ss)+''])*(1/6)\n",
    "    \n",
    "                    RESULTS['v3_3_mean_dryse_ene_deficit']=(\n",
    "                        RESULTS['v3_2_12_M1_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']+\n",
    "                        RESULTS['v3_2_12_M2_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']+\n",
    "                        RESULTS['v3_2_12_M3_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']+\n",
    "                        RESULTS['v3_2_12_M4_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']+\n",
    "                        RESULTS['v3_2_12_M5_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+'']+\n",
    "                        RESULTS['v3_2_12_M6_MES_Mj_d_'+str(b)+'_'+str(c)+'_'+str(ss)+''])*(1/6)\n",
    "                     \n",
    "                \n",
    "                \n",
    "                # emissions per animal per year\n",
    "                RESULTS['v4_1_Enteric_CH4_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+'']=enteric_CH4[ss,b,c]\n",
    "                RESULTS['v4_2_Manure_CH4_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+'']=manure_CH4[ss,b,c] \n",
    "                RESULTS['v4_3_Manure_N2O_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+'']=manure_N2O[ss,b,c] \n",
    "                RESULTS['v4_4_Soil_N2O_cropland_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+''] =soil_N2O_cropland[ss,b,c] \n",
    "                RESULTS['v4_4_Soil_N2O_grassland_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+''] =soil_N2O_grassland[ss,b,c] \n",
    "                RESULTS['v4_5_Feed_CO2_kg_CO2eq_per_head_per_year_'+str(ss)+'_'+str(b)+'_'+str(c)+''] =energy_use_CO2[ss,b,c]\n",
    "\n",
    "                # emissions from soils\n",
    "                RESULTS['v5_Grasslands_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*soil_N2O_feed['grass']/grass_area[(ss,b,c)]\n",
    "                RESULTS['v5_Napier_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*soil_N2O_feed['napier']/napier_area[(ss,b,c)]\n",
    "                RESULTS['v5_Maize_bran_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*(1/allocation_factor_maize_bran)*soil_N2O_feed['maize_bran']/maize_as_bran_area[(ss,b,c)]\n",
    "                RESULTS['v5_Sunflower_cake_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*(1/allocation_factor_sunflower_seed_meal)*soil_N2O_feed['sunflower_seed_meal']/sunflower_as_meal_area[(ss,b,c)]\n",
    "                RESULTS['v5_stover_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*(1/allocation_factor_stover)*soil_N2O_feed['stover']/stover_area[(ss,b,c)]\n",
    "                \n",
    "                \n",
    "                if (pasture_area[(ss,b,c)] != 0): # To avoid divide by zero error if pasture area is zero\n",
    "                    RESULTS['v5_Pasture_N2O_kg_per_ha_per_yr_ss_'+str(ss)+'']=(44/28)*soil_N2O_feed['pasture']/pasture_area[(ss,b,c)]\n",
    " \n",
    "                # total area of feed land\n",
    "                RESULTS['v5_6_Total_feed_primary_area_expansion']+=herd_size[(ss,b,c)]*(herd_size[(ss,b,'herd_pop')]*feed_primary_area[(ss,b,c)]-herd_size_y_prev[(ss,b)]*feed_primary_area[(ss,b,c)])\n",
    "\n",
    "                \n",
    "                #  emissions results aggregated over breed and cohort for each lps\n",
    "                RESULTS['v6_1_Enteric_CH4_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*enteric_CH4[ss,b,c]/1000 \n",
    "                RESULTS['v6_2_Manure_CH4_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*manure_CH4[ss,b,c]/1000 \n",
    "                RESULTS['v6_3_Manure_N2O_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*manure_N2O[ss,b,c]/1000 \n",
    "                RESULTS['v6_4_Soil_N2O_cropland_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*soil_N2O_cropland[ss,b,c]/1000 \n",
    "                RESULTS['v6_4_Soil_N2O_grassland_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*soil_N2O_grassland[ss,b,c]/1000 \n",
    "                RESULTS['v6_4_2_Feed_CO2_aggregate_Mg_CO2eq_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*energy_use_CO2[ss,b,c]/1000 \n",
    "\n",
    "                # luc as part of aLCA\n",
    "                RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']+=(44/12)*herd_size[(ss,b,c)]*(herd_size[(ss,b,'herd_pop')]*crop_primary_area[(ss,b,c)]-herd_size_y_prev[(ss,b)]*crop_primary_area[(ss,b,c)])*(luc_emission_coefficient_grass_to_crop)*(timestep/amort_period)/timestep\n",
    "\n",
    "                RESULTS['v6_5_Direct_emissions_aggregate_Mg_CO2eq_per_lps']=(\n",
    "                    RESULTS['v6_1_Enteric_CH4_aggregate_Mg_CO2eq_per_lps']+\n",
    "                    RESULTS['v6_2_Manure_CH4_aggregate_Mg_CO2eq_per_lps']+\n",
    "                    RESULTS['v6_3_Manure_N2O_aggregate_Mg_CO2eq_per_lps']+\n",
    "                    RESULTS['v6_4_Soil_N2O_grassland_aggregate_Mg_CO2eq_per_lps']+\n",
    "                    RESULTS['v6_4_Soil_N2O_cropland_aggregate_Mg_CO2eq_per_lps']+\n",
    "                    RESULTS['v6_4_2_Feed_CO2_aggregate_Mg_CO2eq_per_lps'])\n",
    "\n",
    "                # herd populations\n",
    "                RESULTS['v8_1_Herd_size_'+str(ss)+'_'+str(b)+'_'+str(c)]=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v8_2_Herd_size_previous_year_'+str(ss)+'_'+str(b)+'_'+str(c)]=herd_size_y_prev[(ss,b)]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v8_3_Herd_size_base_'+str(ss)+'_'+str(b)+'_'+str(c)]=herd_size_base[(l,ss,b,'herd_pop')]*herd_size_base[(l,ss,b,c)]\n",
    "                RESULTS['v8_4_Herd_size_total_1000_head_per_lps']+=herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]/1000                \n",
    "                RESULTS['v8_6_Herd_size_total_TLU_per_lps']+=herd_size[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*TLU_equiv[(b,c)]\n",
    "                RESULTS['v8_5_Herd_size_total_Million_TLU_per_lps']=RESULTS['v8_6_Herd_size_total_TLU_per_lps']/1000000\n",
    "                \n",
    "                # feed areas aggregated categories \n",
    "                RESULTS['v9_1_1_Grasslands_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=total_grasslands_area[(ss,b,c)]\n",
    "                RESULTS['v9_1_2_Feed_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=feed_area[(ss,b,c)]\n",
    "                RESULTS['v9_2_1_Grasslands_area_ha_per_lps'] += herd_size[(ss,b,c)]*(herd_size[(ss,b,'herd_pop')]*total_grasslands_area[(ss,b,c)])\n",
    "                RESULTS['v9_2_2_Feed_area_ha_per_lps'] +=  herd_size[(ss,b,c)]*(herd_size[(ss,b,'herd_pop')]*(feed_area[(ss,b,c)]))\n",
    "                RESULTS['v9_2_3_Crop_primary_area_ha_per_lps']+=herd_size[(ss,b,c)]*(herd_size[(ss,b,'herd_pop')]*(feed_area[(ss,b,c)]-total_grasslands_area[(ss,b,c)]))\n",
    "                   \n",
    "                # feed areas individual (ha per head)\n",
    "                RESULTS['v9_3_1_Grass_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=grass_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_2_Stover_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=stover_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_3_Napier_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=napier_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_4_Napier_hy_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=napier_hy_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_5_Sunflower_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=sunflower_as_meal_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_6_Maize_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=maize_as_bran_area[(ss,b,c)]\n",
    "                RESULTS['v9_3_7_Pasture_area_ha_per_head_'+str(b)+'_'+str(c)+'_'+str(ss)+'']=pasture_area[(ss,b,c)]\n",
    "    \n",
    "                # feed areas  (ha per lps)\n",
    "                RESULTS['v9_4_1_Grass_area_ha_per_lps']+=grass_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_2_Stover_area_ha_per_lps']+=stover_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_3_Napier_area_ha_per_lps']+=napier_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_4_Napier_hy_area_ha_per_lps']+=napier_hy_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_5_Sunflower_area_ha_per_lps']+=sunflower_as_meal_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_6_Maize_area_ha_per_lps']+=maize_as_bran_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_7_Pasture_area_ha_per_lps']+=pasture_area[(ss,b,c)]*herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]\n",
    "                RESULTS['v9_4_8_Total_area_ha_per_lps']=RESULTS['v9_4_1_Grass_area_ha_per_lps']+RESULTS['v9_4_2_Stover_area_ha_per_lps']+RESULTS['v9_4_3_Napier_area_ha_per_lps']+RESULTS['v9_4_4_Napier_hy_area_ha_per_lps']+RESULTS['v9_4_5_Sunflower_area_ha_per_lps']+RESULTS['v9_4_6_Maize_area_ha_per_lps']                \n",
    "                \n",
    "                \n",
    "                \n",
    "    '''\n",
    "    Calculate biophysical feasibility of scenario \n",
    "\n",
    "    '''\n",
    "\n",
    "    if (check_feasibility == 1):\n",
    "        \n",
    "        RESULTS['v6_0_total_native_ecosystems_converted(ha)']=0\n",
    "        RESULTS['v6_0_total_grasslands_now_utilized(ha)']=0\n",
    "        RESULTS['v6_0_Units_limited']=0\n",
    "    \n",
    "        # load data on spatial land use, lives populations\n",
    "        system_feasibility_data=(cwd+str('\\\\land_feasibility.xlsx'))\n",
    "        system_feasibility=xlrd.open_workbook(system_feasibility_data)\n",
    "    \n",
    "        if (l == 'MRT'):\n",
    "            feas_sheet = system_feasibility.sheet_by_index(0)\n",
    "        elif (l == 'MRH'):\n",
    "            feas_sheet = system_feasibility.sheet_by_index(1)\n",
    "\n",
    "        # load parameters from excel\n",
    "        density={} # lives density per pixel (h/pixel)\n",
    "        dgl_available={} #  grassland available as percentage of pixel\n",
    "        frac_pop={} # the fraction of livestock in a pixel relative to total in the lps\n",
    "        \n",
    "        u=feas_sheet.nrows # number of rows in excel sheet\n",
    "\n",
    "        total_popn=0 # total population of cattle in lps (head)\n",
    "        \n",
    "        for i in range(1,u): \n",
    "            total_popn+=feas_sheet.cell_value(i,2)\n",
    "\n",
    "        for i in range(1,u): \n",
    "            \n",
    "            dgl_available[i]=feas_sheet.cell_value(i,1)\n",
    "            frac_pop[i]=feas_sheet.cell_value(i,2)/total_popn\n",
    "            \n",
    "        \n",
    "        # total expansion of primary feed area for this scenario (ha/lps)\n",
    "        total_expansion=RESULTS['v5_6_Total_feed_primary_area_expansion']\n",
    "        area_growth={} # amount of new land needed for conversion per pixel (sq km/unit)\n",
    "        \n",
    "        for i in range(1,u):\n",
    "            \n",
    "            # amount of new land needed for primary feed expansion (weighted by cattle pop'n) (sq km)\n",
    "            # pixel is 10km x 10km = 100 sq km \n",
    "            area_growth[i]=frac_pop[i]*total_expansion*.01   # (sq km per pixel)\n",
    "            \n",
    "            # if amount of new land needed (ha/100 sq km) is less than available grassland to be utilized (ha) (feas_sheet.cell_value(i,1))\n",
    "            \n",
    "            if (area_growth[i] <= dgl_available[i]):\n",
    "                \n",
    "                RESULTS['v6_0_total_grasslands_now_utilized(ha)']+=area_growth[i]\n",
    "                \n",
    "            # if amount of new land needed is greater than available grassland to be utilized\n",
    "            elif (area_growth[i] > dgl_available[i]):\n",
    "                \n",
    "                RESULTS['v6_0_total_native_ecosystems_converted(ha)']+=(area_growth[i]-dgl_available[i])\n",
    "                RESULTS['v6_0_total_grasslands_now_utilized(ha)']+=dgl_available[i]\n",
    "                RESULTS['v6_0_Units_limited']+=1\n",
    "            \n",
    "        fraction_feed_primary_area_expansion_converting_native_ecosystems=(RESULTS['v6_0_total_native_ecosystems_converted(ha)']/(RESULTS['v6_0_total_grasslands_now_utilized(ha)']+RESULTS['v6_0_total_native_ecosystems_converted(ha)']))\n",
    "            \n",
    "        RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']=fraction_feed_primary_area_expansion_converting_native_ecosystems*(44/12)*RESULTS['v5_6_Total_feed_primary_area_expansion']*(luc_emission_coefficient_total_area)*(timestep/amort_period)/timestep\n",
    "                    \n",
    "\n",
    "    # direct non co2 emissions\n",
    "    RESULTS['v7_1_Enteric_CH4_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_1_Enteric_CH4_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total'] \n",
    "    RESULTS['v7_2_Manure_CH4_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_2_Manure_CH4_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_3_Manure_N2O_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_3_Manure_N2O_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_4_Soil_N2O_cropland_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_4_Soil_N2O_cropland_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_4_Soil_N2O_grassland_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_4_Soil_N2O_grassland_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_4_Feed_CO2_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_4_2_Feed_CO2_aggregate_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "   \n",
    "    RESULTS['v7_5_Direct emissions intensity_kg_CO2eq_per_kg_milk_per_lps']=(\n",
    "        RESULTS['v7_1_Enteric_CH4_intensity_kg_CO2eq_per_kg_milk_per_lps']+\n",
    "        RESULTS['v7_2_Manure_CH4_intensity_kg_CO2eq_per_kg_milk_per_lps']+\n",
    "        RESULTS['v7_3_Manure_N2O_intensity_kg_CO2eq_per_kg_milk_per_lps']+\n",
    "        RESULTS['v7_4_Soil_N2O_cropland_intensity_kg_CO2eq_per_kg_milk_per_lps']\n",
    "        +RESULTS['v7_4_Soil_N2O_grassland_intensity_kg_CO2eq_per_kg_milk_per_lps']+\n",
    "        RESULTS['v7_4_Feed_CO2_intensity_kg_CO2eq_per_kg_milk_per_lps'])\n",
    "\n",
    "    \n",
    "    # luc emissions\n",
    "    RESULTS['v7_6_grassland_expansion_emissions_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_8_cropland_expansion_emissions_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v7_9_Indirect_emissions_intensity_kg_CO2eq_per_kg_milk_per_lps']=1000*(RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']+RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps'])/RESULTS['v1_2_Milk_yield_total']\n",
    "    \n",
    "    RESULTS['v7_5_Total emissions intensity_kg_CO2eq_per_kg_milk_per_lps']=RESULTS['v7_5_Direct emissions intensity_kg_CO2eq_per_kg_milk_per_lps']+RESULTS['v7_9_Indirect_emissions_intensity_kg_CO2eq_per_kg_milk_per_lps']\n",
    "    \n",
    "    # total\n",
    "    RESULTS['v7_10_Emissions_total_Mg_CO2eq_per_lps']=RESULTS['v6_5_Direct_emissions_aggregate_Mg_CO2eq_per_lps']+RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']+RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']\n",
    "    \n",
    "    # Emissions per tlu\n",
    "    RESULTS['v7_11_Enteric_CH4_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_1_Enteric_CH4_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_12_Manure_CH4_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_2_Manure_CH4_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_13_Manure_N2O_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_3_Manure_N2O_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_14_Soil_N2O_cropland_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_4_Soil_N2O_cropland_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_14_Soil_N2O_grassland_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_4_Soil_N2O_grassland_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_14_Feed_CO2_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_4_2_Feed_CO2_aggregate_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "   \n",
    "\n",
    "    RESULTS['v7_16_grassland_expansion_emissions_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v7_18_cropland_expansion_emissions_intensity_kg_CO2eq_per_TLU_per_lps']=1000*RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "   \n",
    "    # multiply by 100 to convert to a percentage\n",
    "    RESULTS['v1_5_Feed_intake_Mt_per_LPS']=(1/(1000*1000000))*(RESULTS['v3_2_6_Feed_intake_grass_kg_dm']+RESULTS['v3_2_7_Feed_intake_stover_kg_dm']+RESULTS['v3_2_8_Feed_intake_napier_kg_dm']+RESULTS['v3_2_9_Feed_intake_maize_bran_kg_dm']+RESULTS['v3_2_10_Feed_intake_sfsm_kg_dm']+RESULTS['v3_2_11_Feed_intake_pasture_kg_dm'])\n",
    "    RESULTS['v3_2_1_Feed_intake_grass_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_6_Feed_intake_grass_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_2_Feed_intake_stover_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_7_Feed_intake_stover_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_3_Feed_intake_napier_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_8_Feed_intake_napier_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_3_Feed_intake_napier_hy_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_8_Feed_intake_napier_hy_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_4_Feed_intake_maize_bran_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_9_Feed_intake_maize_bran_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_5_Feed_intake_sfsm_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_10_Feed_intake_sfsm_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    RESULTS['v3_2_6_Feed_intake_pasture_(% dm)']=(100/(1000*1000000))*RESULTS['v3_2_11_Feed_intake_pasture_kg_dm']/RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "\n",
    "    RESULTS['v3_3_1_Feed_intake_grass_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_6_Feed_intake_grass_kg_dm']\n",
    "    RESULTS['v3_3_2_Feed_intake_stover_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_7_Feed_intake_stover_kg_dm']\n",
    "    RESULTS['v3_3_3_Feed_intake_napier_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_8_Feed_intake_napier_kg_dm']\n",
    "    RESULTS['v3_3_3_Feed_intake_napier_hy_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_8_Feed_intake_napier_hy_kg_dm']\n",
    "    RESULTS['v3_3_4_Feed_intake_maize_bran_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_9_Feed_intake_maize_bran_kg_dm']\n",
    "    RESULTS['v3_3_5_Feed_intake_sfsm_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_10_Feed_intake_sfsm_kg_dm']\n",
    "    RESULTS['v3_3_6_Feed_intake_pasture_kg_per_TLU']=(1/RESULTS['v8_6_Herd_size_total_TLU_per_lps'])*RESULTS['v3_2_11_Feed_intake_pasture_kg_dm']\n",
    "\n",
    "    # summation\n",
    "    RESULTS['v9_15_Grass_land_use_ha_per_TLU']=RESULTS['v9_4_1_Grass_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_16_Stover_land_use_ha_per_TLU']=RESULTS['v9_4_2_Stover_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_17_Napier_land_use_ha_per_TLU']=RESULTS['v9_4_3_Napier_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_17_Napier_hy_land_use_ha_per_TLU']=RESULTS['v9_4_4_Napier_hy_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_18_Sunflower_land_use_ha_per_TLU']=RESULTS['v9_4_5_Sunflower_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_19_Maize_land_use_ha_per_TLU']=RESULTS['v9_4_6_Maize_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    RESULTS['v9_20_Pasture_land_use_ha_per_TLU']=RESULTS['v9_4_7_Pasture_area_ha_per_lps']/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "    \n",
    "    # Calculate aggregated productivity, land use variables\n",
    "    RESULTS['v1_3_Milk_productivity_kg_per_ha']= RESULTS['v1_2_Milk_yield_total']/(RESULTS['v9_4_8_Total_area_ha_per_lps'])\n",
    "    RESULTS['v1_4_Land_footprint_Mha_per_LPS']=  (1/1000000)*(RESULTS['v9_4_8_Total_area_ha_per_lps'])\n",
    "    RESULTS['v1_6_Feed_efficiency_kg_FPCM_per_Mg_feed_dm']=(1)*(1000/(1000*1000000))*RESULTS['v1_2_Milk_yield_total']/ RESULTS['v1_5_Feed_intake_Mt_per_LPS']\n",
    "    \n",
    "    # Total emissions \n",
    "    RESULTS['v_11_1_Emissions_total_Mg_CO2eq_per_lps']=RESULTS['v6_5_Direct_emissions_aggregate_Mg_CO2eq_per_lps']+(RESULTS['v6_8_cropland_expansion_emissions_Mg_CO2eq_per_lps']+RESULTS['v6_6_grassland_expansion_emissions_Mg_CO2eq_per_lps'])\n",
    "\n",
    "    # sort results alphabetically/alphanumerically\n",
    "    sorted_list=sorted(RESULTS.keys(), key=lambda x:x.lower())\n",
    "\n",
    "    results_ordered={}\n",
    "\n",
    "    for i in sorted_list:\n",
    "        results_ordered[i]=RESULTS[i]\n",
    "    \n",
    "    RESULTS={}\n",
    "    RESULTS=results_ordered\n",
    "    \n",
    "    \n",
    "    # Calculate uncertainty\n",
    "    for r in range(0,MC_sims):\n",
    "        direct_emission_absolute_rd[r]=0\n",
    "        \n",
    "        grassland_area_expansion_emissions_rd[r]=0\n",
    "        cropland_area_expansion_emissions_rd[r]=0\n",
    "                    \n",
    "                    \n",
    "        for ss in subsectors:\n",
    "             \n",
    "            for b in breed:\n",
    "                if (herd_size[(ss,b,'herd_pop')] ==0):\n",
    "                    continue\n",
    "                for c in cohort:\n",
    "                    direct_emission_absolute_rd[r]+= herd_size[(ss,b,'herd_pop')]*herd_size[(ss,b,c)]*(enteric_CH4_rd[r,ss,b,c]+manure_CH4_rd[r,ss,b,c]+manure_N2O_rd[r,ss,b,c]+soil_N2O_rd[r,ss,b,c]+energy_use_CO2_rd[r,ss,b,c])/1000\n",
    "                    \n",
    "        for ss in subsectors:    \n",
    "            grassland_area_expansion_emissions_rd[r]+=grassland_area_expansion_emissions_rd_ss[r,ss]*fraction_feed_primary_area_expansion_converting_native_ecosystems\n",
    "            cropland_area_expansion_emissions_rd[r]+=cropland_area_expansion_emissions_rd_ss[r,ss]\n",
    "                    \n",
    "                    \n",
    "                    \n",
    "        direct_emission_absolute_rd[r]+=(cropland_area_expansion_emissions_rd[r]+grassland_area_expansion_emissions_rd[r])\n",
    "\n",
    "        \n",
    "    variation_ab_squared=list(range(0, MC_sims))\n",
    "    total_ab_var=list(range(0, MC_sims))\n",
    "    total_ab=0\n",
    "    location_start_error=len(RESULTS)-MC_sims\n",
    "    total_ab_variation=0\n",
    "\n",
    "        \n",
    "    for r in range(0,MC_sims):\n",
    "        variation_ab_squared[r]=(direct_emission_absolute_rd[r]-RESULTS['v_11_1_Emissions_total_Mg_CO2eq_per_lps'])**2\n",
    "        total_ab_variation+=variation_ab_squared[r]\n",
    "            \n",
    "        \n",
    "    total_ab=math.sqrt(total_ab_variation/MC_sims)\n",
    "    ninetyfivepci=1.96*total_ab\n",
    "    \n",
    "    # Final uncertainty results\n",
    "    RESULTS['v_11_2_error_ab']=ninetyfivepci/1000000\n",
    "    RESULTS['v_11_2_error_int']=(1000)*ninetyfivepci/RESULTS['v1_2_Milk_yield_total']\n",
    "    RESULTS['v_11_2_error_per_tlu']=(1000)*ninetyfivepci/RESULTS['v8_6_Herd_size_total_TLU_per_lps']\n",
    "\n",
    "        \n",
    "        \n",
    "    return RESULTS\n",
    "                            "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}