diff --git a/cases.csv b/cases.csv
index 4840a719..ece24ca2 100644
--- a/cases.csv
+++ b/cases.csv
@@ -33,6 +33,11 @@ coalscen,"Coal price scenario (inputs\fuelprices\{coalscen}.csv). The options fo
 demandscen,Demand scenario (see inputs\load),N/A,AEO_2025_reference,
 dr_shedscen,DR shed scenario,demo_data_January_2025,demo_data_January_2025,
 evmcscen,EVMC scenario for load and supply curve and flexibility,Baseline,Baseline,
+GSw_FINITO_Link,Switch to enable linkage with FINITO,0; 1,0,
+finito_dir,Directory with cloned FINITO repository,N/A,/projects/finito/molmezt/FINITO,
+finito_cases_file,"extension of the FINITO cases file including all scenarios to be run, e.g., 'linked' for cases_linked.csv",N/A,linked,
+finito_case,Case name for FINITO in {finito_dir}/cases_linked.csv; use 'same' if the scenario name in cases_linked.csv matches the scenario name in the cases file,N/A,BAU,
+first_year_finito,First year when the FINITO module is enabled,N/A,2018,
 geodiscov,Annual discover rates for new geothermal sites,BAU; TI,BAU,
 geohydrosupplycurve,Geohydro Supply Curve,ATB_2023,ATB_2023,
 egssupplycurve,EGS Supply Curve,ATB_2023; reV,reV,
diff --git a/cases_finito.csv b/cases_finito.csv
new file mode 100644
index 00000000..6e6c732e
--- /dev/null
+++ b/cases_finito.csv
@@ -0,0 +1,44 @@
+,Default Value,ReEDSonly,BAUreg,Reference,Reference-AEO,Reference-HOG,TradeReference,TradeReference-AEO,Trade-HOG,Reference-LowBio,Reference-HighBio,Tax200,Tax200-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level
+ignore,1,1,0,1,1,1,1,1,1,,,,,,,,,,,,,
+GSw_FINITO_Link,1,0,,,,,,,,,,,,,,,,,,,,
+finito_dir,/path/to/FINITO/directory,0,,,,,,,,,,,,,,,,,,,,
+finito_cases_file,linked,,,,,,,,,,,,,,,,,,,,,
+finito_case,BAU,,BAUreg,,BAU-AEO,AEO-HOG,TradeReEDS,TradeReEDS-AEO,TradeReEDS-HOG,BAU-LowBio,BAU-HighBio,,BAU-AEO,NetZero,NetZero-LowBio,NetZero-HighBio,NetZero-SlackLim1500,NetZero-SlackLim1500-CO2Stor-Mid,NetZero-SlackLim1500-CO2Stor-Low,NetZero-SlackLim500-CO2Stor-Low,NetZero-SlackCost2500-CO2Stor-Low,aggreg_54_level
+coalscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,,
+ngscen,AEO_2025_reference,,,,,AEO_2025_HOG,,,AEO_2025_HOG,,,,,,,,,,,,,
+uraniumscen,AEO_2025_reference,,,,,,,,,,,,,,,,,,,,,
+GSw_Region,country/USA,,st/CO,,,,,,,,,,,,,,,,,,,
+GSw_ZoneSet,z54,,,,,,,,,,,,,,,,,,,,,z54
+GSw_GasCurve,1,,2,,,,,,,,,,,,,,,,,,,
+GSw_gopt,2,,,,,,,,,,,,,,,,,,,,,
+GSw_NewValCapShrink,1,,,,,,,,,,,,,,,,,,,,,
+endyear,2050,,2035,,,,,,,,,,,,,,,,,,,
+yearset,2010_2015_2020..2050..5,,,,,,,,,,,,,,,,,,,,,
+first_year_finito,2018,,,,,,,,,,,,,,,,,,,,,
+incentives_suffix,obbba,,,,,,,,,,,,,,,,,,,,,
+GSw_Upgrades,0,,,,,,,,,,,,,,,,,,,,,
+GSw_AnnualCap,0,,,,,,,,,,,,,,,,,,,,,
+GSw_AnnualCapScen,start2023_100pct2035,,,,,,,,,,,,,,,,,,,,,
+GSw_Clean_Air_Act,0,,,,,,,,,,,,,,,,,,,,,
+GSw_CSP,0,,,,,,,,,,,,,,,,,,,,,
+GSw_Geothermal,0,,,,,,,,,,,,,,,,,,,,,
+GSw_H2,2,,,,,,,,,,,,,,,,,,,,,
+GSw_H2_PTC,0,,,,,,,,,,,,,,,,,,,,,
+GSw_H2_Demand_Case,BAU,,,,,,,,,,,,,,,,,,,,,
+GSw_H2_SMR,1,,,,,,,,,,,,,,,,,,,,,
+GSw_H2_Transport,0,,,,,,,,,,,,,,,,,,,,,
+GSw_H2Combustion,1,,,,,,,,,,,,,,,,,,,,,
+GSw_H2Combustionupgrade,1,,,,,,,,,,,,,,,,,,,,,
+GSw_H2CombinedCycle,1,,,,,,,,,,,,,,,,,,,,,
+GSw_DAC,1,,,,,,,,,,,,,,,,,,,,,
+GSw_BECCS,1,,,,,,,,,,,,,,,,,,,,,
+GSw_TCPhaseout,1,,,,,,,,,,,,,,,,,,,,,
+GSw_CO2_Detail,1,,,,,,,,,,,,,,,,,,,,,
+numhintage,4,,,,,,,,,,,,,,,,,,,,,
+GSw_CarbTax,0,,,,,,,,,,,1,1,,,,,,,,,
+GSw_CarbTaxOption,default,,,,,,,,,,,t200_2050,t200_2050,,,,,,,,,
+GSw_PRM_CapCredit,0,,,,,,,,,,,,,,,,,,,,,
+GSw_PRM_StressIterateMax,0,,,,,,,,,,,,,,,,,,,,,
+GSw_HourlyNumClusters,9,,,,,,,,,,,,,,,,,,,,,
+GSw_HourlyChunkLengthRep,4,,,,,,,,,,,,,,,,,,,,,
+GSw_HourlyChunkLengthStress,4,,,,,,,,,,,,,,,,,,,,,
diff --git a/helpers/restart_runs.py b/helpers/restart_runs.py
index f1132e2e..d353d2ae 100644
--- a/helpers/restart_runs.py
+++ b/helpers/restart_runs.py
@@ -4,6 +4,7 @@
 import shutil
 import subprocess
 import argparse
+import re
 from glob import glob
 from pathlib import Path
 sys.path.append(str(Path(__file__).parent.parent))
@@ -95,6 +96,14 @@
 
     #%% Copy additional files if desired
     for f in more_copyfiles:
+        if f.lower().startswith('finito'):
+            # if file starts with 'finito' append the finito directory
+            sw = reeds.io.get_switches(case)
+            if int(sw.GSw_FINITO_Link):
+                f_copy = re.sub("^finito/", "" , f, flags=re.IGNORECASE)
+                shutil.copy(os.path.join(sw.finito_dir,f_copy), os.path.join(case,f))
+        else:
+            shutil.copy(os.path.join(reeds_path,f), os.path.join(case,f))
         shutil.copy(f, Path(case, f))
     if copy_reeds:
         shutil.copytree(Path(reeds_path, 'reeds'), Path(case, 'reeds'), dirs_exist_ok=True)
diff --git a/reeds/core/setup/a_createmodel.gms b/reeds/core/setup/a_createmodel.gms
index 9250f30e..19db6037 100644
--- a/reeds/core/setup/a_createmodel.gms
+++ b/reeds/core/setup/a_createmodel.gms
@@ -4,7 +4,18 @@ $setglobal ds /
 $endif.unix
 
 $include reeds%ds%core%ds%setup%ds%b_inputs.gms
+
+$ifthene.finito_link %GSw_FINITO_Link% == 1
+$include finito%ds%model%ds%finito_input.gms
+$include finito%ds%model%ds%finito_vars_eqs.gms
+$endif.finito_link
+
 $include reeds%ds%core%ds%setup%ds%c_model.gms
+
+$ifthene.finito_link %GSw_FINITO_Link% == 1
+$include finito%ds%model%ds%finito_model.gms
+$endif.finito_link
+
 $include reeds%ds%core%ds%setup%ds%d_objective.gms
 $include reeds%ds%core%ds%setup%ds%d_mga.gms
 $include reeds%ds%core%ds%setup%ds%e_solveprep.gms
diff --git a/reeds/core/setup/b_inputs.gms b/reeds/core/setup/b_inputs.gms
index 990bc9fc..771cc656 100644
--- a/reeds/core/setup/b_inputs.gms
+++ b/reeds/core/setup/b_inputs.gms
@@ -1018,7 +1018,8 @@ set tmodel(t) "years to include in the model",
     tfix(t) "years to fix variables over when summing over previous years",
     tprev(t,tt) "previous modeled tt from year t",
     stfeas(st) "states to include in the model",
-    tsolved(t) "years that have solved" ;
+    tsolved(t) "years that have solved",
+    tfuel(t) "years that use ReEDS fuel supply curve module (otherwise uses supply curves in FINITO)" ;
 
 *following parameters get re-defined when the solve years have been declared
 parameter mindiff(t) "minimum difference between t and all other tt that are in tmodel(t)" ;
