diff --git a/config/config.default.yaml b/config/config.default.yaml
index 1d28223f7..be2e07d88 100644
--- a/config/config.default.yaml
+++ b/config/config.default.yaml
@@ -393,7 +393,12 @@ sector:
   bev_charge_rate: 0.011
   bev_avail_max: 0.95
   bev_avail_mean: 0.8
+  bev_charge_avail: 5
   v2g: true
+  endogenous_transport: false
+  EV_consumption_1car: 0.01 
+  ICE_consumption_1car: 0.033 
+  H2_consumption_1car: 0.02 
   land_transport_fuel_cell_share:
     2020: 0
     2025: 0
@@ -418,7 +423,8 @@ sector:
     2040: 0.3
     2045: 0.15
     2050: 0
-  transport_fuel_cell_efficiency: 0.5
+  transport_electric_vehicle_efficiency: 0.99
+  transport_fuel_cell_efficiency: 0.71 #0.5
   transport_internal_combustion_efficiency: 0.3
   agriculture_machinery_electric_share: 0
   agriculture_machinery_oil_share: 1
@@ -1058,4 +1064,5 @@ plotting:
     DC: "#8a1caf"
     DC-DC: "#8a1caf"
     DC link: "#8a1caf"
+    land transport demand: '#2596be'
     load: "#dd2e23"
diff --git a/doc/configtables/sector.csv b/doc/configtables/sector.csv
index 338cf34e6..7e78cb920 100644
--- a/doc/configtables/sector.csv
+++ b/doc/configtables/sector.csv
@@ -24,6 +24,10 @@ bev_charge_rate,MWh,float,The power consumption for one electric vehicle (EV) in
 bev_avail_max,--,float,The maximum share plugged-in availability for passenger electric vehicles.
 bev_avail_mean,--,float,The average share plugged-in availability for passenger electric vehicles.
 v2g,--,"{true, false}",Allows feed-in to grid from EV battery