@@ -1031,7 +1032,7 @@ tfix(t) = no ;
 stfeas(st) = no ;
 tprev(t,tt) = no ;
 tsolved(t) = no ;
-
+tfuel(t)=no;
 
 *==============================
 * Year specification
@@ -1046,6 +1047,12 @@ tprev(t,tt)$[tmodel_new(t)$tmodel_new(tt)$(tt.val<t.val)] = yes ;
 mindiff(t)$tmodel_new(t) = smin{tt$tprev(t,tt), t.val-tt.val} ;
 tprev(t,tt)$[tmodel_new(t)$tmodel_new(tt)$(t.val-tt.val<>mindiff(t))] = no ;
 
+* If FINITO linkage is on, remove all modeled years from tfuel after first_year_finito
+* as the FINITO supply curves will be used instead of those in ReEDS;
+* otherwise, all modeled years use ReEDS supply curves and are eligible for tfuel 
+tfuel(t)$[tmodel_new(t)]=yes;
+tfuel(t)$[tmodel_new(t)$Sw_FINITO_Link$(t.val>=%first_year_finito%)] = no ;
+
 * In order to fill all necessary dimensions of upgrade techs parameters, we require
 * Sw_UpgradeYear in ban(i) to be a modeled year and thus we compute as either
 * the GSw_UpgradeYear option or the next modeled years after GSw_UpgradeYear
diff --git a/reeds/core/setup/c_model.gms b/reeds/core/setup/c_model.gms
index 005dca1a..1f9e28a8 100644
--- a/reeds/core/setup/c_model.gms
+++ b/reeds/core/setup/c_model.gms
@@ -396,6 +396,14 @@ eq_loadcon(r,h,t)$tmodel(t)..
 *   (the effect is the same but avoiding the h-indexed OP_LOADSITE reduces solve time)
     + OP_LOADSITE(r,h,t)$[Sw_LoadSiteCF$(Sw_LoadSiteCF<1)$val_loadsite(r)]
     + CAP_LOADSITE(r,t)$[(Sw_LoadSiteCF=1)$val_loadsite(r)]
+
+* [plus] load for industrial and converted fuel facilities (FINITO),
+* including the PRM for stress periods
+* USE_ELE_FINITO is end-use, so divide by (1-distloss) to convert it to busbar
+* [MWh/hr = MW]
+$ifthene.linked_load Sw_FINITO_Link==1
+    + (USE_ELE_FINITO(r,h,t) / (1.0 - distloss))$[t_finito(t)] * (1 + prm(r,t)$h_stress(h))
+$endif.linked_load
 ;
 
 * ---------------------------------------------------------------------------
@@ -2881,7 +2889,7 @@ eq_national_gen(t)$[tmodel(t)$national_gen_frac(t)$Sw_GenMandate]..
 * ---------------------------------------------------------------------------
 
 *gas used from each bin is the sum of all gas used
-eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))]..
+eq_gasused(cendiv,h,t)$[tmodel(t)$tfuel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))]..
 
     sum{gb,GASUSED(cendiv,gb,h,t) }
 
@@ -2898,7 +2906,7 @@ eq_gasused(cendiv,h,t)$[tmodel(t)$((Sw_GasCurve=0) or (Sw_GasCurve=3))]..
 * ---------------------------------------------------------------------------
 
 * gas from each bin needs to less than its capacity
-eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)]..
+eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=0)]..
 
     gaslimit(cendiv,gb,t)
 
@@ -2909,7 +2917,7 @@ eq_gasbinlimit(cendiv,gb,t)$[tmodel(t)$(Sw_GasCurve=0)]..
 
 * ---------------------------------------------------------------------------
 
-eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)]..
+eq_gasbinlimit_nat(gb,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=3)]..
 
    gaslimit_nat(gb,t)
 
@@ -2922,7 +2930,7 @@ eq_gasbinlimit_nat(gb,t)$[tmodel(t)$(Sw_GasCurve=3)]..
 
 * ---------------------------------------------------------------------------
 
-eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)]..
+eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)]..
 
     sum{fuelbin, VGASBINQ_REGIONAL(fuelbin,cendiv,t) }
 
@@ -2935,7 +2943,7 @@ eq_gasaccounting_regional(cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)]..
 
 * ---------------------------------------------------------------------------
 
-eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)]..
+eq_gasaccounting_national(t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)]..
 
     sum{fuelbin,VGASBINQ_NATIONAL(fuelbin,t) }
 
@@ -2948,7 +2956,7 @@ eq_gasaccounting_national(t)$[tmodel(t)$(Sw_GasCurve=1)]..
 
 * ---------------------------------------------------------------------------
 
-eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)]..
+eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)]..
 
     Gasbinwidth_regional(fuelbin,cendiv,t)
 
@@ -2959,7 +2967,7 @@ eq_gasbinlimit_regional(fuelbin,cendiv,t)$[tmodel(t)$(Sw_GasCurve=1)]..
 
 * ---------------------------------------------------------------------------
 
-eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)]..
+eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$tfuel(t)$(Sw_GasCurve=1)]..
 
     Gasbinwidth_national(fuelbin,t)
 
@@ -2973,10 +2981,10 @@ eq_gasbinlimit_national(fuelbin,t)$[tmodel(t)$(Sw_GasCurve=1)]..
 *==============================
 * -- Bioenergy Supply Curve --
 *==============================
+* defer to FINITO representation when models are linked
 
 * ---------------------------------------------------------------------------
-
-eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)]..
+eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)$tfuel(t)]..
 
     sum{bioclass, BIOUSED(bioclass,r,t) }
 
@@ -2995,7 +3003,7 @@ eq_bioused(r,t)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }$tmodel(t)]..
 * ---------------------------------------------------------------------------
 
 * biomass consumption limit is annual
-eq_biousedlimit(bioclass,usda_region,t)$tmodel(t)..
+eq_biousedlimit(bioclass,usda_region,t)$[tmodel(t)$tfuel(t)]..
 
     biosupply(usda_region,bioclass,"cap")
 
@@ -3565,6 +3573,11 @@ eq_h2_demand(p,t)$[(sameas(p,"H2"))$tmodel(t)$(yeart(t)>=h2_demand_start)$(Sw_H2
     + sum{(i,v,r,h)$[valgen(i,v,r,t)$h2_combustion(i)$h_rep(h)],
             GEN(i,v,r,h,t) * hours(h) * h2_combustion_intensity * heat_rate(i,v,r,t)
     }
+
+* hydrogen demand from indusrty (FINITO): FINITO demand [MMBtu/yr] * conversion [metric tons-H2/MMBtu-H2]
+$ifthene.linked_h2_nat Sw_FINITO_Link==1
+    + [sum{(r,h)$h_rep(h), hours(h) * USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu }]$t_finito(t)
+$endif.linked_h2_nat
 ;
 
 * ---------------------------------------------------------------------------
@@ -3596,6 +3609,11 @@ eq_h2_demand_regional(r,h,t)
     + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)],
             GEN(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t)
        }
+
+* regional hydrogen demand for industry (FINITO)
+$ifthene.linked_h2_reg Sw_FINITO_Link==1
+    + [ USE_H2_FINITO(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t)
+$endif.linked_h2_reg
 ;
 
 * ---------------------------------------------------------------------------
@@ -3838,6 +3856,11 @@ eq_co2_capture(r,h,t)
 
 * capture from DAC
     + sum{(i,v)$[dac(i)$valcap(i,v,r,t)$i_p(i,"DAC")], PRODUCE("DAC",i,v,r,h,t) }$Sw_DAC
+* (ReEDS-FINITO) capture from industry [metric_tons-CO2/hr]: 
+* calculation: hours_per_year [yrs/hr] * capture [scaled_metric_tons-CO2/yr] / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]
+$ifthene.linked_co2_capture Sw_FINITO_Link==1
+    + [CAPTURE_CO2EM(r,h,t) / co2_scale]$[t_finito(t)$h_rep(h)]
+$endif.linked_co2_capture
 ;
 
 * ---------------------------------------------------------------------------
@@ -3867,6 +3890,11 @@ eq_co2_spurline_caplimit(r,cs,h,t)$[Sw_CO2_Detail$r_cs(r,cs)$tmodel(t)$(yeart(t)
     =g=
 
     CO2_STORED(r,cs,h,t)
+* (ReEDS-FINITO) extraction of CO2 [metric_tons-CO2/hr] 
+* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] 
+$ifthene.linked_co2_spurline_caplimit Sw_FINITO_Link==1
+    + [ (1 / co2_scale) * EXTRACT_CO2_CS(cs,r,h,t) ]$[h_rep(h)$t_finito(t)]
+$endif.linked_co2_spurline_caplimit
 ;
 
 * ---------------------------------------------------------------------------
@@ -3886,6 +3914,18 @@ eq_co2_sink(r,h,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)]..
 
 * net trade
     + sum{rr$co2_routes(r,rr), CO2_FLOW(rr,r,h,t) - CO2_FLOW(r,rr,h,t) }
+
+* (ReEDS-FINITO) extraction - use of CO2 [metric_tons-CO2/hr] 
+*calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * extraction/use [scaled_metric_tons-CO2/yr]
+$ifthene.linked_co2_sink Sw_FINITO_Link==1
+    + (1 / co2_scale) * [
+*       extraction from all cs sites in r 
+        + sum{cs$[csfeas(cs)$r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t)}
+*       total use of CO2, equivalent to supply 
+        - USE_CO2(r,h,t) 
+    ]$[h_rep(h)$t_finito(t)]
+$endif.linked_co2_sink
+
 ;
 
 * ---------------------------------------------------------------------------
@@ -3899,6 +3939,14 @@ eq_co2_injection_limit(cs,h,t)$[Sw_CO2_Detail$tmodel(t)$(yeart(t)>=co2_detail_st
 
 * must exceed metric tons per hour entering storage
     sum{r$r_cs(r,cs), CO2_STORED(r,cs,h,t) }
+
+
+
+* (ReEDS-FINITO) extraction of CO2 for use [metric_tons-CO2/hr] 
+* calculation: hours_per_year [yrs/hr] * (1 / co2_scale [scaled_metric_tons-CO2/metric_tons-CO2]) * use [scaled_metric_tons-CO2/yr] 
+$ifthene.linked_co2_injection_limit Sw_FINITO_Link==1
+    + (1 / co2_scale) * sum{r$[r_cs(r,cs)], EXTRACT_CO2_CS(cs,r,h,t) }$[h_rep(h)$t_finito(t)]
+$endif.linked_co2_injection_limit
 ;
 
 * ---------------------------------------------------------------------------
@@ -3915,6 +3963,14 @@ eq_co2_cumul_limit(cs,t)$[tmodel(t)$Sw_CO2_Detail$(yeart(t)>=co2_detail_startyr)
         $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr)
         $r_cs(r,cs)$h_rep(h)],
         yearweight(tt) * hours(h) * CO2_STORED(r,cs,h,tt) }
+
+* (ReEDS-FINITO) cumulative amount extracted over time
+$ifthene.linked_co2_storage_cumul_limit Sw_FINITO_Link==1
+    - sum{(r,h,tt)
+        $[(yeart(tt)<=yeart(t))$(tmodel(tt) or tfix(tt))$(yeart(tt)>=co2_detail_startyr)
+        $r_cs(r,cs)$(tfuel(tt))],
+        yearweight(tt) * hours(h) * EXTRACT_CO2_CS(cs,r,h,tt) }$[t_finito(t)]
+$endif.linked_co2_storage_cumul_limit 
 ;
 * ---------------------------------------------------------------------------
 
diff --git a/reeds/core/setup/d_objective.gms b/reeds/core/setup/d_objective.gms
index fae5e5ca..af8672bf 100644
--- a/reeds/core/setup/d_objective.gms
+++ b/reeds/core/setup/d_objective.gms
@@ -19,7 +19,17 @@ Variable    Z        "--$-- total cost of operations and investment, scale varie
 
 * objective function is the sum over modeled years of the investment
 * and operations components
-eq_ObjFn.. Z =e= cost_scale * sum{t$tmodel(t), Z_inv(t) + Z_op(t) } ;
+eq_ObjFn.. Z =e= cost_scale * (
+* electricity and H2 costs
+                    sum{t$tmodel(t), Z_inv(t) + Z_op(t) } 
+* economy-wide costs from FINITO: deflate from $2018 to $2004 
+* and remove any FINITO scaling
+$ifthene.linked_objective Sw_FINITO_Link==1
+                    + deflator('2018') / obj_scale 
+                     * sum{t$[t_finito(t)], Z_finito(t) }
+$endif.linked_objective 
+     )
+;
 
 *=======================================================
 * -- Investment component of the objective function --
@@ -221,7 +231,7 @@ eq_Objfn_op(t)$tmodel(t)..
 * plus cost of H2 fuel when using fixed price (Sw_H2=0) or during stress periods.
 * When using endogenous H2 price (Sw_H2=1 or Sw_H2=2), H2 fuel cost is captured elsewhere
 * via the capex + opex costs of H2 production and its associated electricity demand.
-              + sum{(i,v,r,h)$[valgen(i,v,r,t)$heat_rate(i,v,r,t)
+              + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$heat_rate(i,v,r,t)
                              $(not gas(i))$(not bio(i))$(not cofire(i))
                              $((not h2_combustion(i)) or h2_combustion(i)$[(Sw_H2=0) or h_stress(h)])],
                    hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) }
@@ -232,47 +242,47 @@ eq_Objfn_op(t)$tmodel(t)..
 
 * --cofire coal consumption---
 * cofire bio consumption already accounted for in accounting of BIOUSED
-              + sum{(i,v,r,h)$[valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)],
+              + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$cofire(i)$heat_rate(i,v,r,t)],
                    (1-bio_cofire_perc) * hours(h) * heat_rate(i,v,r,t)
                    * fuel_price("coal-new",r,t) * GEN(i,v,r,h,t) }
 
 * --- cost of natural gas---
 *Sw_GasCurve = 2 (static natural gas prices)
 *first - gas consumed for electricity generation
-              + sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)],
+              + sum{(i,v,r,h)$[tfuel(t)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)$(Sw_GasCurve = 2)],
                    hours(h) * heat_rate(i,v,r,t) * fuel_price(i,r,t) * GEN(i,v,r,h,t) }
 
 *second - gas consumed by gas-powered DAC
-              + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)],
-                   hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$Sw_DAC_Gas
+              + sum{(v,r,h)$[tfuel(t)$valcap("dac_gas",v,r,t)$(Sw_GasCurve = 2)],
+                   hours(h) * dac_gas_cons_rate("dac_gas",v,t) * PRODUCE("DAC","dac_gas",v,r,h,t) }$[Sw_DAC_Gas]
 
 *Sw_GasCurve = 0 (census division supply curves natural gas prices)
-              + sum{(cendiv,gb), sum{h, hours(h) * GASUSED(cendiv,gb,h,t) }
+              + sum{(cendiv,gb)$[tfuel(t)], sum{h, hours(h) * GASUSED(cendiv,gb,h,t) }
                    * gasprice(cendiv,gb,t)
                    }$(Sw_GasCurve = 0)
 
 *Sw_GasCurve = 3 (national supply curve for natural gas prices with census division multipliers)
-              + sum{(h,cendiv,gb), hours(h) * GASUSED(cendiv,gb,h,t)
+              + sum{(h,cendiv,gb)$[tfuel(t)], hours(h) * GASUSED(cendiv,gb,h,t)
                    * gasadder_cd(cendiv,t,h) + gasprice_nat_bin(gb,t)
                    }$(Sw_GasCurve = 3)
 
 *Sw_GasCurve = 1 (national and census division supply curves for natural gas prices)
 *first - anticipated costs of gas consumption given last year's amount
-              + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)],
+              + (sum{(i,r,v,cendiv,h)$[valgen(i,v,r,t)$gas(i)$tfuel(t)],
                    gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) *
                    hours(h) * heat_rate(i,v,r,t) * GEN(i,v,r,h,t) }
 
 *second - adjustments based on changes from last year's consumption at the regional and national level