+endogenous transport,--,"{true, false}",Allows to optimize land transport shares endogenously
+EV_consumption_1car,MWh_elec/hour,float,assuming 0.2 kWh/km https://github.com/PyPSA/pypsa-eur/blob/1fbe971ab8dab60d972d3a7b905b9cec7171c0ad/config/config.default.yaml#L382 and velocity of 50km/h
+ICE_consumption_1car,Mwh_oil/hour,float,"assuming 0.66 kWh_oil/km and 50 km/h (with the link efficiency 0.3, they are equivalent to 0.01 MWh_elec/hour"
+H2_consumption_1car,Mwh_H2/hour,float,"assuming 0.4 kWh_H2/km and 50km/h) (with efficiency 0.5, they are equivalent to 0.01 MWh_elec/hour"
 land_transport_fuel_cell _share,--,Dictionary with planning horizons as keys.,The share of vehicles that uses fuel cells in a given year
 land_transport_electric _share,--,Dictionary with planning horizons as keys.,The share of vehicles that uses electric vehicles (EV) in a given year
 land_transport_ice _share,--,Dictionary with planning horizons as keys.,The share of vehicles that uses internal combustion engines (ICE) in a given year. What is not EV or FCEV is oil-fuelled ICE.
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 52b2ddc27..c776b8f7b 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -212,6 +212,8 @@ PyPSA-Eur 0.9.0 (5th January 2024)
   different directory using the new ``root_dir`` argument
   (https://github.com/PyPSA/pypsa-eur/pull/771).
 
+* Add option to optimize fuel type shares of land transport endogenously in configuration file (https://github.com/PyPSA/pypsa-eur/pull/734)
+
 * Rule ``purge`` now initiates a dialog to confirm if purge is desired
   (https://github.com/PyPSA/pypsa-eur/pull/745).
 
diff --git a/rules/build_electricity.smk b/rules/build_electricity.smk
index 6b0926386..62819cf77 100644
--- a/rules/build_electricity.smk
+++ b/rules/build_electricity.smk
@@ -306,7 +306,7 @@ rule build_renewable_profiles:
         BENCHMARKS + "build_renewable_profiles_{technology}"
     threads: ATLITE_NPROCESSES
     resources:
-        mem_mb=ATLITE_NPROCESSES * 5000,
+        mem_mb=ATLITE_NPROCESSES * 9000,
     wildcard_constraints:
         technology="(?!hydro).*",  # Any technology other than hydro
     conda:
@@ -460,7 +460,7 @@ rule simplify_network:
         BENCHMARKS + "simplify_network/elec_s{simpl}"
     threads: 1
     resources:
-        mem_mb=12000,
+        mem_mb=24000,
     conda:
         "../envs/environment.yaml"
     script:
diff --git a/rules/build_sector.smk b/rules/build_sector.smk
index 4744aa250..a2bfc9874 100644
--- a/rules/build_sector.smk
+++ b/rules/build_sector.smk
@@ -424,7 +424,7 @@ rule build_industrial_production_per_country:
         + "industrial_production_per_country.csv",
     threads: 8
     resources:
-        mem_mb=1000,
+        mem_mb=5000,
     log:
         LOGS + "build_industrial_production_per_country.log",
     benchmark:
diff --git a/rules/retrieve.smk b/rules/retrieve.smk
index 5e1e3e598..d26d81341 100644
--- a/rules/retrieve.smk
+++ b/rules/retrieve.smk
@@ -378,4 +378,4 @@ if config["enable"]["retrieve"]:
         conda:
             "../envs/environment.yaml"
         script:
-            "../scripts/retrieve_monthly_fuel_prices.py"
+            "../scripts/retrieve_monthly_fuel_prices.py"
\ No newline at end of file
diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk
index 7ca8857db..b20a8877b 100644
--- a/rules/solve_myopic.smk
+++ b/rules/solve_myopic.smk
@@ -23,6 +23,8 @@ rule add_existing_baseyear:
         existing_solar="data/existing_infrastructure/solar_capacity_IRENA.csv",
         existing_onwind="data/existing_infrastructure/onwind_capacity_IRENA.csv",
         existing_offwind="data/existing_infrastructure/offwind_capacity_IRENA.csv",
+        existing_transport=RESOURCES + "transport_data_s{simpl}_{clusters}.csv",
+        temp_air_total=RESOURCES + "temp_air_total_elec_s{simpl}_{clusters}.nc",
     output:
         RESULTS
         + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
@@ -108,6 +110,7 @@ rule solve_sector_network_myopic:
     resources:
         mem_mb=config["solving"]["mem"],
         walltime=config["solving"].get("walltime", "12:00:00"),
+        disk_mb=180000 #500000
     benchmark:
         (
             BENCHMARKS
diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py
index cb1f51c8b..ab9b4d6a2 100644
--- a/scripts/add_brownfield.py
+++ b/scripts/add_brownfield.py
@@ -67,6 +67,16 @@ def add_brownfield(n, n_p, year):
                 & (c.df[f"{attr}_nom_opt"] < threshold)
             ],
         )
+        check_transport = ( (snakemake.config["sector"]["land_transport_electric_share"][year] is not None ) 
+        or (snakemake.config["sector"]["land_transport_fuel_cell_share"][year] is not None)
+        or (snakemake.config["sector"]["land_transport_ice_share"][year] is not None))
+        
+        if not snakemake.config["sector"]["endogenous_transport"] or check_transport:
+            n_p.mremove(
+                c.name,
+                c.df.index[c.df.carrier.str.contains("land transport" or "V2G" or "EV battery storage" or "BEV charger")
+                ],
+            )
 
         # copy over assets but fix their capacity
         c.df[f"{attr}_nom"] = c.df[f"{attr}_nom_opt"]
@@ -119,6 +129,92 @@ def add_brownfield(n, n_p, year):
             n.links.loc[new_pipes, "p_nom"] = 0.0
             n.links.loc[new_pipes, "p_nom_min"] = 0.0
 
+def adjust_EVs(n, n_p, year):
+    # set p_min_pu and p_max_pu for solved network, so that only the EV link for the current time horizon
+    # is not constraint (constraining all land transport link with p_min_pu/p_max_pu while p_nom_extentable=True leads to infeasible by gurobi
+
+    for car in ["EV","fuel cell","oil"]:
+        cartype = "land transport " + car
+        if not n.links[(n.links.carrier==cartype) ].index.empty:
+            lifetime_EV = n.links.lifetime[n.links[(n.links.carrier==cartype) ].index[0]]
+            i = 0
+            while (year-lifetime_EV+i) < year:
+                if not n.links_t.p_min_pu[(n.links.filter(like=cartype +"-"+str(int(year-lifetime_EV+i)),axis=0)).index].empty:
+                    #p_set =  n_p.loads_t.p[n.loads[n.loads.carrier.str.contains('land transport demand')].index]
+                    #eff = n.links_t.efficiency[(n.links.filter(like="land transport EV-"+str(int(year-lifetime_EV+i)),axis=0)).index]
+                    #p_set = p_set.add_suffix(' EV-'+str(int(year-lifetime_EV+i)))
+                    #p_set = p_set.drop([col for col in p_set.columns if col in p_set.columns and col not in eff.columns], axis=1)
+                    #pnom = (p_set.divide(eff)).max()
+                    pu=n.links_t.p_min_pu[(n.links.filter(like=cartype +"-"+str(int(year-lifetime_EV+i)),axis=0)).index] #p_set.divide(eff)/pnom
+                    n.links_t.p_max_pu[(n.links.filter(like=cartype +"-"+str(int(year-lifetime_EV+i)),axis=0)).index] = pu.values
+                    #n.links_t.p_max_pu[(n.links.filter(like="land transpo-"+str(int(year-lifetime_EV+i)),axis=0)).index] = pu.values
+                i = i+1   
+    
+    #Get existing capacity of EV batteries
+    ev_battery_storage = n.stores[n.stores.carrier.str.contains("Li ion|EV battery storage")]
+    sum_ev_battery = ev_battery_storage.groupby("bus")["e_nom"].sum()
+    print("EV battery storage",ev_battery_storage.groupby("bus")["e_nom"].sum())
+    sum_ev_battery.index = sum_ev_battery.index.str.replace(' EV battery','',regex=False)
+
+    #Set e_nom_max to remaining capacity in relation to maximum capacity of fully electric car fleet of 2015
+    ev_store = n.stores.carrier.str.contains('Li ion|EV battery storage') 
+    ev_store_ext = n.stores[ev_store].query("e_nom_extendable").index
+
+    #Caclulate existing EV land transport capacity 
+    ev_land_transport = n.links[n.links.carrier.str.contains("land transport EV")]
+    sum_ev_transport = ev_land_transport.groupby("bus0")["p_nom"].sum()
+    sum_ev_transport.index = sum_ev_transport.index.str.replace(' EV battery', '', regex=False)
+    print("sum ev transport",sum_ev_transport)
+
+    bev_availability = snakemake.config['sector']["bev_availability"]
+    for o in opts:
+
+        if "bevavail" not in o:
+            continue
+        oo = o.split("+")
+        bev_availability = float(oo[1])
+    ev_transport_to_store = sum_ev_transport/snakemake.config['sector']['EV_consumption_1car']*snakemake.config['sector']['bev_energy']*snakemake.config['sector']['bev_charge_avail']*bev_availability
+    print(sum_ev_battery)
+    
+    if sum(sum_ev_battery):
+        enommax = n.stores.loc[ev_store_ext,'e_nom_max']
+        enommax.index = enommax.index.str.replace(' EV battery storage-'+str(year), '', regex=False) 
+        print(enommax)
+        #subtract existing battery capacity but add capacity if more EVs are used in land transport than is reflected in EV battery capacity
+        add_enommax = ev_transport_to_store-sum_ev_battery
+        add_enommax[add_enommax<0] = 0
+        newenommax = enommax-sum_ev_battery+add_enommax
+        newenommax[newenommax < 0] = 0
+        n.stores.loc[ev_store_ext,'e_nom_max'] = newenommax.values
+
+    BEV_charger = n.links[n.links.carrier.str.contains("BEV charger")]
+    sum_bev_charger = BEV_charger.groupby("bus0")["p_nom"].sum()
+    print("sum bev charger",sum_bev_charger)
+   
+    sum_bev_charger.index = sum_bev_charger.index.str.replace(' low voltage', '', regex=False)
+
+    print("sum bev charger",sum_bev_charger)
+   
+    
+
+    ev_transport_to_bev = sum_ev_transport/snakemake.config['sector']['EV_consumption_1car']*snakemake.config['sector']['bev_charge_rate']*snakemake.config['sector']['bev_charge_avail']
+    print("ev transport to bev",ev_transport_to_bev)
+    bev_charger_ext = n.links[n.links.carrier.str.contains("BEV charger")].query("p_nom_extendable").index
+    pnommax = n.links.loc[bev_charger_ext,"p_nom_max"]
+    pnommax.index = pnommax.index.str.replace(' BEV charger-'+str(year), '', regex=False)
+    add_pnommax = ev_transport_to_bev-sum_bev_charger
+    print("add pnomax",add_pnommax)
+    add_pnommax[add_pnommax < 0] = 0
+    newpnommax = pnommax - sum_bev_charger + add_pnommax
+    print(pnommax, sum_bev_charger, add_pnommax)
+    newpnommax[newpnommax < 0] = 0
+    print(newpnommax)
+    print(len(newpnommax.values))
+    print(newpnommax.values.shape)
+    n.links.loc[bev_charger_ext,'p_nom_max'] = newpnommax.values
+    print(newpnommax)
+
+
 
 def disable_grid_expansion_if_LV_limit_hit(n):
     if not "lv_limit" in n.global_constraints.index:
@@ -177,5 +273,9 @@ def disable_grid_expansion_if_LV_limit_hit(n):
 
     disable_grid_expansion_if_LV_limit_hit(n)
 
+    opts = snakemake.wildcards.sector_opts.split("-")
+    if "T" in opts: # and snakemake.config["sector"]["endogenous_transport"]:
+        adjust_EVs(n, n_p, year)
+
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output[0])
diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py
index e626f4563..d0fc02e4a 100755
--- a/scripts/add_electricity.py
+++ b/scripts/add_electricity.py
@@ -127,6 +127,7 @@ def add_missing_carriers(n, carriers):
     Function to add missing carriers to the network without raising errors.
     """
     missing_carriers = set(carriers) - set(n.carriers.index)
+
     if len(missing_carriers) > 0:
         n.madd("Carrier", missing_carriers)
 
@@ -473,22 +474,38 @@ def attach_conventional_generators(
     # Define generators using modified ppl DataFrame
     caps = ppl.groupby("carrier").p_nom.sum().div(1e3).round(2)
     logger.info(f"Adding {len(ppl)} generators with capacities [GW] \n{caps}")
-
-    n.madd(
-        "Generator",
-        ppl.index,
-        carrier=ppl.carrier,
-        bus=ppl.bus,
-        p_nom_min=ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
-        p_nom=ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
-        p_nom_extendable=ppl.carrier.isin(extendable_carriers["Generator"]),
-        efficiency=ppl.efficiency,
-        marginal_cost=marginal_cost,
-        capital_cost=ppl.capital_cost,
-        build_year=ppl.datein.fillna(0).astype(int),
-        lifetime=(ppl.dateout - ppl.datein).fillna(np.inf),
-        **committable_attrs,
-    )
+    if snakemake.config["foresight"]== "myopic":
+        n.madd(
+            "Generator",
+            ppl.index,
+            carrier=ppl.carrier,
+            bus=ppl.bus,
+            p_nom_min=0, #ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
+            p_nom=ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
+            p_nom_extendable=ppl.carrier.isin(extendable_carriers["Generator"]),
+            efficiency=ppl.efficiency,
+            marginal_cost=marginal_cost,
+            capital_cost=ppl.capital_cost,
+            build_year=ppl.datein.fillna(0).astype(int),
+            lifetime=(ppl.dateout - ppl.datein).fillna(np.inf),
+            **committable_attrs,
+        )
+    else:
+        n.madd(
+            "Generator",
+            ppl.index,
+            carrier=ppl.carrier,
+            bus=ppl.bus,
+            p_nom_min=ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
+            p_nom=ppl.p_nom.where(ppl.carrier.isin(conventional_carriers), 0),
+            p_nom_extendable=ppl.carrier.isin(extendable_carriers["Generator"]),
+            efficiency=ppl.efficiency,
+            marginal_cost=marginal_cost,
+            capital_cost=ppl.capital_cost,
+            build_year=ppl.datein.fillna(0).astype(int),
+            lifetime=(ppl.dateout - ppl.datein).fillna(np.inf),
+            **committable_attrs,
+        )
 
     for carrier in set(conventional_params) & set(carriers):
         # Generators with technology affected
@@ -878,6 +895,10 @@ def attach_line_rating(
     conventional_inputs = {
         k: v for k, v in snakemake.input.items() if k.startswith("conventional_")
     }
+    conventional_inputs2 = {
+        k: v for k, v in snakemake.input.items() if k.startswith("extendable_")
+    }
+    conventional_inputs.update(conventional_inputs2)
 
     if params.conventional["unit_commitment"]:
         unit_commitment = pd.read_csv(snakemake.input.unit_commitment, index_col=0)
@@ -958,6 +979,6 @@ def attach_line_rating(
         )
 
     sanitize_carriers(n, snakemake.config)
-
+   
     n.meta = snakemake.config
     n.export_to_netcdf(snakemake.output[0])
diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py
index e7894324d..6b1233c9f 100644
--- a/scripts/add_existing_baseyear.py
+++ b/scripts/add_existing_baseyear.py
@@ -24,6 +24,7 @@
 from _helpers import update_config_with_sector_opts
 from add_electricity import sanitize_carriers
 from prepare_sector_network import cluster_heat_buses, define_spatial, prepare_costs
+from build_transport_demand import transport_degree_factor
 
 cc = coco.CountryConverter()
 
@@ -123,7 +124,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
     )
 
     df_agg = pd.read_csv(snakemake.input.powerplants, index_col=0)
-
+    
     rename_fuel = {
         "Hard Coal": "coal",
         "Lignite": "lignite",
@@ -221,6 +222,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
     }
 
     for grouping_year, generator in df.index:
+        print(generator)
         # capacity is the capacity in MW at each node for this
         capacity = df.loc[grouping_year, generator]
         capacity = capacity[~capacity.isna()]
@@ -632,6 +634,220 @@ def add_heating_capacities_installed_before_baseyear(
             n.mremove("Link", links_i)
 
 
+def add_land_transport_installed_before_baseyear(
+        n,
+        grouping_year,       
+):
+    airtemp_fn = snakemake.input.temp_air_total
+    temperature = xr.open_dataarray(airtemp_fn).to_pandas()
+    dd_ICE = transport_degree_factor(
+        temperature,
+        options["transport_heating_deadband_lower"],
+        options["transport_heating_deadband_upper"],
+        options["ICE_lower_degree_factor"],
+        options["ICE_upper_degree_factor"],
+    )
+    # transport demand is fulfilled by ICEs for first planning_horizon in myopic pathway
+    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
+    nodes = pop_layout.index
+    # calculate p_nom and profile
+    land_transport_demand_index = n.loads.carrier=='land transport demand'
+    p_set = n.loads_t.p_set.loc[:,land_transport_demand_index]
+    ice_index = n.links.carrier=='land transport oil'
+    #eff_ICE = n.links_t.efficiency.loc[:, ice_index]
+    #eff_ICE.columns = eff_ICE.columns.str.rstrip('-'+str(grouping_year))
+    eff_ICE = n.links.efficiency.loc[ice_index]
+    p_set = p_set.add_suffix(' oil')
+    eff_ICE.index = eff_ICE.index.str.rstrip('-'+str(grouping_year))
+    split_years = costs.at['Liquid fuels ICE (passenger cars)', 'lifetime'] - 1
+    print(grouping_year)
+    year = range(int(grouping_year-split_years), int(grouping_year),1)
+    print(year)
+    pnom = (p_set/eff_ICE).max()/(len(year)+1)
+    #pnom = ((p_set/(eff_ICE/(1+dd_ICE.add_suffix(" land transport oil")))).max()) / len(year) 
+    pnom.index = pnom.index.str.rstrip('-'+str(grouping_year))
+    set_p_nom_ICE = pnom * 0.96
+    p_set_year = p_set/(len(year))
+    #profile = p_set_year.divide(eff_ICE/(1+dd_ICE.add_suffix(" land transport oil")))/pnom/(1+dd_ICE.add_suffix(" land transport oil"))
+    profile = n.links_t.p_min_pu.loc[:,ice_index]
+    
+    set_p_nom_ICE.index = set_p_nom_ICE.index.str.rstrip('land transport oil')
+    profile.columns = profile.columns.str.rstrip('-'+str(grouping_year))
+    profile.columns = profile.columns.str.rstrip('land transport oil')
+    eff_ICE.index = eff_ICE.index.str.rstrip('land transport oil')
+    #eff_ICE.columns = eff_ICE.columns.str.rstrip('land transport oil')
+    #divide ICE build year linearly 
+    for year in year:
+        n.madd(
+        "Link",
+        nodes,
+        suffix=f" land transport oil-{year}",
+        bus0=spatial.oil.nodes,
+        bus1=nodes + " land transport",
+        bus2="co2 atmosphere",
+        carrier="land transport oil",
+        capital_cost = costs.at["Liquid fuels ICE (passenger cars)", "fixed"]/options['ICE_consumption_1car'],
+        efficiency = eff_ICE,
+        efficiency2 = costs.at['oil', 'CO2 intensity'], 
+        lifetime = costs.at['Liquid fuels ICE (passenger cars)', 'lifetime'], 
+        p_nom = set_p_nom_ICE, 
+        p_min_pu = profile, 
+        p_max_pu = profile, 
+        build_year = year
+        )
+        print(year)
+
+
+
+    p_set.columns = p_set.columns.str.rstrip(' oil')
+
+    split_years = costs.at['Liquid fuels ICE (passenger cars)', 'lifetime'] - 1 
+    split_years = 4
+    year = range(int(grouping_year-split_years), int(grouping_year),1)
+
+    ev_index = n.links.carrier=='land transport EV'
+    #eff = n.links_t.efficiency.loc[:, ev_index]
+    eff_EV = n.links.efficiency.loc[ev_index]
+    eff_EV.index = eff_EV.index.str.rstrip('-'+str(grouping_year))
+    #eff = eff/(1+dd_EV.add_suffix(" land transport EV"))
+    
+    pset = p_set.add_suffix(' EV')
+    pnom = (pset.divide(eff_EV)).max()/(len(year)+1)
+    ev_share = 0.04
+    set_p_nom = pnom *ev_share
+    set_p_nom.index = pnom.index.str.rstrip('-'+str(grouping_year))
+    set_p_nom.index = set_p_nom.index.str.rstrip('land transport EV')
+
+    eff_EV.index = eff_EV.index.str.rstrip('land transport EV')
+    bev_availability = options["bev_availability"]
+    for o in opts:
+
+        if "bevavail" not in o:
+            continue
+        oo = o.split("+")
+
+        bev_availability = float(oo[1])
+    pnom_BEV_charger = set_p_nom*options['bev_charge_avail']/options['EV_consumption_1car']*options['bev_charge_rate']
+    pnom_V2G = pnom_BEV_charger*bev_availability
+    enom_EV_storage = set_p_nom*options['bev_charge_avail']/options['EV_consumption_1car']*(options['bev_energy'])*bev_availability
+    availprofile = n.links_t.p_max_pu.loc[:,n.links.carrier=='land transport EV']
+    ev_battery_storage = n.stores[n.stores.carrier.str.contains("Li ion|EV battery storage")]
+    ev_battery_storage_index = ev_battery_storage.index
+    dsmprofile = n.stores_t.e_min_pu.loc[:,ev_battery_storage_index]
+    dsmprofile.columns = dsmprofile.columns.str.rstrip('-'+str(grouping_year))
+    dsmmax = n.stores.e_nom_max.loc[ev_battery_storage_index]
+    dsmmax.index = dsmmax.index.str.rstrip('-'+str(grouping_year))
+    p_nom_max = n.links.p_nom_max.loc[n.links.carrier=='BEV charger']
+    p_nom_max.index = p_nom_max.index.str.replace(' BEV charger-'+str(grouping_year),"",regex=False)
+    p_nom_max_V2G = p_nom_max
+
+    for o in opts:
+        if o=="V2G+p0":
+            p_nom_max_V2G = p_nom_max_V2G*0
+            pnom_V2G = pnom_V2G*0
+        if o=="EV battery storage+x0":
+            enom_EV_storage = enom_EV_storage*0
+    for year in year:
+
+        n.madd(
+            "Link",
+            nodes,
+            suffix=f" land transport EV-{year}",
+            bus0=nodes + " EV battery",                              
+            bus1 =nodes + " land transport",
+            carrier="land transport EV",
+            lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+            capital_cost = costs.at["Battery electric (passenger cars)", "fixed"]/options['EV_consumption_1car'],
+            efficiency = eff_EV,     
+            p_nom = set_p_nom, 
+            p_min_pu = profile, 
+            p_max_pu = profile, 
+            build_year = year
+        )
+
+        n.madd(
+            "Link",
+            nodes,
+            suffix=f" BEV charger-{year}",
+            bus0=nodes + " low voltage",
+            bus1=nodes + " EV battery",
+            p_nom= pnom_BEV_charger,
+            p_nom_max = p_nom_max,
+            carrier="BEV charger",
+            p_max_pu=availprofile,
+            efficiency=options.get("bev_charge_efficiency", 0.9),
+            lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+            build_year = year
+        )
+
+        if options["bev_dsm"]:
+            
+
+            n.madd(
+                "Store",
+                nodes,
+                suffix=f" EV battery storage-{year}",
+                bus=nodes + " EV battery",
+                carrier="EV battery storage", #"Li ion",
+                e_cyclic=True,
+                e_max_pu=1,
+                e_min_pu=dsmprofile.values,
+                e_nom = enom_EV_storage,
+                e_nom_max = dsmmax,
+                lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+                build_year = year
+            )
+
+        if options["v2g"]:
+            n.madd(
+                "Link",
+                nodes,
+                suffix=f" V2G-{year}",
+                bus1=nodes + " low voltage",
+                bus0=nodes + " EV battery",
+                carrier="V2G",
+                p_nom = pnom_V2G,
+                p_nom_max =p_nom_max_V2G*bev_availability,
+                p_max_pu = availprofile,
+                efficiency=options.get("bev_charge_efficiency", 0.9),
+                lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+                build_year = year
+            )
+
+
+    print(n.links.loc[(n.links.filter(like="land transport oil-"+str(year+1),axis=0)).index,'p_nom'])
+
+    n.links.loc[(n.links.filter(like="land transport oil-"+str(year+1),axis=0)).index,'p_nom'] = set_p_nom_ICE.values
+    n.links.loc[(n.links.filter(like="land transport EV-"+str(year+1),axis=0)).index,'p_nom'] = set_p_nom.values
+    n.links.loc[(n.links.filter(like="BEV charger-"+str(year+1),axis=0)).index,'p_nom'] = pnom_BEV_charger.values
+    n.links.loc[(n.links.filter(like="V2G-"+str(year+1),axis=0)).index,'p_nom'] = pnom_V2G.values
+    n.stores.loc[n.stores.filter(like="EV battery storage"+str(year+1),axis=0).index, 'e_nom'] = enom_EV_storage
+    print(n.links.loc[(n.links.filter(like="land transport oil-"+str(year+1),axis=0)).index,'p_nom'])
+    n.links.loc[(n.links.filter(like="land transport oil-"+str(year+1),axis=0)).index,'p_nom_extendable'] = False
+    n.links.loc[(n.links.filter(like="land transport EV-"+str(year+1),axis=0)).index, 'p_nom_extendable'] = False
+    n.links.loc[(n.links.filter(like="BEV charger-"+str(year+1),axis=0)).index, 'p_nom_extendable'] = False
+    n.links.loc[(n.links.filter(like="V2G-"+str(year+1),axis=0)).index, 'p_nom_extendable'] = False
+
+    ev_battery_storage = n.stores[n.stores.carrier.str.contains("Li ion|EV battery storage")]
+    ev_battery_storage_index = ev_battery_storage.index
+    n.stores.loc[ev_battery_storage_index, "e_nom_extendable"] = False
+    print("4",n.stores.loc[n.stores.filter(like="EV battery storage",axis=0).index, 'e_nom_extendable'])
+
+    number_cars = pd.read_csv(snakemake.input.existing_transport, index_col=0)[
+        "number cars"
+    ]
+
+    print(n.stores.loc[n.stores.filter(like="EV battery storage",axis=0).index, 'e_nom'])
+    print(n.links.p_nom_max.loc[n.links.carrier=="V2G"])
+    print(n.links.p_nom[n.links.carrier.str.contains("transport oil")].sum()*0.3)
+    print(n.links.p_nom[n.links.carrier.str.contains("transport EV")].sum()*0.99)
+    print(n.loads_t.p_set[n.loads[n.loads.carrier.str.contains("land ")].index].sum(axis=1))
+
+    n.links.loc[(n.links.filter(like="land transport fuel cell-"+str(year+1),axis=0)).index, 'p_nom_extendable'] = False
+    
+    logger.info(f"The model only assumes the minimum number of ICE vehicles to supply the transport demand which represents {round(int(snakemake.config['scenario']['clusters'][0])*p_set.max(axis=0).sum()/number_cars.sum()*100,2)}% of the existing ICE vehicle fleet in 2020")
+
+
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
@@ -699,6 +915,12 @@ def add_heating_capacities_installed_before_baseyear(
             default_lifetime,
         )
 
+    endo_transport = ( (options["land_transport_electric_share"][baseyear] is None ) 
+        and (options["land_transport_fuel_cell_share"][baseyear] is None)
+        and (options["land_transport_ice_share"][baseyear] is None))
+    if "T" in opts and options["endogenous_transport"] and endo_transport and snakemake.config["foresight"]=="myopic":
+        add_land_transport_installed_before_baseyear(n, baseyear)
+        
     if options.get("cluster_heat_buses", False):
         cluster_heat_buses(n)
 
diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py
index f9db92711..5ac6ca321 100644
--- a/scripts/build_renewable_profiles.py
+++ b/scripts/build_renewable_profiles.py
@@ -296,6 +296,14 @@
     duration = time.time() - start
     logger.info(f"Completed landuse availability calculation ({duration:2.2f}s)")
 
+    # For Moldova and Ukraine: Overwrite parts not covered by Corine with
+    # externally determined available areas
+    if "availability_matrix_MD_UA" in snakemake.input.keys():
+        availability_MDUA = xr.open_dataarray(
+            snakemake.input["availability_matrix_MD_UA"]
+        )
+        availability.loc[availability_MDUA.coords] = availability_MDUA
+
     # For Moldova and Ukraine: Overwrite parts not covered by Corine with
     # externally determined available areas
     if "availability_matrix_MD_UA" in snakemake.input.keys():
diff --git a/scripts/build_retro_cost.py b/scripts/build_retro_cost.py
index d2aae1404..cca3c99c9 100755
--- a/scripts/build_retro_cost.py
+++ b/scripts/build_retro_cost.py
@@ -208,8 +208,7 @@ def prepare_building_stock_data():
 
     # heated floor area ----------------------------------------------------------
     area = building_data[
-        (building_data.type == "Heated area [Mm²]")
-        & (building_data.subsector != "Total")
+        (building_data.type == "Heated area [Mm²]") & (building_data.detail != "Total")
     ]
     area_tot = area[["country", "sector", "value"]].groupby(["country", "sector"]).sum()
     area = pd.concat(
diff --git a/scripts/build_transport_demand.py b/scripts/build_transport_demand.py
index 0bcfb7ed6..f77bafed7 100644
--- a/scripts/build_transport_demand.py
+++ b/scripts/build_transport_demand.py
@@ -15,13 +15,16 @@
 
 
 def build_nodal_transport_data(fn, pop_layout):
+    # get numbers of car and fuel efficieny per country
     transport_data = pd.read_csv(fn, index_col=0)
 
+    # break number of cars down to nodal level based on population density
     nodal_transport_data = transport_data.loc[pop_layout.ct].fillna(0.0)
     nodal_transport_data.index = pop_layout.index
     nodal_transport_data["number cars"] = (
         pop_layout["fraction"] * nodal_transport_data["number cars"]
     )
+    # fill missing fuel efficiency with average data
     nodal_transport_data.loc[
         nodal_transport_data["average fuel efficiency"] == 0.0,
         "average fuel efficiency",
@@ -33,8 +36,10 @@ def build_nodal_transport_data(fn, pop_layout):
 def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data):
     ## Get overall demand curve for all vehicles
 
+    # averaged weekly counts from the year 2010-2015
     traffic = pd.read_csv(traffic_fn, skiprows=2, usecols=["count"]).squeeze("columns")
 
+    # create annual profile take account time zone + summer time
     transport_shape = generate_periodic_profiles(
         dt_index=snapshots,
         nodes=nodes,
@@ -42,15 +47,6 @@ def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data):
     )
     transport_shape = transport_shape / transport_shape.sum()
 
-    # electric motors are more efficient, so alter transport demand
-
-    plug_to_wheels_eta = options["bev_plug_to_wheel_efficiency"]
-    battery_to_wheels_eta = plug_to_wheels_eta * options["bev_charge_efficiency"]
-
-    efficiency_gain = (
-        nodal_transport_data["average fuel efficiency"] / battery_to_wheels_eta
-    )
-
     # get heating demand for correction to demand time series
     temperature = xr.open_dataarray(airtemp_fn).to_pandas()
 
@@ -74,17 +70,49 @@ def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data):
     # divide out the heating/cooling demand from ICE totals
     # and multiply back in the heating/cooling demand for EVs
     ice_correction = (transport_shape * (1 + dd_ICE)).sum() / transport_shape.sum()
-
+    print("ice correction",ice_correction)
     energy_totals_transport = (
         pop_weighted_energy_totals["total road"]
         + pop_weighted_energy_totals["total rail"]
         - pop_weighted_energy_totals["electricity rail"]
     )
+    plug_to_wheels_eta = options["bev_plug_to_wheel_efficiency"]
+    battery_to_wheels_eta = plug_to_wheels_eta * options["bev_charge_efficiency"]
 
-    return (
-        (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears)
+    efficiency_gain = (
+        nodal_transport_data["average fuel efficiency"] / battery_to_wheels_eta
+    )
+    print("only transport",(energy_totals_transport).sum().sum())
+    print("transport",(transport_shape.multiply(energy_totals_transport)).sum().sum() * 1e6 * nyears)
+    print("previous", (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears)
         .divide(efficiency_gain * ice_correction)
-        .multiply(1 + dd_EV)
+        .multiply(1 + dd_EV))
+    print(efficiency_gain)
+    print(nodal_transport_data["average fuel efficiency"])
+    print("current",(transport_shape.multiply(energy_totals_transport) * 1e9 * nyears)
+       .multiply(nodal_transport_data["average fuel efficiency"]*(1+dd_ICE)))
+    
+    print("try",(transport_shape.multiply(energy_totals_transport) * nyears)
+       .divide(1e3 * nodal_transport_data["average fuel efficiency"]*(1+dd_ICE)))
+    
+    print("try",(transport_shape.multiply(energy_totals_transport) * 1e6 * nyears * options["transport_internal_combustion_efficiency"])
+       .divide((ice_correction)))
+    
+
+    A = (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears * options["transport_internal_combustion_efficiency"]).divide((ice_correction))
+
+    print(A.sum().sum())
+
+    print("years", nyears)
+    # calculate actual energy demand independent of ICE
+    # return (
+    #     (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears)
+    #    .divide(nodal_transport_data["average fuel efficiency"]*ice_correction)
+    # )
+    return (
+        (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears * options["transport_internal_combustion_efficiency"])
+       #.divide((ice_correction))
+   
     )
 
 
@@ -121,11 +149,15 @@ def bev_availability_profile(fn, snapshots, nodes, options):
     """
     Derive plugged-in availability for passenger electric vehicles.
     """
+    # car count in typical week
     traffic = pd.read_csv(fn, skiprows=2, usecols=["count"]).squeeze("columns")
 
+    # maximum share plugged-in availability for passenger electric vehicles
     avail_max = options["bev_avail_max"]
+    # average share plugged-in availability for passenger electric vehicles
     avail_mean = options["bev_avail_mean"]
 
+    # linear scaling, highest when traffic is lowest, decreases if traffic increases
     avail = avail_max - (avail_max - avail_mean) * (traffic - traffic.min()) / (
         traffic.mean() - traffic.min()
     )
@@ -140,6 +172,8 @@ def bev_availability_profile(fn, snapshots, nodes, options):
 def bev_dsm_profile(snapshots, nodes, options):
     dsm_week = np.zeros((24 * 7,))
 
+    # assuming that at a certain time ("bev_dsm_restriction_time") EVs have to
+    # be charged to a minimum value (defined in bev_dsm_restriction_value)
     dsm_week[(np.arange(0, 7, 1) * 24 + options["bev_dsm_restriction_time"])] = options[
         "bev_dsm_restriction_value"
     ]
@@ -158,7 +192,7 @@ def bev_dsm_profile(snapshots, nodes, options):
         snakemake = mock_snakemake(
             "build_transport_demand",
             simpl="",
-            clusters=48,
+            clusters=37,
         )
 
     pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py
index 28f08396e..b90a7ca50 100644
--- a/scripts/cluster_network.py
+++ b/scripts/cluster_network.py
@@ -234,7 +234,7 @@ def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="cbc"):
 
     assert (
         n_clusters >= len(N) and n_clusters <= N.sum()
-    ), f"Number of clusters must be {len(N)} <= n_clusters <= {N.sum()} for this selection of countries."
+    ), f"Number of clusters {n_clusters} must be {len(N)} <= n_clusters <= {N.sum()} for this selection of countries."
 
     if isinstance(focus_weights, dict):
         total_focus = sum(list(focus_weights.values()))
@@ -502,7 +502,6 @@ def plot_busmap_for_n_clusters(n, n_clusters, fn=None):
                 carriers += [f"{c} {label} efficiency" for label in labels]
                 n.generators.carrier.update(gens.carrier + " " + suffix + " efficiency")
         aggregate_carriers = carriers
-
     if n_clusters == len(n.buses):
         # Fast-path if no clustering is necessary
         busmap = n.buses.index.to_series()
diff --git a/scripts/plot_network.py b/scripts/plot_network.py
index 674811208..cac56d0f3 100644
--- a/scripts/plot_network.py
+++ b/scripts/plot_network.py
@@ -1098,8 +1098,8 @@ def plot_map_perfect(
             transmission=False,
         )
 
-        plot_h2_map(n, regions)
-        plot_ch4_map(n)
+        #plot_h2_map(n, regions)
+        #plot_ch4_map(n)
         plot_map_without(n)
 
     # plot_series(n, carrier="AC", name=suffix)
diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py
index 67ac9b553..02dd68f9f 100644
--- a/scripts/plot_summary.py
+++ b/scripts/plot_summary.py
@@ -443,7 +443,7 @@ def historical_emissions(countries):
     return emissions
 
 
-def plot_carbon_budget_distribution(input_eurostat):
+def plot_carbon_budget_distribution(input_eurostat,opts):
     """
     Plot historical carbon emissions in the EU and decarbonization path.
     """
@@ -479,8 +479,10 @@ def plot_carbon_budget_distribution(input_eurostat):
 
     if snakemake.config["foresight"] == "myopic":
         path_cb = "results/" + snakemake.params.RDIR + "/csvs/"
-        co2_cap = pd.read_csv(path_cb + "carbon_budget_distribution.csv", index_col=0)[
-            ["cb"]
+        print(opts[-1])
+        print(path_cb + "carbon_budget_distribution"+ str(opts[-1]) +".csv")
+        co2_cap = pd.read_csv(path_cb + "carbon_budget_distribution"+ str(opts[-1]) +".csv", index_col=0)[
+            [opts[-1]]
         ]
         co2_cap *= e_1990
     else:
@@ -581,4 +583,4 @@ def plot_carbon_budget_distribution(input_eurostat):
     for sector_opts in snakemake.params.sector_opts:
         opts = sector_opts.split("-")
         if any("cb" in o for o in opts) or snakemake.config["foresight"] == "perfect":
-            plot_carbon_budget_distribution(snakemake.input.eurostat)
+            plot_carbon_budget_distribution(snakemake.input.eurostat,opts)
diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py
index b57b832c4..f925eb2d1 100755
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@@ -20,6 +20,7 @@
 from _helpers import generate_periodic_profiles, update_config_with_sector_opts
 from add_electricity import calculate_annuity, sanitize_carriers
 from build_energy_totals import build_co2_totals, build_eea_co2, build_eurostat_co2
+from build_transport_demand import transport_degree_factor
 from networkx.algorithms import complement
 from networkx.algorithms.connectivity.edge_augmentation import k_edge_augmentation
 from pypsa.geo import haversine_pts
@@ -288,10 +289,14 @@ def build_carbon_budget(o, input_eurostat, fn, emissions_scope, report_year):
         year=2018,
     )
 
+
     planning_horizons = snakemake.params.planning_horizons
     if not isinstance(planning_horizons, list):
         planning_horizons = [planning_horizons]
     t_0 = planning_horizons[0]
+    if  t_0 > 2020 :
+        print('changing emissions at beginning of path from', e_0, ' to', 3.2)    
+        e_0 = 3.2
 
     if "be" in o:
         # final year in the path
@@ -309,12 +314,18 @@ def beta_decay(t):
         m = (1 + np.sqrt(1 + r * T)) / T
 
         def exponential_decay(t):
-            return (e_0 / e_1990) * (1 + (m + r) * (t - t_0)) * np.exp(-m * (t - t_0))
+            if t == planning_horizons[-1] and snakemake.params.foresight in ["myopic", "perfect"]:
+                return 0
+            else:
+                return (e_0 / e_1990) * (1 + (m + r) * (t - t_0)) * np.exp(-m * (t - t_0))
 
         co2_cap = pd.Series(
             {t: exponential_decay(t) for t in planning_horizons}, name=o
         )
 
+    if 2050 in co2_cap.index:
+        co2_cap.loc[2050] = 0
+
     # TODO log in Snakefile
     csvs_folder = fn.rsplit("/", 1)[0]
     if not os.path.exists(csvs_folder):
@@ -501,14 +512,14 @@ def add_carrier_buses(n, carrier, nodes=None):
         carrier=carrier,
         capital_cost=capital_cost,
     )
-
+    
     n.madd(
         "Generator",
         nodes,
         bus=nodes,
         p_nom_extendable=True,
         carrier=carrier,
-        marginal_cost=costs.at[carrier, "fuel"],
+        marginal_cost= costs.at[carrier, "fuel"]*margoil if carrier=="oil" else costs.at[carrier, "fuel"],
     )
 
 
@@ -1476,20 +1487,81 @@ def add_land_transport(n, costs):
     fuel_cell_share = get(options["land_transport_fuel_cell_share"], investment_year)
     electric_share = get(options["land_transport_electric_share"], investment_year)
     ice_share = get(options["land_transport_ice_share"], investment_year)
+    
+    if options["endogenous_transport"]==False:
+        if any([fuel_cell_share == None, electric_share == None, ice_share == None]):
+            logger.warning(
+                f"Exogenous land transport selected, but not all land transport shares are defined."
+            )
+        total_share = fuel_cell_share + electric_share + ice_share
+        if total_share != 1:
+            logger.warning(
+                f"Total land transport shares sum up to {total_share:.2%}, corresponding to increased or decreased demand assumptions."
+            )
 
-    total_share = fuel_cell_share + electric_share + ice_share
-    if total_share != 1:
-        logger.warning(
-            f"Total land transport shares sum up to {total_share:.2%}, corresponding to increased or decreased demand assumptions."
-        )
 
-    logger.info(f"FCEV share: {fuel_cell_share*100}%")
-    logger.info(f"EV share: {electric_share*100}%")
-    logger.info(f"ICEV share: {ice_share*100}%")
+    if options["endogenous_transport"]==True:
+        if all([fuel_cell_share != None, electric_share != None, ice_share != None]):
+            logger.warning(
+                f"Endogenous land transport is activated, but land transport shares are exogenously defined for the three types of land transport. The exogenous shares are enforced."
+            )
+        elif any([fuel_cell_share != None, electric_share != None, ice_share != None]):
+            logger.warning(
+                f"Endogenous land transport is activated, but land transport shares are exogenously defined. The exogenous shares are enforced, limiting the endogenous calculation of optimal shares"
+            )
+
+    if fuel_cell_share != None:
+        logger.info(f"FCEV share: {fuel_cell_share*100}%")
+    if electric_share != None:
+        logger.info(f"EV share: {electric_share*100}%")
+    if ice_share != None:
+        logger.info(f"ICEV share: {ice_share*100}%")
 
     nodes = pop_layout.index
 
-    if electric_share > 0:
+    # Add load for transport demand
+    n.add("Carrier", "land transport demand")
+
+    n.madd("Bus", 
+            nodes,
+            location=nodes,
+            suffix=" land transport",
+            carrier="land transport demand")
+    
+    p_set = transport[nodes]
+
+
+    n.madd(
+        "Load",
+        nodes,
+        suffix=" land transport",
+        bus=nodes + " land transport",
+        carrier="land transport demand",
+        p_set=p_set, 
+    )
+
+    # get heating demand for correction to demand time series
+    airtemp_fn = snakemake.input.temp_air_total
+    temperature = xr.open_dataarray(airtemp_fn).to_pandas()
+    dd_ICE = transport_degree_factor(
+        temperature,
+        options["transport_heating_deadband_lower"],
+        options["transport_heating_deadband_upper"],
+        options["ICE_lower_degree_factor"],
+        options["ICE_upper_degree_factor"],
+    )
+
+    dd_EV = transport_degree_factor(
+        temperature,
+        options["transport_heating_deadband_lower"],
+        options["transport_heating_deadband_upper"],
+        options["EV_lower_degree_factor"],
+        options["EV_upper_degree_factor"],
+    )
+
+    # Add EV links
+    if electric_share != 0:
+        
         n.add("Carrier", "Li ion")
 
         n.madd(
@@ -1500,131 +1572,235 @@ def add_land_transport(n, costs):
             carrier="Li ion",
             unit="MWh_el",
         )
-
-        p_set = (
-            electric_share
-            * (
-                transport[nodes]
-                + cycling_shift(transport[nodes], 1)
-                + cycling_shift(transport[nodes], 2)
-            )
-            / 3
-        )
-
-        n.madd(
-            "Load",
-            nodes,
-            suffix=" land transport EV",
-            bus=nodes + " EV battery",
-            carrier="land transport EV",
-            p_set=p_set,
-        )
-
-        p_nom = number_cars * options.get("bev_charge_rate", 0.011) * electric_share
-
+        
+        p_availEV = number_cars * options.get("bev_charge_rate", 0.011)
         n.madd(
             "Link",
             nodes,
             suffix=" BEV charger",
             bus0=nodes,
             bus1=nodes + " EV battery",
-            p_nom=p_nom,
+            p_nom_extendable = True,
+            p_nom_max = p_availEV*3.3,
+            #p_nom = p_availEV,
             carrier="BEV charger",
             p_max_pu=avail_profile[nodes],
             efficiency=options.get("bev_charge_efficiency", 0.9),
-            # These were set non-zero to find LU infeasibility when availability = 0.25
-            # p_nom_extendable=True,
-            # p_nom_min=p_nom,
-            # capital_cost=1e6,  #i.e. so high it only gets built where necessary
+            lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
         )
 
-    if electric_share > 0 and options["v2g"]:
-        n.madd(
-            "Link",
-            nodes,
-            suffix=" V2G",
-            bus1=nodes,
-            bus0=nodes + " EV battery",
-            p_nom=p_nom,
-            carrier="V2G",
-            p_max_pu=avail_profile[nodes],
-            efficiency=options.get("bev_charge_efficiency", 0.9),
-        )
+        # Only when v2g option is True, add the v2g link
+        if options["v2g"]:
+            n.madd(
+                "Link",
+                nodes,
+                suffix=" V2G",
+                bus1=nodes,
+                bus0=nodes + " EV battery",
+                p_nom_extendable = True,
+                carrier="V2G",
+                p_max_pu=avail_profile[nodes],
+                p_nom_max = p_availEV*3.3*bev_availability,
+                #p_nom = p_availEV,
+                efficiency=options.get("bev_charge_efficiency", 0.9),
+                lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+            )
 
-    if electric_share > 0 and options["bev_dsm"]:
-        e_nom = (
-            number_cars
-            * options.get("bev_energy", 0.05)
-            * options["bev_availability"]
-            * electric_share
-        )
+        # add bev management when enabled
+        if options["bev_dsm"]:
+            
+            n.add("Carrier", "EV battery storage")
 
+            n.madd(
+                "Store",
+                nodes,
+                suffix=" EV battery storage",
+                bus=nodes + " EV battery",
+                carrier="EV battery storage", #"Li ion",
+                e_cyclic=True,
+                e_nom_extendable=True,
+                e_max_pu=1,
+                e_min_pu=dsm_profile[nodes],
+                e_nom_max = number_cars.values*options.get('bev_energy')*3.3,  #options.get("bev_availability"),
+                lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+            )
+        
+        eff = options["transport_electric_vehicle_efficiency"] #costs.at['Battery electric (passenger cars)', 'efficiency'] #/(1+dd_EV)
         n.madd(
-            "Store",
+            "Link",
             nodes,
-            suffix=" battery storage",
-            bus=nodes + " EV battery",
-            carrier="battery storage",
-            e_cyclic=True,
-            e_nom=e_nom,
-            e_max_pu=1,
-            e_min_pu=dsm_profile[nodes],
+            suffix=" land transport EV",
+            bus0=nodes + " EV battery",                              
+            bus1 =nodes + " land transport",
+            carrier="land transport EV",
+            lifetime = costs.at['Battery electric (passenger cars)', 'lifetime'], 
+            capital_cost = costs.at["Battery electric (passenger cars)", "fixed"]/options['EV_consumption_1car'],
+            efficiency = eff,     
+            p_nom_extendable=True,
         )
 
-    if fuel_cell_share > 0:
+        # exogenous pathway if electric share is defined
+        if float(electric_share or 0) > 0:
+            # set p_nom to fulfill electric share, capital cost to 0
+            ev_index = n.links.carrier=='land transport EV'
+            n.links.loc[ev_index, "p_nom_extendable"] = False
+
+            # available EVs for land transport demand
+            pnom = (p_set.divide(eff)).max()
+            pnom = pnom.add_suffix(' land transport EV')
+            n.links.loc[ev_index,'p_nom'] = pnom * electric_share
+            n.links.loc[ev_index, "capital_cost"] = 0
+            
+            # available EVs for charging/discharging
+            p_availEV = number_cars * options.get("bev_charge_rate", 0.011)
+            bev_charger_index = n.links.carrier=="BEV charger"
+            n.links.loc[bev_charger_index, "p_nom_extendable"] = False
+            n.links.loc[bev_charger_index, "p_nom"] = p_availEV.add_suffix(' BEV charger')
+
+                  
+            
+            if options["v2g"]:
+                v2g_index = n.links.carrier=="V2G"
+                n.links.loc[v2g_index, "p_nom_extendable"] = False
+                n.links.loc[v2g_index, "p_nom"] = p_availEV.add_suffix(' V2G') * electric_share 
+                            
+            if options["bev_dsm"]:
+                ev_battery_storage_index = n.stores.carrier==("Li ion" or "EV battery storage")
+                n.stores.loc[ev_battery_storage_index, "e_nom_extendable"] = False
+                e_nom = number_cars* options.get("bev_energy", 0.05) * bev_availability * electric_share #options["bev_availability"] 
+                n.stores.loc[ev_battery_storage_index, "e_nom"] = e_nom.add_suffix(" EV battery storage")
+
+
+    # Add hydrogen vehicle links
+    if fuel_cell_share != 0:
+        eff = options["transport_fuel_cell_efficiency"] #/(1+dd_ICE)
+        
         n.madd(
-            "Load",
+            "Link",
             nodes,
             suffix=" land transport fuel cell",
-            bus=nodes + " H2",
+            bus0=nodes + " H2",
+            bus1=nodes + " land transport",
             carrier="land transport fuel cell",
-            p_set=fuel_cell_share
-            / options["transport_fuel_cell_efficiency"]
-            * transport[nodes],
-        )
-
-    if ice_share > 0:
-        add_carrier_buses(n, "oil")
-
-        ice_efficiency = options["transport_internal_combustion_efficiency"]
+            lifetime = costs.at['Hydrogen fuel cell (passenger cars)', 'lifetime'], 
+            efficiency=  eff, 
+            capital_cost = costs.at["Hydrogen fuel cell (passenger cars)", "fixed"]/options["H2_consumption_1car"],
+            p_nom_extendable=True,
 
-        p_set_land_transport_oil = (
-            ice_share
-            / ice_efficiency
-            * transport[nodes].rename(columns=lambda x: x + " land transport oil")
         )
+        # exogenous pathway if fuel_cell_share is defined
+        if float(fuel_cell_share or 0) > 0:
+            pnom = (p_set.divide(eff)).max()
+            pnom=pnom.add_suffix(' land transport fuel cell')
+            fuel_cell_index = n.links.carrier=='land transport fuel cell'
+            n.links.loc[fuel_cell_index, "p_nom_extendable"] = False
+            n.links.loc[fuel_cell_index,'p_nom'] = fuel_cell_share * pnom
+            n.links.loc[fuel_cell_index, "capital_cost"] = 0
+           
 
-        if not options["regional_oil_demand"]:
-            p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame(
-                name="EU land transport oil"
-            )
-
-        n.madd(
-            "Bus",
-            spatial.oil.land_transport,
-            location=spatial.oil.demand_locations,
-            carrier="land transport oil",
-            unit="land transport",
-        )
-
-        n.madd(
-            "Load",
-            spatial.oil.land_transport,
-            bus=spatial.oil.land_transport,
-            carrier="land transport oil",
-            p_set=p_set_land_transport_oil,
-        )
+    if ice_share != 0:
+        add_carrier_buses(n, "oil")
 
+        eff = options["transport_internal_combustion_efficiency"] #/(1+dd_ICE)
+        
         n.madd(
             "Link",
-            spatial.oil.land_transport,
+            nodes,
+            suffix=" land transport oil",
             bus0=spatial.oil.nodes,
-            bus1=spatial.oil.land_transport,
+            bus1=nodes + " land transport",
             bus2="co2 atmosphere",
             carrier="land transport oil",
-            efficiency2=costs.at["oil", "CO2 intensity"],
-            p_nom_extendable=True,
-        )
+            lifetime = costs.at['Liquid fuels ICE (passenger cars)', 'lifetime'], 
+            capital_cost = costs.at["Liquid fuels ICE (passenger cars)", "fixed"]/options['ICE_consumption_1car'],
+            efficiency = eff,  
+            efficiency2 = costs.at["oil", "CO2 intensity"],
+            p_nom_extendable=True, 
+        )
+
+        # exogenous pathway if ice_share is defined
+        if float(ice_share or 0) > 0:
+            pnom = (p_set.divide(eff)).max()
+            pnom = pnom.add_suffix(' land transport oil')
+            ice_index = n.links.carrier=='land transport oil'
+            n.links.loc[ice_index, "p_nom_extendable"] = False
+            n.links.loc[ice_index ,'p_nom'] = pnom
+            n.links.loc[ice_index, "capital_cost"] = 0
+
+
+def adjust_land_transport(n):
+
+    # get heating demand for correction to demand time series
+    airtemp_fn = snakemake.input.temp_air_total
+    temperature = xr.open_dataarray(airtemp_fn).to_pandas()
+    dd_ICE = transport_degree_factor(
+        temperature,
+        options["transport_heating_deadband_lower"],
+        options["transport_heating_deadband_upper"],
+        options["ICE_lower_degree_factor"],
+        options["ICE_upper_degree_factor"],
+    )
+
+    dd_EV = transport_degree_factor(
+        temperature,
+        options["transport_heating_deadband_lower"],
+        options["transport_heating_deadband_upper"],
+        options["EV_lower_degree_factor"],
+        options["EV_upper_degree_factor"],
+    )
+    # Adjust p_min_pu and p_max_pu after temporal aggregation, because the temporal changing efficiencies 
+    # lead to wrong calculations otherwise
+    land_transport_demand_index = n.loads.carrier=='land transport demand'
+    p_set = n.loads_t.p_set.loc[:,land_transport_demand_index]
+    fuel_cell_share = get(options["land_transport_fuel_cell_share"], investment_year)
+    electric_share = get(options["land_transport_electric_share"], investment_year)
+    ice_share = get(options["land_transport_ice_share"], investment_year)
+    if ice_share != 0:
+        name = "land transport oil"
+        ice_index = n.links.carrier==name
+        #eff = n.links_t.efficiency.loc[:, ice_index]
+        #print("eff", eff)
+        eff = n.links.efficiency.loc[ice_index]
+        #eff = eff/(1+dd_ICE.add_suffix(" land transport oil"))
+        pset = p_set.add_suffix(' oil')
+        pnom = (pset.divide(eff)).max()
+        profile=pset.divide(eff)/pnom #/(1+dd_ICE.add_suffix(" land transport oil"))
+        n.links_t.p_min_pu[n.links[ice_index].index] = profile
+        #n.links_t.p_max_pu[n.links[ice_index].index] = profile
+        if float(ice_share or 0) > 0:
+            n.links.loc[ice_index,'p_nom'] = ice_share * pnom 
+            #n.links_t.p_max_pu[n.links[ice_index].index] = profile
+            
+    if fuel_cell_share != 0:
+        fuel_cell_index = n.links.carrier=='land transport fuel cell'
+        #eff = n.links_t.efficiency.loc[:, fuel_cell_index]
+        eff = n.links.efficiency.loc[fuel_cell_index]
+        #eff = eff/(1+dd_ICE.add_suffix(" land transport fuel cell"))
+        pset = p_set.add_suffix(' fuel cell')
+        pnom = (pset.divide(eff)).max()
+        pu=pset.divide(eff)/pnom #/(1+dd_ICE.add_suffix(" land transport fuel cell"))
+        n.links_t.p_min_pu[n.links[fuel_cell_index].index] = pu
+        #n.links_t.p_max_pu[n.links[fuel_cell_index].index] = pu
+        if float(fuel_cell_share or 0) > 0:
+            n.links.loc[fuel_cell_index,'p_nom'] = fuel_cell_share * pnom 
+            #n.links_t.p_max_pu[n.links[fuel_cell_index].index] = pu
+
+    # EV link for endogenous transport is not constraint because constraining all land transport link 
+    # with p_min_pu/p_max_pu while p_nom_extentable=True leads to solving infeasibility by gurobi      
+    if electric_share!= 0:
+        ev_index = n.links.carrier=='land transport EV'
+        #eff = n.links_t.efficiency.loc[:, ev_index]
+        eff = n.links.efficiency.loc[ev_index]
+        #eff = eff/(1+dd_EV.add_suffix(" land transport EV"))
+        pset = p_set.add_suffix(' EV')
+        pnom = (pset.divide(eff)).max()
+        pu=pset.divide(eff)/pnom #/(1+dd_EV.add_suffix(" land transport EV"))
+        n.links_t.p_min_pu[n.links[ev_index].index] = pu
+        #n.links_t.p_max_pu[n.links[ev_index].index] = pu 
+        if float(electric_share or 0) > 0:
+            n.links.loc[ev_index,'p_nom'] = electric_share * pnom 
+            #n.links_t.p_min_pu[n.links[ev_index].index] = pu
+            #n.links_t.p_max_pu[n.links[ev_index].index] = pu   
 
 
 def build_heat_demand(n):
@@ -3305,10 +3481,11 @@ def remove_h2_network(n):
 
 def maybe_adjust_costs_and_potentials(n, opts):
     for o in opts:
-        flags = ["+e", "+p", "+m", "+c"]
+        flags = ["+e", "+p", "+m", "+c", "+n", "+l", "+x", "+k","+q"]
         if all(flag not in o for flag in flags):
             continue
         oo = o.split("+")
+        
         carrier_list = np.hstack(
             (
                 n.generators.carrier.unique(),
@@ -3325,17 +3502,27 @@ def maybe_adjust_costs_and_potentials(n, opts):
                 "e": "e_nom_max",
                 "c": "capital_cost",
                 "m": "marginal_cost",
+                "k": "efficiency",
+                "n": "efficiency2",
+                "l": "efficiency3",
+                "q": "p_nom",
+                "x": "e_nom_extendable"
             }
             attr = attr_lookup[oo[1][0]]
+
             factor = float(oo[1][1:])
             # beware if factor is 0 and p_nom_max is np.inf, 0*np.inf is nan
             if carrier == "AC":  # lines do not have carrier
                 n.lines[attr] *= factor
             else:
-                if attr == "p_nom_max":
+                if attr == "p_nom_max" or attr=="p_nom":
                     comps = {"Generator", "Link", "StorageUnit"}
-                elif attr == "e_nom_max":
+                elif attr == "e_nom_max" or attr == "e_nom_extendable":
                     comps = {"Store"}
+                    if attr == "e_nom_extendable":
+                        factor = bool(factor)
+                elif attr == "efficiency2" or attr == "efficiency3" or attr == "efficiency":
+                    comps = {"Link"}
                 else:
                     comps = {"Generator", "Link", "StorageUnit", "Store"}
                 for c in n.iterate_components(comps):
@@ -3346,6 +3533,9 @@ def maybe_adjust_costs_and_potentials(n, opts):
                     else:
                         sel = c.df.carrier.str.contains(carrier)
                     c.df.loc[sel, attr] *= factor
+                   
+     
+            
             logger.info(f"changing {attr} for {carrier} by factor {factor}")
 
 
@@ -3619,7 +3809,11 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
     options = snakemake.params.sector
 
     opts = snakemake.wildcards.sector_opts.split("-")
-
+    margoil = 1
+    for o in opts:
+        if "oilmarg" in o:
+            oo = o.split("+")
+            margoil = float(oo[1])
     investment_year = int(snakemake.wildcards.planning_horizons[-4:])
 
     n = pypsa.Network(snakemake.input.network)
@@ -3678,7 +3872,20 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
         options["use_electrolysis_waste_heat"] = False
 
     if "T" in opts:
+        
+        bev_availability = options.get("bev_availability")
+
+        for o in opts:
+
+            if "bevavail" not in o:
+                continue
+            oo = o.split("+")
+
+            bev_availability = float(oo[1])
+
+
         add_land_transport(n, costs)
+        #adjust_land_transport(n)
 
     if "H" in opts:
         add_heat(n, costs)
@@ -3716,13 +3923,16 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
     solver_name = snakemake.config["solving"]["solver"]["name"]
     n = set_temporal_aggregation(n, opts, solver_name)
 
+    if "T" in opts:
+        adjust_land_transport(n)
+
     limit_type = "config"
     limit = get(snakemake.params.co2_budget, investment_year)
     for o in opts:
         if "cb" not in o:
             continue
         limit_type = "carbon budget"
-        fn = "results/" + snakemake.params.RDIR + "/csvs/carbon_budget_distribution.csv"
+        fn = "results/" + snakemake.params.RDIR + "/csvs/carbon_budget_distribution" + o +".csv"
         if not os.path.exists(fn):
             emissions_scope = snakemake.params.emissions_scope
             report_year = snakemake.params.eurostat_report_year
@@ -3733,7 +3943,7 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
                 fn,
                 emissions_scope,
                 report_year,
-                input_co2,
+                #input_co2,
             )
         co2_cap = pd.read_csv(fn, index_col=0).squeeze()
         limit = co2_cap.loc[investment_year]
@@ -3759,7 +3969,7 @@ def lossy_bidirectional_links(n, carrier, efficiencies={}):
         insert_electricity_distribution_grid(n, costs)
 
     maybe_adjust_costs_and_potentials(n, opts)
-
+    print(n.links.capital_cost[n.links.carrier.str.contains("land transport oil")])
     if options["gas_distribution_grid"]:
         insert_gas_distribution_costs(n, costs)
 
diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py
index f88d10d48..158984776 100644
--- a/scripts/simplify_network.py
+++ b/scripts/simplify_network.py
@@ -251,7 +251,6 @@ def replace_components(n, c, df, pnl):
     _adjust_capital_costs_using_connection_costs(n, connection_costs_to_bus, output)
 
     generator_strategies = aggregation_strategies["generators"]
-
     carriers = set(n.generators.carrier) - set(exclude_carriers)
     generators, generators_pnl = aggregateoneport(
         n,
@@ -262,11 +261,10 @@ def replace_components(n, c, df, pnl):
     )
 
     replace_components(n, "Generator", generators, generators_pnl)
-
     for one_port in aggregate_one_ports:
         df, pnl = aggregateoneport(n, busmap, component=one_port)
         replace_components(n, one_port, df, pnl)
-
+    
     buses_to_del = n.buses.index.difference(busmap)
     n.mremove("Bus", buses_to_del)
     for c in n.branch_components:
@@ -300,7 +298,7 @@ def simplify_links(
     labels = pd.Series(labels, n.buses.index)
 
     G = n.graph()
-
+    
     def split_links(nodes):
         nodes = frozenset(nodes)
 
@@ -339,7 +337,7 @@ def split_links(nodes):
     connection_costs_to_bus = pd.DataFrame(
         0.0, index=n.buses.index, columns=list(connection_costs_per_link)
     )
-
+    
     for lbl in labels.value_counts().loc[lambda s: s > 2].index:
         for b, buses, links in split_links(labels.index[labels == lbl]):
             if len(buses) <= 2:
@@ -401,7 +399,7 @@ def split_links(nodes):
             # n.add("Link", **params)
 
     logger.debug("Collecting all components using the busmap")
-
+    
     _aggregate_and_move_components(
         n,
         busmap,
@@ -410,6 +408,7 @@ def split_links(nodes):
         aggregation_strategies=aggregation_strategies,
         exclude_carriers=exclude_carriers,
     )
+    
     return n, busmap
 
 
@@ -529,25 +528,26 @@ def cluster(
 
         snakemake = mock_snakemake("simplify_network", simpl="")
     configure_logging(snakemake)
-
+    
     params = snakemake.params
     solver_name = snakemake.config["solving"]["solver"]["name"]
 
     n = pypsa.Network(snakemake.input.network)
+    
     Nyears = n.snapshot_weightings.objective.sum() / 8760
 
     # remove integer outputs for compatibility with PyPSA v0.26.0
     n.generators.drop("n_mod", axis=1, inplace=True, errors="ignore")
-
+    
     n, trafo_map = simplify_network_to_380(n)
-
+    
     technology_costs = load_costs(
         snakemake.input.tech_costs,
         params.costs,
         params.max_hours,
         Nyears,
     )
-
+    
     n, simplify_links_map = simplify_links(
         n,
         technology_costs,
@@ -560,7 +560,7 @@ def cluster(
     )
 
     busmaps = [trafo_map, simplify_links_map]
-
+    
     if params.simplify_network["remove_stubs"]:
         n, stub_map = remove_stubs(
             n,
@@ -572,11 +572,11 @@ def cluster(
             aggregation_strategies=params.aggregation_strategies,
         )
         busmaps.append(stub_map)
-
+    
     if params.simplify_network["to_substations"]:
         n, substation_map = aggregate_to_substations(n, params.aggregation_strategies)
         busmaps.append(substation_map)
-
+    
     # treatment of outliers (nodes without a profile for considered carrier):
     # all nodes that have no profile of the given carrier are being aggregated to closest neighbor
     if params.simplify_network["algorithm"] == "hac":
@@ -592,7 +592,7 @@ def cluster(
                 n, params.aggregation_strategies, buses_i
             )
             busmaps.append(busmap_hac)
-
+    
     if snakemake.wildcards.simpl:
         n, cluster_map = cluster(
             n,
@@ -620,7 +620,7 @@ def cluster(
     n.lines.drop(remove, axis=1, errors="ignore", inplace=True)
 
     update_p_nom_max(n)
-
+    
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output.network)
 
@@ -628,3 +628,4 @@ def cluster(
     busmap_s.to_csv(snakemake.output.busmap)
 
     cluster_regions(busmaps, snakemake.input, snakemake.output)
+    
diff --git a/scripts/solve_network.py b/scripts/solve_network.py
index 1c0f6505c..bcd6a6d74 100644
--- a/scripts/solve_network.py
+++ b/scripts/solve_network.py
@@ -44,6 +44,14 @@
 pypsa.pf.logger.setLevel(logging.WARNING)
 from pypsa.descriptors import get_switchable_as_dense as get_as_dense
 
+from pathlib import Path
+import os
+
+tmpdir = '/scratch/' + os.environ['SLURM_JOB_ID']
+
+if tmpdir is not None:
+    Path(tmpdir).mkdir(parents=True, exist_ok=True)
+
 
 def add_land_use_constraint(n, planning_horizons, config):
     if "m" in snakemake.wildcards.clusters:
@@ -794,6 +802,216 @@ def add_pipe_retrofit_constraint(n):
 
     n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit")
 
+def add_battery_constraints(n):
+    """
+    Add constraint ensuring that charger = discharger, i.e.
+    1 * charger_size - efficiency * discharger_size = 0
+    """
+    if not n.links.p_nom_extendable.any():
+        return
+
+    discharger_bool = n.links.index.str.contains("battery discharger")
+    charger_bool = n.links.index.str.contains("battery charger")
+
+    dischargers_ext = n.links[discharger_bool].query("p_nom_extendable").index
+    chargers_ext = n.links[charger_bool].query("p_nom_extendable").index
+    eff = n.links.efficiency[dischargers_ext].values
+    lhs = (
+        n.model["Link-p_nom"].loc[chargers_ext]
+        - n.model["Link-p_nom"].loc[dischargers_ext] * eff
+    )
+
+    n.model.add_constraints(lhs == 0, name="Link-charger_ratio")
+
+
+def add_chp_constraints(n):
+    electric = (
+        n.links.index.str.contains("urban central")
+        & n.links.index.str.contains("CHP")
+        & n.links.index.str.contains("electric")
+    )
+    heat = (
+        n.links.index.str.contains("urban central")
+        & n.links.index.str.contains("CHP")
+        & n.links.index.str.contains("heat")
+    )
+
+    electric_ext = n.links[electric].query("p_nom_extendable").index
+    heat_ext = n.links[heat].query("p_nom_extendable").index
+
+    electric_fix = n.links[electric].query("~p_nom_extendable").index
+    heat_fix = n.links[heat].query("~p_nom_extendable").index
+
+    p = n.model["Link-p"]  # dimension: [time, link]
+
+    # output ratio between heat and electricity and top_iso_fuel_line for extendable
+    if not electric_ext.empty:
+        p_nom = n.model["Link-p_nom"]
+
+        lhs = (
+            p_nom.loc[electric_ext]
+            * (n.links.p_nom_ratio * n.links.efficiency)[electric_ext].values
+            - p_nom.loc[heat_ext] * n.links.efficiency[heat_ext].values
+        )
+        n.model.add_constraints(lhs == 0, name="chplink-fix_p_nom_ratio")
+
+        rename = {"Link-ext": "Link"}
+        lhs = (
+            p.loc[:, electric_ext]
+            + p.loc[:, heat_ext]
+            - p_nom.rename(rename).loc[electric_ext]
+        )
+        n.model.add_constraints(lhs <= 0, name="chplink-top_iso_fuel_line_ext")
+
+    # top_iso_fuel_line for fixed
+    if not electric_fix.empty:
+        lhs = p.loc[:, electric_fix] + p.loc[:, heat_fix]
+        rhs = n.links.p_nom[electric_fix]
+        n.model.add_constraints(lhs <= rhs, name="chplink-top_iso_fuel_line_fix")
+
+    # back-pressure
+    if not electric.empty:
+        lhs = (
+            p.loc[:, heat] * (n.links.efficiency[heat] * n.links.c_b[electric].values)
+            - p.loc[:, electric] * n.links.efficiency[electric]
+        )
+        n.model.add_constraints(lhs <= rhs, name="chplink-backpressure")
+
+
+def add_pipe_retrofit_constraint(n):
+    """
+    Add constraint for retrofitting existing CH4 pipelines to H2 pipelines.if
+    """
+    gas_pipes_i = n.links.query("carrier == 'gas pipeline' and p_nom_extendable").index
+    h2_retrofitted_i = n.links.query(
+        "carrier == 'H2 pipeline retrofitted' and p_nom_extendable"
+    ).index
+
+    if h2_retrofitted_i.empty or gas_pipes_i.empty:
+        return
+
+    p_nom = n.model["Link-p_nom"]
+
+    CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"]
+    lhs = p_nom.loc[gas_pipes_i] + CH4_per_H2 * p_nom.loc[h2_retrofitted_i]
+    rhs = n.links.p_nom[gas_pipes_i].rename_axis("Link-ext")
+
+    n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit")
+
+
+def add_v2g_constraint(n, bev_availability):
+    bev_charger = n.links.carrier.str.contains('BEV charger')
+    bev_charger_ext = n.links[bev_charger].query("p_nom_extendable").index
+    lhs1 = n.model["Link-p_nom"].loc[bev_charger_ext]*bev_availability
+    v2g = n.links.carrier.str.contains('V2G')
+    v2g_ext = n.links[v2g].query("p_nom_extendable").index
+    lhs2 = n.model["Link-p_nom"].loc[v2g_ext]
+    n.model.add_constraints(lhs1==lhs2, name="constraint_v2g")
+
+
+def add_test_storage_constraint(n):
+    ev_store = n.stores.carrier.str.contains("urban central")
+    ev_store_ext = n.stores[ev_store].query("e_nom_extendable").index[12]
+    print(ev_store_ext)
+    lhs2 = n.model.variables['Store-e_nom'].loc[ev_store_ext]
+
+    ev = n.links.carrier.str.contains('land transport EV')
+    ev_ext = n.links[ev].query("p_nom_extendable").index
+    lhs1 = n.model["Link-p_nom"].loc[ev_ext]/n.config['sector']['EV_consumption_1car']
+    n.model.add_constraints(lhs2==1000, name="test_store")
+
+#constrains EV battery in relation to parking cars based on land transport EV installed
+def add_EV_storage_constraint(n, bev_availability, ev_ext, ev_store_ext, i):
+
+
+    lhs1 = n.model["Link-p_nom"].loc[ev_ext]
+
+    lhs2 = n.model.variables['Store-e_nom'].loc[ev_store_ext]
+
+    factor = bev_availability*n.config["sector"]["bev_avail_max"]/(1-n.config["sector"]["bev_avail_max"])*(n.config['sector']['bev_energy'])
+
+    #factor = n.config["sector"]["bev_availability"]*n.config["sector"]["bev_avail_max"]/(1-n.config["sector"]["bev_avail_max"])*(n.config['sector']['bev_energy'])
+    rhs = lhs1*factor/n.config['sector']['EV_consumption_1car']
+
+    n.model.add_constraints(lhs2==rhs, name="constraint_EV_storage"+str(i))
+
+
+def add_BEV_charger_constraint(n, i):
+
+    bev = n.links.carrier.str.contains('BEV charger') 
+    bev_links = n.links[bev].index[i]
+    bus = n.links.bus0[bev_links]
+    indexBEV = n.links[(n.links.bus0==bus)].index.intersection(n.links[bev].index)
+    existing_link_BEV = n.links.p_nom[indexBEV].sum()
+    bev_charger = n.links.carrier.str.contains('BEV charger')
+    bev_charger_ext = n.links[bev_charger].query("p_nom_extendable").index[i]
+    lhs = (n.model["Link-p_nom"].loc[bev_charger_ext]+existing_link_BEV)/n.config['sector']['bev_charge_rate']
+
+
+
+    ev = n.links.carrier.str.contains('land transport EV') 
+    ev_links = n.links[ev].index[i]
+    bus = n.links.bus0[ev_links]
+    indexEV = n.links[(n.links.bus0==bus)].index.intersection(n.links[ev].index)
+    existing_link_EV = n.links.p_nom[indexEV].sum()
+    ev = n.links.carrier.str.contains('land transport EV')
+    ev_ext = n.links[ev].query("p_nom_extendable").index[i]
+    rhs2 = (n.model["Link-p_nom"].loc[ev_ext]+existing_link_EV)/n.config['sector']['EV_consumption_1car']
+
+    #factor = n.config["sector"]["bev_avail_max"]/(1-n.config["sector"]["bev_avail_max"])*bev_availability
+    factor  = n.config["sector"]["bev_charge_avail"]
+
+    #factor = n.config["sector"]["bev_availability"]*n.config["sector"]["bev_avail_max"]/(1-n.config["sector"]["bev_avail_max"])*(n.config['sector']['bev_energy'])
+    rhs = rhs2*factor
+
+    n.model.add_constraints(lhs==rhs, name="constraint_BEV_charger"+str(i))
+
+#important for myopic: makes sure that the overall capacity of EV batteries is at least as much as land transport used
+def add_EV_storage_constraint2(n, bev_availability, ev_ext, ev_store_ext, i):
+
+
+    rhs1 = n.model["Link-p_nom"].loc[ev_ext]
+
+    ev = n.links.carrier.str.contains('land transport EV') 
+    ev_links = n.links[ev].index[i]
+    bus = n.links.bus0[ev_links]
+    indexEV = n.links[(n.links.bus0==bus)].index.intersection(n.links[ev].index)
+    existing_link = n.links.p_nom[indexEV].sum()
+
+    lhs2 = n.model.variables['Store-e_nom'].loc[ev_store_ext]
+
+    store = n.stores.carrier.str.contains('Li ion|EV battery storage') 
+    ev_store = n.stores[store].index[i]
+    bus = n.stores.bus[ev_store]
+    indexstore = n.stores[(n.stores.bus==bus)].index
+    existing_store = n.stores.e_nom[indexstore].sum()
+
+    lhs = (lhs2 + existing_store)/(n.config['sector']['bev_energy'])
+
+    #factor = n.config["sector"]["bev_avail_max"]/(1-n.config["sector"]["bev_avail_max"])
+    factor = 5
+    rhs = (rhs1+existing_link)/n.config['sector']['EV_consumption_1car']*bev_availability*factor
+
+    n.model.add_constraints(lhs==rhs, name="constraint_EV_storage2_"+str(i))
+
+def add_EV_number_constraint(n):
+    bev_charger = n.links.carrier.str.contains('BEV charger')
+    max_bev = n.links.p_nom_max[n.links[bev_charger].query("p_nom_extendable").index]/n.config['sector']['bev_charge_rate']
+
+    bev_charger_ext = n.links[bev_charger].query("p_nom_extendable").index
+    lhs1 = n.model["Link-p_nom"].loc[bev_charger_ext]/n.config['sector']['bev_charge_rate']
+
+    v2g = n.links.carrier.str.contains('V2G')
+    v2g_ext = n.links[v2g].query("p_nom_extendable").index
+    lhs1 = n.model["Link-p_nom"].loc[v2g_ext]/n.config['sector']['bev_charge_rate']
+    
+    ev = n.links.carrier.str.contains('land transport EV')
+    ev_ext = n.links[ev].query("p_nom_extendable").index
+    lhs2 = n.model["Link-p_nom"].loc[ev_ext]/n.config['sector']['EV_consumption_1car']
+
+    rhs = lhs2*1E6
+
+    n.model.add_constraints(lhs1<=rhs, name="constraint_EV_number3")
 
 def extra_functionality(n, snapshots):
     """
@@ -837,7 +1055,47 @@ def extra_functionality(n, snapshots):
         add_carbon_constraint(n, snapshots)
         add_carbon_budget_constraint(n, snapshots)
         add_retrofit_gas_boiler_constraint(n, snapshots)
-
+    current_horizon = int(snakemake.wildcards.planning_horizons)
+
+
+    if config['sector']["land_transport_electric_share"][current_horizon] is None and any(n.links.loc[(n.links[n.links.carrier.str.contains("BEV charger")].index),'p_nom_extendable']):   
+        print(any(n.links.loc[(n.links[n.links.carrier.str.contains("BEV charger")].index),'p_nom_extendable']))
+        bev_availability = n.config["sector"]["bev_availability"]
+        for o in opts:
+
+            if "bevavail" not in o:
+                continue
+            oo = o.split("+")
+
+            bev_availability = float(oo[1])
+        BEV_charger_ext = n.links[n.links.carrier.str.contains("BEV charger")].query("p_nom_extendable").index
+        i = 0
+        while i < len(BEV_charger_ext):
+            BEV_charger_max = n.links.p_nom_max[BEV_charger_ext[i]]
+            if BEV_charger_max > 0.0001:
+                add_BEV_charger_constraint(n, i)
+            i = i+1
+        print(bev_availability)
+        if config['sector']["bev_dsm"]:
+
+          ev_store = n.stores.carrier.str.contains('Li ion|EV battery storage') 
+          ev_store_ext = n.stores[ev_store].query("e_nom_extendable").index
+          ev = n.links.carrier.str.contains('land transport EV') 
+          ev_ext = n.links[ev].query("p_nom_extendable").index
+          i=0
+          while i < len(ev_store_ext):
+               print(n.stores.e_nom_max[ev_store_ext[i]])
+               if n.stores.e_nom_max[ev_store_ext[i]] > 0.0001:
+                print(n.stores.e_nom_max[ev_store_ext[i]])
+                add_EV_storage_constraint2(n, bev_availability, ev_ext[i], ev_store_ext[i], i)
+                #add_EV_storage_constraint2(n, bev_availability, ev_ext[i], ev_store_ext[i], i)
+               i = i+1
+
+
+          #add_EV_storage_constraint(n, bev_availability)
+        if config['sector']["v2g"] and sum(n.links.p_nom_max[n.links[n.links.carrier.str.contains("V2G")].index]):
+          #add_EV_number_constraint(n)
+          add_v2g_constraint(n, bev_availability) 
     if snakemake.params.custom_extra_functionality:
         source_path = snakemake.params.custom_extra_functionality
         assert os.path.exists(source_path), f"{source_path} does not exist"
@@ -863,6 +1121,7 @@ def solve_network(n, config, solving, opts="", **kwargs):
         "linearized_unit_commitment", False
     )
     kwargs["assign_all_duals"] = cf_solving.get("assign_all_duals", False)
+    #kwargs["solver_dir"] = tmpdir
 
     if kwargs["solver_name"] == "gurobi":
         logging.getLogger("gurobipy").setLevel(logging.CRITICAL)
@@ -889,16 +1148,17 @@ def solve_network(n, config, solving, opts="", **kwargs):
         kwargs["min_iterations"] = (cf_solving.get("min_iterations", 4),)
         kwargs["max_iterations"] = (cf_solving.get("max_iterations", 6),)
         status, condition = n.optimize.optimize_transmission_expansion_iteratively(
-            **kwargs
+            **kwargs,
         )
-
+    print(status, condition)
     if status != "ok" and not rolling_horizon:
         logger.warning(
             f"Solving status '{status}' with termination condition '{condition}'"
         )
     if "infeasible" in condition:
+    #if "warning" in status:
         labels = n.model.compute_infeasibilities()
-        logger.info("Labels:\n" + labels)
+        logger.info("Labels:\n{labels}") # + ' '.join([str(elem) for i, elem in enumerate(labels)])) #labels)
         n.model.print_infeasibilities()
         raise RuntimeError("Solving status 'infeasible'")