-              + sum{(fuelbin,cendiv),
+              + sum{(fuelbin,cendiv)$tfuel(t),
                    gasbinp_regional(fuelbin,cendiv,t) * VGASBINQ_REGIONAL(fuelbin,cendiv,t) }
 
-              + sum{(fuelbin),
+              + sum{(fuelbin)$tfuel(t),
                    gasbinp_national(fuelbin,t) * VGASBINQ_NATIONAL(fuelbin,t) }
 
               )$[Sw_GasCurve = 1]
 
 * ---cost of biofuel consumption and biomass transport---
-              + sum{(r,bioclass)$[sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }],
+              + sum{(r,bioclass)$[tfuel(t)$sum{(i,v)$(bio(i) or cofire(i)), valgen(i,v,r,t) }],
                    BIOUSED(bioclass,r,t) *
                    (sum{usda_region$r_usda(r,usda_region), biosupply(usda_region, bioclass, "price") } + bio_transport_cost) }
 
diff --git a/reeds/core/solve/2_temporal_params.gms b/reeds/core/solve/2_temporal_params.gms
index 5c0c7ab5..6e0b2a81 100644
--- a/reeds/core/solve/2_temporal_params.gms
+++ b/reeds/core/solve/2_temporal_params.gms
@@ -302,9 +302,6 @@ can_exports_h(r,h,t)$[hours(h)] = can_exports(r,t) * can_exports_h_frac(h) / hou
 
 $endif.Canada
 
-* zero Canada exports when Canada is not modeled
-can_imports_szn(r,szn,t)$[Sw_Canada=0] = 0 ;
-can_exports_h(r,h,t)$[Sw_Canada=0] = 0 ;
 
 $onempty
 parameter canmexload(r,allh) "load for canadian and mexican regions"
@@ -890,12 +887,20 @@ szn_adj_gas(h)$frac_h_quarter_weights(h,"wint") =
     szn_adj_gas(h) + frac_h_quarter_weights(h,"wint") * szn_adj_gas_winter ;
 
 
+*=============================================
+* ----- ReEDS-FINITO temporal parameters -----
+*=============================================
+* remove focused sectors load from load_exog_static
+$ifthene.linked_finito_temporal_params %GSw_FINITO_Link% ==1
+$include finito/model/linked_finito_temporal_params.gms
+$endif.linked_finito_temporal_params
+
+
 *=============================================
 * -- Round parameters for GAMS --
 *=============================================
 
 avail(i,r,h)$avail(i,r,h) = round(avail(i,r,h),3) ;
-can_imports_szn(r,szn,t)$can_imports_szn(r,szn,t) = round(can_imports_szn(r,szn,t),3) ;
 can_exports_h(r,h,t)$can_exports_h(r,h,t) = round(can_exports_h(r,h,t),3) ;
 frac_h_ccseason_weights(h,ccseason)$frac_h_ccseason_weights(h,ccseason) = round(frac_h_ccseason_weights(h,ccseason),3) ;
 gasadder_cd(cendiv,t,h)$gasadder_cd(cendiv,t,h) = round(gasadder_cd(cendiv,t,h),3) ;
diff --git a/reeds/core/solve/3_solve_oneyear.gms b/reeds/core/solve/3_solve_oneyear.gms
index dde80524..d059ff3a 100644
--- a/reeds/core/solve/3_solve_oneyear.gms
+++ b/reeds/core/solve/3_solve_oneyear.gms
@@ -23,6 +23,13 @@ tload(t) = no ;
 tmodel(t) = no ;
 tmodel("%cur_year%") = yes ;
 
+* (ReEDS-FINITO) also reset t_finito
+$ifthene.linked_finito_time %GSw_FINITO_Link% == 1
+t_finito(t) = no ; 
+t_finito(t)$[tmodel(t)$(t.val>=first_year_finito)$(t.val<=endyear)] = yes ;
+$endif.linked_finito_time
+
+
 $log 'Solving sequential case for...'
 $log '  Case: %case%'
 $log '  Year: %cur_year%'
@@ -275,11 +282,21 @@ $endif.mga
 
 *** Adjust some parameters based on the solution for this solve year
 $include reeds%ds%core%ds%solve%ds%4_post_solve_adjustments.gms
+$ifthene.linked_finito_shrink %GSw_FINITO_Link% == 1
+* Industrial variable fix
+$include finito/model/finito_post_solve_adjustments.gms
+$endif.linked_finito_shrink
 
 *** Fix decision variables to their optimized levels for this solve year
 tfix("%cur_year%") = yes ;
 $include reeds%ds%core%ds%solve%ds%5_varfix.gms
 
+*** Fix FINITO decision variables to their optimized levels for this solve year
+$ifthene.linked_finito_varfix %GSw_FINITO_Link% == 1
+* Industrial variable fix
+$include finito/model/finito_varfix.gms
+$endif.linked_finito_varfix
+
 *** Dump data used in calculations between solve years
 $include reeds%ds%core%ds%solve%ds%6_data_dump.gms
 
diff --git a/reeds/core/solve/6_data_dump.gms b/reeds/core/solve/6_data_dump.gms
index 5e804c1d..cf534766 100644
--- a/reeds/core/solve/6_data_dump.gms
+++ b/reeds/core/solve/6_data_dump.gms
@@ -8,6 +8,7 @@ This file creates a gdx file with all of the data necessary for the resource ade
     - availability rates (1 - outage rates)
     - transmission capacities and loss rates
     - technology sets
+    - industrial electricity demand from FINITO when running linked model 
 $offtext
 
 $if not set start_year $setglobal start_year %startyear%
@@ -62,6 +63,7 @@ h2_usage_regional(r,allh,t)        "--metric tons-- H2 usage by region"
 inv_cond_filt(i,v,t)               "--set-- vintage-year mapping for investments by technology"
 inv_ivrt(i,v,r,t)                  "--MW-- investments in power generation capacity"
 inv_energy_ivrt(i,v,r,t)           "--MWh-- investments in energy generation capacity"
+load_finito_rt(r,allh,t)           "--MWh-- total load from industrial sectors and converted fuels modeled by FINITO"
 m_cf_filt(i,v,r,allh)              "--fraction-- capacity factor used in the model"
 m_cf_szn_filt(i,v,r,allszn)        "--fraction-- modelled capacity factors filtered for hydro resources to set seasonal energy constraints"
 minloadfrac_filt(r,i,allszn)       "--fraction-- modelled mingen fraction filtered for hydro resources to set mingen constraints"
@@ -161,6 +163,10 @@ h2_usage_regional(r,h,t)$tcur(t) =
         h2_exogenous_demand_regional(r,'h2',h,t)
         + sum{(i,v)$[valgen(i,v,r,t)$h2_combustion(i)],
             GEN.l(i,v,r,h,t) * h2_combustion_intensity * heat_rate(i,v,r,t)}
+*       regional hydrogen demand for industry (FINITO)
+$ifthene.linked_h2_reg_report Sw_FINITO_Link==1
+        + [USE_H2_FINITO.l(r,h,t) * h2_metric_tons_per_mmbtu ]$t_finito(t)  
+$endif.linked_h2_reg_report
     )
 ;
 
@@ -284,15 +290,24 @@ cap_converter_filt(r) = sum{t$tcur(t), CAP_CONVERTER.l(r,t) } ;
 routes_filt(r,rr,trtype) = sum{t$tcur(t), routes(r,rr,trtype,t) } ;
 
 *============================
-* Flexible load data
+* Load data
 *============================
 
+* flexible load
 flex_load(r,h) = sum{(flex_type,t)$tcur(t), load_exog_flex(flex_type,r,h,t) } ;
 
 flex_load_opt(r,h) = sum{(flex_type,t)$tcur(t), FLEX.l(flex_type,r,h,t) } ;
 
 ra_cap_loadsite(r,t)$[Sw_LoadSiteCF$val_loadsite(r)] = CAP_LOADSITE.l(r,t) ;
 
+* FINITO load
+$ifthene.linked_load Sw_FINITO_Link==1
+    load_finito_rt(r,h,t) = USE_ELE_FINITO.l(r,h,t) ;
+$else.linked_load
+    load_finito_rt(r,h,t) = 0 ; 
+$endif.linked_load
+
+
 *============================
 * Extra consumption data
 *============================
@@ -382,6 +397,7 @@ execute_unload 'handoff%ds%reeds_data%ds%reeds_data_%cur_year%.gdx'
     inv_ivrt
     inv_energy_ivrt
     ivt_num
+    load_finito_rt
     m_cf_filt
     m_cf_szn_filt
     maxage
diff --git a/reeds/core/terminus/report.gms b/reeds/core/terminus/report.gms
index bb729367..5a0f5e8d 100644
--- a/reeds/core/terminus/report.gms
+++ b/reeds/core/terminus/report.gms
@@ -107,7 +107,7 @@ rev_cat "categories for renvenue streams" /load, res_marg, oper_res, rps, charge
 
 lcoe_cat "categories for LCOE calculation" /capcost, upgradecost, rsccost, fomcost, vomcost, gen /
 
-loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac /
+loadtype "categories for types of load" / end_use, dist_loss, trans_loss, stor_charge, h2_prod, h2_network, dac, finito /
 
 h2_demand_type / "electricity", "cross-sector"/
 
@@ -495,18 +495,18 @@ ptc_out(i,v,t)$[tmodel_new(t)$ptc_value_scaled(i,v,t)] = ptc_value_scaled(i,v,t)
 * Case 2: the resource of one or more biomass classes ARE exhausted, i.e., BIOUSED.l(bioclass) = biosupply(bioclass)
 *    Marginal Biomass Price = maximum difference between eq_bioused.m and eq_biousedlimit.m(bioclass) across all biomass classes in a region
 
-repbioprice(r,t)$tmodel_new(t) = max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) -
+repbioprice(r,t)$[tmodel_new(t)$tfuel(t)] = max{0, smax{bioclass$BIOUSED.l(bioclass,r,t), eq_bioused.m(r,t) -
                                               sum{usda_region$r_usda(r,usda_region), eq_biousedlimit.m(bioclass,usda_region,t) } } } / pvf_onm(t) ;
 
 * quantity of biomass used (convert from mmBTU to dry tons using biomass energy content)
-bioused_out(bioclass,r,t)$tmodel_new(t) = BIOUSED.l(bioclass,r,t) / bio_energy_content ;
-bioused_usda(bioclass,usda_region,t)$tmodel_new(t) = sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ;
+bioused_out(bioclass,r,t)$[tmodel_new(t)$tfuel(t)] = BIOUSED.l(bioclass,r,t) / bio_energy_content ;
+bioused_usda(bioclass,usda_region,t)$[tmodel_new(t)$tfuel(t)] = sum{r$r_usda(r,usda_region), bioused_out(bioclass,r,t) } ;
 
 * 1e9 converts from MMBtu to Quads
-repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)] =
+repgasquant(cendiv,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 3)$tmodel_new(t)$tfuel(t)] =
     sum{(gb,h), GASUSED.l(cendiv,gb,h,t) * hours(h) } * gas_scale/ 1e9 ;
 
-repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2)$tmodel_new(t)] =
+repgasquant(cendiv,t)$[(Sw_GasCurve = 1 or Sw_GasCurve = 2 or Sw_FINITO_Link = 1)$tmodel_new(t)] =
     ( sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)],
           hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t)}
     + sum{(v,r,h)$[valcap("dac_gas",v,r,t)$r_cendiv(r,cendiv)],
@@ -528,27 +528,52 @@ repgasquant_nat(t)$tmodel_new(t) = sum{cendiv, repgasquant(cendiv,t) } ;
 
 *for reported gasprice (not that used to compute system costs)
 *scale back to $ / mmbtu
-repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)] =
+repgasprice(cendiv,t)$[(Sw_GasCurve = 0)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] =
     smax{gb$[sum{h, GASUSED.l(cendiv,gb,h,t) }], gasprice(cendiv,gb,t) } / gas_scale ;
 
-repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)] =
+repgasprice(cendiv,t)$[(Sw_GasCurve = 2)$tmodel_new(t)$repgasquant(cendiv,t)$tfuel(t)] =
     sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)],
           hours(h)*heat_rate(i,v,r,t)*fuel_price(i,r,t)*GEN.l(i,v,r,h,t)
        } / (repgasquant(cendiv,t) * 1e9) ;
 
-repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)] = sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ;
 
-repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)] =
-              ( sum{(h,cendiv),
-                   gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) *
-                   hours(h) } / sum{h, hours(h) }
+* gas price by timeslice when linked with FINITO (see similar calculation in finito_report.gms) [$2004/MMBtu] 
+* We maybe should be taking weighted averages instead of max across regions and categories (pool vs roi etc.)
+$ifthene.finitogasprice Sw_FINITO_Link == 1
+repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] =
+    ( 1/obj_scale * 1/pvf_onm(t) * deflator('2018') * 
+        [ max{  
+            smax{(eus,usep,r)$[map_roe_ec_usep("NG",eus,usep)$r_cendiv(r,cendiv)], eq_supplydemand_fe_pool.m("NG",eus,usep,r,h,t)},
+            smax{(usep,r)$[map_roi_ec_usep("NG",usep)$r_cendiv(r,cendiv)],eq_supplydemand_fe_roi.m("NG",usep,r,h,t)},
+            smax{(cfp2,cfvin2,fi,r)$[map_ec_ei("NG",fi)$valei_cf(cfp2,cfvin2,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_cf.m(cfp2,cfvin2,fi,r,h,t)},
+            smax{(tech,vin,fi,r)$[map_ec_ei("NG",fi)$valei_ind(tech,vin,fi,r,t)$r_cendiv(r,cendiv)],eq_use_fi_ind.m(tech,vin,fi,r,h,t)}
+            } / hours(h)
+        ] 
+    )
+;
+$else.finitogasprice
+    repgasprice_finito(cendiv,h,t)$[tmodel_new(t)$(not tfuel(t))] = 0 ; 
+$endif.finitogasprice
+
+repgasprice(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] =
+* average price over hours when linked with FINITO
+    sum{h, hours(h) * repgasprice_finito(cendiv,h,t) }
+;
+
+repgasprice_r(r,t)$[(Sw_GasCurve = 0 or Sw_GasCurve = 2)$tmodel_new(t)$tfuel(t)] = 
+  [sum{cendiv$r_cendiv(r,cendiv), repgasprice(cendiv,t) } ];
+
+repgasprice_r(r,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$tfuel(t)] =
+   sum{(h,cendiv),
+       gasmultterm(cendiv,t) * szn_adj_gas(h) * cendiv_weights(r,cendiv) *
+       hours(h) } / sum{h, hours(h) }
 
-              + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) )
+    + smax((fuelbin,cendiv)$[VGASBINQ_REGIONAL.l(fuelbin,cendiv,t)$r_cendiv(r,cendiv)], gasbinp_regional(fuelbin,cendiv,t) )
 
-              + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) )
-              ) ;
+    + smax(fuelbin$VGASBINQ_NATIONAL.l(fuelbin,t), gasbinp_national(fuelbin,t) )
+    ;
 
-repgasprice(cendiv,t)$[(Sw_GasCurve = 1)$tmodel_new(t)$repgasquant(cendiv,t)] =
+repgasprice(cendiv,t)$[((Sw_GasCurve = 1) or (Sw_FINITO_Link = 1))$tmodel_new(t)$repgasquant(cendiv,t)] =
     sum{(i,r)$r_cendiv(r,cendiv), repgasprice_r(r,t) * repgasquant_irt(i,r,t) } / repgasquant(cendiv,t) ;
 
 repgasprice_nat(t)$[tmodel_new(t)$sum{cendiv, repgasquant(cendiv,t) }] =
@@ -560,14 +585,16 @@ repgasprice_nat(t)$[tmodel_new(t)$sum{cendiv, repgasquant(cendiv,t) }] =
 *========================================
 
 gasshare_ba(r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)] =
-                 sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ;
+  sum{i$[valgen_irt(i,r,t)$gas(i)],repgasquant_irt(i,r,t) / repgasquant(cendiv,t) } ;
 
 gasshare_techba(i,r,cendiv,t)$[r_cendiv(r,cendiv)$tmodel_new(t)$repgasquant(cendiv,t)$gas(i)] =
-                 repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ;
+  repgasquant_irt(i,r,t) / repgasquant(cendiv,t) ;
 
-gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ;
+gasshare_cendiv(cendiv,t)$[sum{cendiv2,repgasquant(cendiv2,t)}] = 
+  repgasquant(cendiv,t) / sum{cendiv2,repgasquant(cendiv2,t)} ;
 
-gascost_cendiv(cendiv,t)$tmodel_new(t) =
+* cost of natural gas - standalone ReEDS
+gascost_cendiv(cendiv,t)$[tmodel_new(t)$tfuel(t)] =
 *cost of natural gas for Sw_GasCurve = 2 (static natural gas prices)
               + sum{(i,v,r,h)$[r_cendiv(r,cendiv)$valgen(i,v,r,t)$gas(i)$heat_rate(i,v,r,t)
                               $[not bio(i)]$[not cofire(i)]$[Sw_GasCurve = 2]],
@@ -595,11 +622,20 @@ gascost_cendiv(cendiv,t)$tmodel_new(t) =
 
               )$[Sw_GasCurve = 1];
 
+* cost of natural gas - linked with FINITO ('not tfuel' indicates years using FINITO supply curves)
+gascost_cendiv(cendiv,t)$[tmodel_new(t)$(not tfuel(t))] =
+* cost = gas price multiplied by gas usage [$ = $/MMBtu * MMBtu/MWh * MW * h]
+    sum{(i,v,r,h)$[valgen(i,v,r,t)$gas(i)$r_cendiv(r, cendiv)], 
+        repgasprice_finito(cendiv,h,t) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) * hours(h) 
+    }
+;
+
+
 *========================================
 * BIOFUEL COSTS
 *========================================
 
-bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)] =
+bioshare_techba(i,r,t)$[(cofire(i) or bio(i))$tmodel_new(t)$tfuel(t)] =
 *  biofuel-based generation of tech i in the BA (biopower + cofire)
                 ((   sum{(v,h)$[valgen(i,v,r,t)$bio(i)], hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) }
                    + sum{(v,h)$[cofire(i)$valgen(i,v,r,t)], bio_cofire_perc * hours(h) * heat_rate(i,v,r,t) * GEN.l(i,v,r,h,t) }
@@ -1740,6 +1776,11 @@ error_check('z') = (
         )
 * account for penalty paid to deploy capacity beyond interconnection queue limits        
         + sum{(tg,r), cap_penalty(tg) * CAP_ABOVE_LIM.l(tg,r,t) }  
+* account for costs from FINITO: deflate from $2018 to $2004,
+* remove any FINITO scaling, and then apply ReEDS scaling
+$ifthene.linked_objective Sw_FINITO_Link==1
+        + cost_scale * ( Z_finito.l(t)$t_finito(t) * deflator('2018') / obj_scale )
+$endif.linked_objective 
     }
 ) / z.l ;
 
@@ -2033,6 +2074,10 @@ load_cat("h2_network",r,t)$tmodel_new(t) =
 ;
 load_cat("dac",r,t)$tmodel_new(t) = sum{i$dac(i), prod_load_ann(i,r,t) } ;
 
+$ifthene.finitoloadcat Sw_FINITO_Link == 1
+load_cat("finito",r,t)$tmodel_new(t) = sum{h, hours(h) * USE_ELE_FINITO.l(r,h,t) } ;
+$endif.finitoloadcat
+
 *========================================
 * H2 NETWORK
 *========================================
diff --git a/reeds/core/terminus/report_dump.py b/reeds/core/terminus/report_dump.py
index 7e3ef747..617433c9 100644
--- a/reeds/core/terminus/report_dump.py
+++ b/reeds/core/terminus/report_dump.py
@@ -66,6 +66,8 @@ def dfdict_to_h5(
     """
     Write dictionary of dataframes to one .h5 file
     """
+    print(f"Saving results to {filepath}")
+
     ### unless a subset is specified, iterate over all keys in dict
     _symbol_list = dfdict.keys() if symbol_list is None else symbol_list
 
@@ -107,6 +109,8 @@ def dfdict_to_excel(
     """
     Write dictionary of dataframes to one .xlsx file
     """
+    print(f"Saving results to {filepath}")
+
     ### unless a subset is specified, iterate over all keys in dict
     _symbol_list = dfdict.keys() if symbol_list is None else symbol_list
 
@@ -275,6 +279,19 @@ def postprocess_outputs(case, outputs_path=None, verbose=0):
     )
     print("Finished loading outputs gdx")
 
+    ### FINITO outputs when running linked model
+    if int(sw.GSw_FINITO_Link):
+        print("Loading FINITO outputs gdx")
+        finito_gdx = os.path.join(outputs_path, f"finito_reeds_outputs_{os.path.basename(case)}.gdx")
+        try:
+            finito_outputs = gdxpds.to_dataframes(finito_gdx)
+            dict_out.update(finito_outputs)
+            print("Finished loading FINITO outputs gdx")
+        except Exception as err:
+            print(err)
+            print(f"Error loading FINITO outputs in {finito_gdx}, skipping.")
+
+
     write_dfdict(
         dfdict=dict_out,
         outputs_path=outputs_path,
diff --git a/reeds/core/terminus/report_params.csv b/reeds/core/terminus/report_params.csv
index 92d20d84..33809d46 100644
--- a/reeds/core/terminus/report_params.csv
+++ b/reeds/core/terminus/report_params.csv
@@ -142,6 +142,7 @@ raw_op_cost(t),2004$,sum of operational costs from systemcost,,,
 RE_gen_price_nat(t),2004$/MWh,marginal cost of the national RE generation constraint,,,
 "repbioprice(r,t)",2004$/MMBtu,highest marginal bioprice of utilized bins for each region,,,
 "repgasprice(cendiv,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each census division,,,
+"repgasprice_finito(cendiv,allh,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins by hour for each census division when using linked model,,,
 "repgasprice_r(r,t)",2004$/MMBtu,highest marginal gas price of utilized gas bins for each region,1,,
 repgasprice_nat(t),2004$/MMBtu,weighted-average national natural gas price assuming that plants pay the marginal price,,,
 "repgasquant(cendiv,t)",Quads,quantity of gas consumed in each census division,,,
diff --git a/reeds/input_processing/copy_files.py b/reeds/input_processing/copy_files.py
index 5a74bdfa..bf573e79 100644
--- a/reeds/input_processing/copy_files.py
+++ b/reeds/input_processing/copy_files.py
@@ -1528,6 +1528,11 @@ def main(reeds_path, inputs_case):
     
     source_deflator_map = get_source_deflator_map(reeds_path)
 
+    ### Rewrite the  finito scalar as GAMS-readable definitions
+    if int(sw['GSw_FINITO_Link']) == 1:
+        pd.read_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv')).iloc[:, :-1].to_csv(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv'),index=False)
+        scalar_csv_to_txt(os.path.join(os.path.dirname(inputs_case),'finito','inputs','scalars.csv'))
+    
     # Copy non-region files
     write_non_region_files(non_region_files, sw, inputs_case, regions_and_agglevel, source_deflator_map)
     
diff --git a/reeds/input_processing/hourly_load.py b/reeds/input_processing/hourly_load.py
index ce4b8e0a..8d44e8ca 100644
--- a/reeds/input_processing/hourly_load.py
+++ b/reeds/input_processing/hourly_load.py
@@ -596,6 +596,44 @@ def reaggregate_to_model_regions(
 
     return regional_load_hourly
 
+#TODO: add docstrings
+def get_hourly_finito_load(
+    inputs_case: str, 
+) -> pd.DataFrame:
+    # load reference FINITO load
+    inputs_case_finito = Path(inputs_case).parent / 'finito' / Path(inputs_case).name
+    load_finito = pd.read_csv(inputs_case_finito / "load_finito.csv")
+    
+    # reshape to match load data. with year in index and region in columns
+    load_finito = load_finito.melt(id_vars='r', var_name='year', value_name='load_MWh')
+    load_finito.year = load_finito.year.astype(int)    
+    load_finito = load_finito.pivot(index='year', columns='r', values='load_MWh')
+    load_finito.columns.name = None
+
+    # allocate annual load to hours, assuming flat demand
+    # TODO: should this use h_weight_finito?
+    hours_per_year = 8760
+    load_hourly_finito = load_finito / hours_per_year
+    load_hourly_finito = load_hourly_finito.astype(np.float32)
+
+    return load_hourly_finito
+
+def remove_finito_load(
+    load_hourly: pd.DataFrame,
+    inputs_case: str, 
+) -> pd.DataFrame:
+
+    # get FINITO reference load
+    load_hourly_finito = get_hourly_finito_load(inputs_case)
+
+    # subtract FINITO reference load from ReEDS load data, 
+    # aligning by model year (index) and region (columns)
+    result = load_hourly - load_hourly_finito
+
+    #TODO: add this
+    # Validation check: confirm that total load - FINITO load > 0
+    
+    return load_hourly
 
 #%% ===========================================================================
 ### --- MAIN FUNCTION ---
@@ -697,6 +735,15 @@ def main(reeds_path, inputs_case):
     peakload = calculate_peak_load(regional_load_hourly, hierarchy)
 
     #%%%#########################################
+    #    -- FINITO Load Adjustment --           #
+    #############################################
+
+    # note that this step occurs after peakload calculation so that the latter
+    # includes a baseline estimate of industrial load captured by FINITO
+    if int(sw.GSw_FINITO_Link):
+        regional_load_hourly = remove_finito_load(regional_load_hourly, inputs_case)
+
+    #############################################
     #    -- DR Shed Load Modifications --    #
     #############################################
 
diff --git a/reeds/input_processing/hourly_repperiods.py b/reeds/input_processing/hourly_repperiods.py
index d45142b8..276f21e4 100644
--- a/reeds/input_processing/hourly_repperiods.py
+++ b/reeds/input_processing/hourly_repperiods.py
@@ -36,6 +36,8 @@
 import reeds
 from reeds.input_processing import hourly_writetimeseries
 from reeds.input_processing import hourly_plots
+from reeds.input_processing import hourly_load
+
 ## Time the operation of this script
 tic = datetime.datetime.now()
 
@@ -50,11 +52,18 @@
 techs_min_vre = ['upv', 'wind-ons']
 
 #%%### Functions
-def get_load(inputs_case, keep_modelyear=None, keep_weatheryears=[2012]):
+def get_load(inputs_case, sw, keep_modelyear=None, keep_weatheryears=[2012]):
     """
     """
     ### Subset to modeled regions
     load = reeds.io.read_file(os.path.join(inputs_case,'load.h5'), parse_timestamps=True)
+
+    ### When running the linked model (GSw_FINITO_Link=1) we add reference load estimates 
+    ### for FINITO load back in to use when selecting representative periods
+    if int(sw.GSw_FINITO_Link):
+        load_hourly_finito = hourly_load.get_hourly_finito_load(inputs_case)
+        load = load + load_hourly_finito
+
     ### Subset to keep_modelyear if provided
     if keep_modelyear:
         load = load.loc[keep_modelyear].copy()
@@ -339,6 +348,7 @@ def main(
     print("Collecting 8760 load data")
     load = get_load(
         inputs_case=inputs_case,
+        sw=sw,
         keep_modelyear=(int(sw['GSw_HourlyClusterYear'])
                   if int(sw['GSw_HourlyClusterYear']) in modelyears
                   else max(modelyears)),
@@ -497,7 +507,7 @@ def main(
         and (sw['GSw_PRM_StressSeedLoadLevel'].lower() not in ['false','none'])
     ):
         ## Get load for all model and weather years
-        load_allyears = get_load(inputs_case, keep_weatheryears='all').loc[modelyears]
+        load_allyears = get_load(inputs_case, sw, keep_weatheryears='all').loc[modelyears]
         ## Add descriptive index
         load_allyears = load_allyears.merge(
             timestamps[['year', 'yperiod', 'h_of_period']], left_on='datetime', right_index=True)
@@ -692,10 +702,13 @@ def main(
         description='Create the necessary 8760 and capacity factor data for hourly resolution')
     parser.add_argument('reeds_path', help='ReEDS directory')
     parser.add_argument('inputs_case', help='ReEDS/runs/{case}/inputs_case directory')
+    parser.add_argument('--nolog', '-n', default=False, action='store_true', help='turn off logging for debugging')
 
     args = parser.parse_args()
     reeds_path = args.reeds_path
     inputs_case = args.inputs_case
+    logging = not args.nolog
+
 
     # #%% Settings for testing
     # reeds_path = reeds.io.reeds_path
@@ -705,10 +718,11 @@ def main(
     # interactive = True
 
     #%% Set up logger
-    log = reeds.log.makelog(
-        scriptname=__file__,
-        logpath=os.path.join(inputs_case,'..','gamslog.txt'),
-    )
+    if logging:
+        log = reeds.log.makelog(
+            scriptname=__file__,
+            logpath=os.path.join(inputs_case,'..','gamslog.txt'),
+        )
     print('Starting hourly_repperiods.py')
     #%% Inputs from switches
     sw = reeds.io.get_switches(inputs_case)
@@ -723,6 +737,7 @@ def main(
         sw=sw, reeds_path=reeds_path, inputs_case=inputs_case,
         periodtype='rep',
         make_plots=1,
+        logging=logging
     )
 
     ############################################
diff --git a/reeds/input_processing/hourly_writetimeseries.py b/reeds/input_processing/hourly_writetimeseries.py
index 8ac35efb..ea50306f 100644
--- a/reeds/input_processing/hourly_writetimeseries.py
+++ b/reeds/input_processing/hourly_writetimeseries.py
@@ -200,6 +200,9 @@ def get_yearly_demand(sw, hmap_myr, hmap_allyrs, inputs_case, periodtype='rep'):
     reload the raw demand and extract the demand on the modeled days for each year.
     """
     ### Get original demand data, subset to cluster year
+    ### Note that this does not include FINITO demand when running with the linked model
+    ### (GSw_FINITO_Link=1) even though the reference estimates for that load is used when
+    ### identifying the representative days
     load_in = reeds.io.read_file(
         os.path.join(inputs_case,'load.h5'), parse_timestamps=True).unstack(level=0)
     load_in.columns = load_in.columns.rename(['r','t'])
diff --git a/reeds/inputs.py b/reeds/inputs.py
index e598dc26..f7cee045 100644
--- a/reeds/inputs.py
+++ b/reeds/inputs.py
@@ -5,6 +5,7 @@
 import yaml
 import hashlib
 import shapely
+import shutil
 import numpy as np
 import pandas as pd
 import sklearn.cluster
@@ -371,6 +372,7 @@ def solvestring_sequential(
             'GSw_StateCO2ImportLevel',
             'GSw_StartMarkets',
             'GSw_ValStr',
+            'GSw_FINITO_Link',
             'solver',
             'debug',
             'startyear',
@@ -912,3 +914,57 @@ def validate_zoneset(GSw_ZoneSet):
                 "to ensure each aggreg is only assigned to a single hierarchy level."
             )
             raise ValueError(err)
+
+def setup_finito(casedir, caseSwitches, BatchName):
+    #%% Copy FINITO code folders and inputs into [casedir]/finito
+    # ... finito directory within the case directory
+    casedir_finito = os.path.join(casedir,'finito')
+    # ... define the inputs case directory for FINITO
+    inputs_case_finito = os.path.join(casedir,'finito','inputs_case')
+
+    # copy directories
+    os.makedirs(casedir_finito, exist_ok=True)
+    shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'inputs'),os.path.join(casedir,'finito', 'inputs'))
+    shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'input_processing'),os.path.join(casedir,'finito', 'input_processing'))
+    shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'model'),os.path.join(casedir,'finito', 'model'))
+    shutil.copytree(os.path.join(caseSwitches['finito_dir'], 'visualization'),os.path.join(casedir, 'finito', 'visualization'))
+    
+    # copy over the FINITO cases files
+    shutil.copy2(os.path.join(caseSwitches['finito_dir'], 'cases.csv'), os.path.join(casedir, 'finito'))
+    shutil.copy2(os.path.join(caseSwitches['finito_dir'], f"cases_{caseSwitches['finito_cases_file']}.csv"), os.path.join(casedir, 'finito'))
+
+    #%% (GSw_Trade_PriceResponse > 0) If doing a price-responsive trade run, retrieve the reference exports/imports prices
+    if int(caseSwitches['GSw_Trade_PriceResponse']) > 0:
+        initialize_price_response_path = os.path.join(casedir_finito, 'input_processing', 'processing', 'initialize_price_response.py')
+        # Collect all arguments for initialize_price_response.py
+        initialize_price_response_args = f" -c {casedir} -b {BatchName} -cr {caseSwitches['GSw_Trade_PriceResponse_RefScen']} -l {caseSwitches['GSw_FINITO_Link']}"
+        # Call initialize_price_response.py file before starting the runs
+        os.system('python ' + initialize_price_response_path + initialize_price_response_args)
+
+    #%% Filter and copy all input files for each scenario
+    # Call FINITO copy_files.py file before starting the runs
+    copy_files_run = os.system(
+        'python ' + os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'copy_files.py') + 
+        f" -c {casedir_finito} -d {inputs_case_finito} --link"
+        )
+
+    # Print an error message if copy_files.py encounters any issue 
+    if copy_files_run != 0:
+        print("\nERROR: FINITO copy_files.py encountered an issue and did not complete successfully.")
+        print("Please check the console output above for details.")
+        print("The issue could be due to regionality, focus sector filtering, or file reading errors.") 
+        os._exit(1)   
+    
+    ## Populate sets for each linked run using autopop_set.py
+    #autopop_args = f" -c {casedir_finito} -d {inputs_case_finito} -a {aeo_year} -s {focus_sectors} -f {caseSwitches['GSw_FixedCostSupply']} -rwf {caseSwitches['GSw_Trade_Partners']} -e {caseSwitches['GSw_ROE_EndUses']}"
+    os.system(
+        'python ' + os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'autopop_set.py') + 
+        f" -c {casedir_finito} -d {inputs_case_finito} --link"
+    )
+
+    ## Call read_mecs_heat.py to generate heat/nonheat/feedstock ratios for FINITO Rest of Industry (ROI)
+    mecs_sectors = caseSwitches['focus_sectors'].replace('.', ' ') 
+    read_mecs_path = os.path.join(caseSwitches['finito_dir'], 'input_processing', 'processing', 'mecs', 'read_mecs_heat.py')
+    # Collect all arguments for read_mecs_heat.py
+    read_mecs_args = f' -s {mecs_sectors} -d {inputs_case_finito}'
+    os.system('python ' + read_mecs_path + read_mecs_args)
\ No newline at end of file
diff --git a/reeds/reedsplots.py b/reeds/reedsplots.py
index cfb3db9c..956c22ca 100644
--- a/reeds/reedsplots.py
+++ b/reeds/reedsplots.py
@@ -5023,6 +5023,7 @@ def plot_seed_stressperiods(
     ### Recalculate peak load days since we dropped duplicates above
     load_allyears = hourly_repperiods.get_load(
         os.path.join(case, 'inputs_case'),
+        sw,
         keep_weatheryears='all').loc[years]
     timestamps = pd.read_csv(os.path.join(case, 'inputs_case', 'rep', 'timestamps.csv'))
     resource_adequacy_years = [int(y) for y in sw.resource_adequacy_years.split('_')]
diff --git a/reeds/resource_adequacy/prep_data.py b/reeds/resource_adequacy/prep_data.py
index eebb2f0a..f351c117 100644
--- a/reeds/resource_adequacy/prep_data.py
+++ b/reeds/resource_adequacy/prep_data.py
@@ -147,8 +147,38 @@ def main(t, casedir, iteration=0):
         h_dt_szn.index.map(hmap_allyrs.set_index(['year', 'hour'])['*timestamp']))
     h_dt_szn = h_dt_szn.reset_index().set_index('timestamp')
 
+    # load exogenous demand seen by ReEDS
     load = reeds.io.read_file(os.path.join(inputs_case, 'load.h5'), parse_timestamps=True)
 
+    load_year = load.loc[t]
+
+    # when running linked model with FINITO (GSw_FINITO_Link=1), 
+    # we also need to add in the load from FINITO
+    if int(sw.GSw_FINITO_Link):
+        load_finito_rt = gdxreeds['load_finito_rt'].rename(columns={'allh':'h', 'Value':'load_MW'})
+
+        # down select to relevant model year
+        load_finito = load_finito_rt.loc[load_finito_rt.t==t].drop('t', axis=1).copy()
+        #TODO: why so much shifting? validate the dynamics here with FINITO team
+        #TODO: also check against reference quantity
+
+        # map from rep day to actual hour
+        # since we don't have multi-year profiles for FINITO load 
+        # we assume they repeat across all weather years
+        # TODO: test that this works with multiple regions
+        load_finito = pd.merge(load_finito, h_dt_szn.reset_index(), on='h', how='outer')[['timestamp', 'r', 'load_MW']]
+        load_finito = load_finito.rename(columns={'timestamp':'datetime'})
+        # add model year back in 
+        #load_finito = load_finito.assign(year=t)
+        load_finito = load_finito.pivot(index=['datetime'], columns='r', values='load_MW')
+
+        # convert timezone and fill any missing columns
+        load_finito.index = load_finito.index.tz_convert(load_year.index.tz)
+        load_finito = load_finito.reindex(columns=load_year.columns, fill_value=0)
+
+        # add to load
+        load_new = load_year + load_finito
+
     resources = pd.read_csv(os.path.join(inputs_case, 'resources.csv'))
     recf = reeds.io.read_file(os.path.join(inputs_case, 'recf.h5'), parse_timestamps=True)
     recf.columns = pd.MultiIndex.from_tuples([tuple(x.split('|')) for x in recf.columns],
@@ -364,7 +394,7 @@ def intify(v):
         .merge(h_dt_szn_load_years[['h']].reset_index(), left_on='allh', right_on='h')
         .pivot(index='timestamp', columns='r', values='Value')
     )
-    load_year = load.loc[t].add(can_exports, fill_value=0)
+    load_year = load_year.add(can_exports, fill_value=0)
 
     ### PRAS doesn't yet handle flexible load, so include all H2/DAC load in the
     ### version we write for PRAS
diff --git a/runreeds.py b/runreeds.py
index 1cbfc517..0b07da7f 100644
--- a/runreeds.py
+++ b/runreeds.py
@@ -4,6 +4,7 @@
 
 import reeds
 import os
+import sys
 import git
 import queue
 import threading
@@ -99,12 +100,23 @@ def create_case_lists(df_cases:pd.DataFrame, BatchName:str, single:str=''):
                 continue
         # Add switch settings to list of options passed to GAMS
         shcom = f' --case={BatchName}_{case}'
-        for i,v in df_cases[case].items():
+        # (ReEDS-FINITO) Combine the cases files for the linked model
+        if int(df_cases[case].loc['GSw_FINITO_Link']) == 1 :
+            # add the FINITO switches to the caseSwitches
+            df_cases_combine=linked_cases(df_cases,case)
+            for i,v in df_cases_combine.items():
+                #exclude certain switches that don't need to be passed to GAMS
+                if i not in ['file_replacements','keep_run_terminal']:
+                    shcom = shcom + ' --' + i + '=' + v
+            caseList.append(shcom)
+            caseSwitches.append(df_cases_combine.to_dict())
+        else:
+            for i,v in df_cases[case].items():
             #exclude certain switches that don't need to be passed to GAMS
-            if i not in ['file_replacements','keep_run_terminal']:
-                shcom += f' --{i}={v}'
-        caseList.append(shcom)
-        caseSwitches.append(df_cases[case].to_dict())
+                if i not in ['file_replacements','keep_run_terminal']:
+                    shcom += f' --{i}={v}'
+            caseList.append(shcom)
+            caseSwitches.append(df_cases[case].to_dict())
 
     return caseSwitches, casenames, caseList
 
@@ -755,6 +767,71 @@ def setup_window(
             OPATH.writelines( "python valuestreams.py" + '\n')
 
 
+# (ReEDS-FINITO) update 'df_cases' to include FINITO switches
+def linked_cases(df_cases,case):
+    """
+    Updates the cases dataframe to include FINITO-specific switches.
+    When a switch is duplicated in FINITO and ReEDS, then we default to 
+    the ReEDS value.
+    
+    For the FINITO switches, the combined cases file defaults to the case-specific 
+    'Default Value' in cases_linked.csv first, before using the universal FINITO 
+    'Default Value' in cases.csv for any un-assigned switches.    
+    """
+    # check for valid finito_dir
+    if not os.path.isdir(df_cases[case]['finito_dir']):
+        raise ValueError(f"finito_dir = {df_cases[case]['finito_dir']} is not a valid path.")
+
+    # define path to and read the FINITO check_inputs function
+    finito_check_inputs_path = os.path.join(df_cases[case]['finito_dir'], 'input_processing', 'processing')
+    sys.path.append(finito_check_inputs_path)
+    from check_inputs import check_inputs
+        
+    ### load the default values for all FINITO switches from ~\FINITO\cases.csv
+    df_cases_finito = pd.read_csv(os.path.join(df_cases[case]['finito_dir'],'cases.csv'), dtype=object, index_col=0)
+    df_cases_finito = df_cases_finito[['Choices', 'Default Value']]
+
+    ### load the scenario-specific switches from ~\FINITO\cases_linked.csv
+    cases_linked_path = os.path.join(df_cases[case]['finito_dir'],f"cases_{df_cases[case]['finito_cases_file']}.csv")
+    df_cases_suf_finito = pd.read_csv(cases_linked_path, dtype=object, index_col=0)
+    ## check that case names are unique in cases_linked.csv
+    # grab the scenario names **exactly** as they appear in the csv file
+    with open(cases_linked_path, 'r', newline='') as csvfile:
+        reader = csv.reader(csvfile)
+        header = next(reader)
+    # find the duplicate column names and raise an error if any are found
+    duplicate_columns = set([x for x in header if header.count(x) > 1])
+    if duplicate_columns:
+        raise ValueError(f"The FINITO cases_{df_cases[case]['finito_cases_file']}.csv has the following duplicate column names: {duplicate_columns}")
+    ### identify the FINITO case 
+    if df_cases[case]['finito_case'] == 'same':
+        finito_case=case
+    else:
+        finito_case=df_cases[case]['finito_case']
+    # ensures **exact** match of names between the ReEDS cases_{}.csv and the FINITO cases_linked.csv
+    if finito_case not in (df_cases_suf_finito.columns):
+        raise ValueError(f"The 'finito_case' input '{finito_case}' in the ReEDS cases file does not exist in FINITO's cases_{df_cases[case]['finito_cases_file']}.csv.")
+
+    ### first use 'Default Value' from the FINITO cases_linked.csv to fill missing switches
+    if 'Default Value' in df_cases_suf_finito.columns:
+        df_cases_suf_finito[finito_case] = df_cases_suf_finito[finito_case].fillna(df_cases_suf_finito['Default Value'])
+    ### then, use 'Default Value' from the FINITO cases.csv to fill un-assigned switches
+    df_cases_suf_finito.drop(['Choices','Default Value'], axis='columns',inplace=True, errors='ignore')
+    df_cases_finito = df_cases_finito.join(df_cases_suf_finito, how='outer')
+    df_cases_finito[finito_case] = df_cases_finito[finito_case].fillna(df_cases_finito['Default Value'])
+  
+    #### create new dataframe for the combined ReEDS and FINITO switches
+    df_cases_combine = pd.concat([df_cases[case],df_cases_finito[finito_case]])
+    ### drop duplicated switches, defaulting to reeds
+    df_cases_combine = df_cases_combine[~df_cases_combine.index.duplicated(keep='first')]
+
+    #%% Check for incompatibility of FINITO switches
+    model_sectors = df_cases_finito['Default Value']['focus_sectors'].split('.')
+    check_inputs(case = case, df_case = df_cases_combine, model_sectors=model_sectors)
+
+    return df_cases_combine
+
+
 #%% ===========================================================================
 ### --- PROCEDURE ---
 ### ===========================================================================
@@ -1237,6 +1314,10 @@ def write_batch_script(
     os.makedirs(os.path.join(casedir,'handoff','reeds_data'), exist_ok=True)
     os.makedirs(os.path.join(casedir,'handoff','PRAS'), exist_ok=True)
 
+    ## Set up FINITO if running linked model
+    if int(caseSwitches['GSw_FINITO_Link'])==1:
+        reeds.inputs.setup_finito(casedir, caseSwitches, BatchName)
+
     ###### Replace files according to 'file_replacements' in cases. Ignore quotes in input text.
     # << is used to separate the file that is to be replaced from the file that is used
     # || is used to separate multiple replacements.
@@ -1383,12 +1464,43 @@ def write_batch_script(
             + ' gdxcompress=1'
             + toLogGamsString
             + f"--fname={batch_case}"
-            + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']} \n"
+            + f" --GSw_calc_powfrac={caseSwitches['GSw_calc_powfrac']}"
+            + f" --first_year_finito={caseSwitches['first_year_finito']} \n"
         )
         OPATH.writelines(writescripterrorcheck("report.gms"))
-        if not LINUXORMAC:
+        if not LINUXORMAC and int(caseSwitches['GSw_FINITO_Link']) != 1:
             OPATH.writelines("endlocal\n")
         OPATH.writelines(f'python {logger}\n')
+
+        ### (ReEDS-FINITO) calls FINITO reporting
+        if int(caseSwitches['GSw_FINITO_Link'])==1:
+            OPATH.writelines(
+            'gams '
+            + f"{os.path.join(casedir,'finito', 'model', 'finito_report.gms')}"
+            + f" o={os.path.join('lstfiles',f'finito_report_{batch_case}.lst')}"
+            + (' license=gamslice.txt' if hpc else '')
+            + (' r=$r' if LINUXORMAC else ' r=!r!')
+            + ' gdxcompress=1'
+            + toLogGamsString
+            + f"--case={batch_case}"
+            + f" --casedir={casedir}"
+            + f" --GSw_FINITO_Link={caseSwitches['GSw_FINITO_Link']}"
+            + f" --GSw_RetailAdder={caseSwitches['GSw_RetailAdder']}"
+            + f" --finito_inputs_dir={os.path.join('finito','inputs')}"
+            + f" --inputs_case_finito_dir={os.path.join('finito','inputs_case')}"
+            + f" --linked_report_dir={caseSwitches['linked_report_dir']} \n"
+            )
+            # (ReEDS-FINITO) calls FINITO postprocessing
+            OPATH.writelines(
+            f"python {os.path.join(casedir, 'finito', 'visualization', 'postprocessing.py')}"
+            + " -b 0"            
+            + f" -l {caseSwitches['GSw_FINITO_Link']}"
+            + f' -c {batch_case} \n\n'
+            )  
+            if not LINUXORMAC:
+                OPATH.writelines("endlocal\n")
+
+        ### save reporting outputs to h5 and/or csv files
         OPATH.writelines(f"python {Path('reeds','core','terminus','report_dump.py')} {casedir} -c\n\n")
         if int(caseSwitches['diagnose']):
              OPATH.writelines(