diff --git a/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_double_loop_usc.py b/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_double_loop_usc.py
index 2846b94a..d4070aab 100644
--- a/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_double_loop_usc.py
+++ b/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_double_loop_usc.py
@@ -1,7 +1,7 @@
 #################################################################################
-# DISPATCHES was produced under the DOE Design Integration and Synthesis
-# Platform to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES),
-# and is copyright (c) 2022 by the software owners: The Regents of the University
+# DISPATCHES was produced under the DOE Design Integration and Synthesis Platform
+# to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES), and is
+# copyright (c) 2020-2023 by the software owners: The Regents of the University
 # of California, through Lawrence Berkeley National Laboratory, National
 # Technology & Engineering Solutions of Sandia, LLC, Alliance for Sustainable
 # Energy, LLC, Battelle Energy Alliance, LLC, University of Notre Dame du Lac, et
@@ -10,12 +10,11 @@
 # Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
 # information, respectively. Both files are also available online at the URL:
 # "https://github.com/gmlc-dispatches/dispatches".
-#
 #################################################################################
 
 """
 This script describes a multiperiod class to build an object for the
-integrated ultra-supercritical power plant and molten-salt based
+integrated ultra-supercritical power plant and motlen-salt based
 thermal energy storage model.
 """
 
@@ -24,8 +23,46 @@
 import pyomo.environ as pyo
 import pandas as pd
 from collections import deque
+import idaes.logger as idaeslog
+
+# IDAES imports
+from idaes.apps.grid_integration.multiperiod.multiperiod import MultiPeriodModel
+
+# DISPATCHES imports
 from dispatches.case_studies.fossil_case.ultra_supercritical_plant \
-    .storage.multiperiod_integrated_storage_usc import create_multiperiod_usc_model
+    .storage.multiperiod_integrated_storage_usc import (create_usc_model,
+                                                        usc_unfix_dof,
+                                                        usc_custom_init,
+                                                        get_usc_link_variable_pairs)
+
+
+def create_multiperiod_usc_model(n_time_points=4, pmin=None, pmax=None):
+    """
+    Create a multi-period usc_mp cycle object. This object contains a pyomo
+    model with a block for each time instance.
+
+    n_time_points: Number of time blocks to create
+    """
+    multiperiod_usc = MultiPeriodModel(
+        n_time_points=n_time_points,
+        process_model_func=create_usc_model,
+        initialization_func=usc_custom_init,
+        unfix_dof_func=usc_unfix_dof,
+        linking_variable_func=get_usc_link_variable_pairs,
+        use_stochastic_build=False,
+        outlvl=idaeslog.INFO,
+        )
+
+    flowsheet_options={"pmin": pmin,
+                        "pmax": pmax}
+
+    # create the multiperiod object
+    multiperiod_usc.build_multi_period_model(
+        model_data_kwargs={t: flowsheet_options for t in range(n_time_points)},
+        flowsheet_options=flowsheet_options,
+    )
+
+    return multiperiod_usc
 
 
 class MultiPeriodUsc:
@@ -34,7 +71,7 @@ def __init__(
     ):
         """
         Arguments:
-            horizon: Int64 - number of time points to use for associated multi-period model
+            horizon::Int64 - number of time points to use for associated multi-period model
 
         Returns:
             Float64: Value of power output in last time step
@@ -44,13 +81,13 @@ def __init__(
         self.result_listimp = []
         self.model_data = model_data
 
-    def populate_model(self, b, horizon):
+    def populate_model(self, blk, horizon):
         """
-        Create an integrated ultra-supercritical power plant and molten salt
+        Create a integrated ultra-supercritical power plant and molten salt
         thermal energy storage model using the `MultiPeriod` package.
 
         Arguments:
-            blk: this is an empty block passed in from either a bidder or tracker
+            blk: this is an empty block passed in from eithe a bidder or tracker
 
         Returns:
              None
@@ -58,7 +95,6 @@ def populate_model(self, b, horizon):
         tank_min = 76000  # in kg
         tank_max = 6739292  # in kg
 
-        blk = b
         if not blk.is_constructed():
             blk.construct()
 
@@ -69,9 +105,9 @@ def populate_model(self, b, horizon):
         blk.usc_mp = multiperiod_usc
 
         active_blks = multiperiod_usc.get_active_process_blocks()
-        active_blks[0].usc_mp.previous_salt_inventory_hot.fix(tank_min)
-        active_blks[0].usc_mp.previous_salt_inventory_cold.fix(tank_max-tank_min)
-        active_blks[0].usc_mp.previous_power.fix(380)
+        active_blks[0].fs.previous_salt_inventory_hot.fix(tank_min)
+        active_blks[0].fs.previous_salt_inventory_cold.fix(tank_max-tank_min)
+        active_blks[0].fs.previous_power.fix(380)
 
         # create expression that references underlying power variables
         blk.HOUR = pyo.Set(initialize=range(horizon))
@@ -92,33 +128,34 @@ def populate_model(self, b, horizon):
         blk.hxd_steam_Tout = pyo.Expression(blk.HOUR)
         blk.hxd_steam_vfrac = pyo.Expression(blk.HOUR)
         for (t, b) in enumerate(active_blks):
-            blk.P_T[t] = b.usc_mp.fs.net_power
-            blk.hot_level[t] = b.usc_mp.salt_inventory_hot
+            blk.P_T[t] = b.fs.net_power
+            blk.hot_level[t] = b.fs.salt_inventory_hot
             blk.storage_power[t] = ((-1e-6)
-                                 * b.usc_mp.fs.es_turbine.work_mechanical[0])
-            blk.plant_duty[t] = b.usc_mp.fs.plant_heat_duty[0]
+                                 * b.fs.es_turbine.work_mechanical[0])
+            blk.plant_duty[t] = b.fs.plant_heat_duty[0]
             blk.tot_cost[t] = (
-                     b.usc_mp.fs.operating_cost
-                     + (b.usc_mp.fs.plant_fixed_operating_cost
-                     + b.usc_mp.fs.plant_variable_operating_cost) / (365 * 24)
+                     b.fs.operating_cost
+                     + (b.fs.plant_fixed_operating_cost
+                     + b.fs.plant_variable_operating_cost) / (365 * 24)
                 )
-            blk.plant_power[t] = b.usc_mp.fs.plant_power_out[0]
-            blk.hxc_salt[t] = b.usc_mp.fs.hxc.tube_inlet.flow_mass[0]
-            blk.hxc_duty[t] = b.usc_mp.fs.hxc.heat_duty[0]
-            blk.hxc_salt_Tin[t] = b.usc_mp.fs.hxc.tube_inlet.temperature[0]
-            blk.hxc_salt_Tout[t] = b.usc_mp.fs.hxc.tube_outlet.temperature[0]
-            blk.hxd_salt[t] = b.usc_mp.fs.hxd.shell_inlet.flow_mass[0]
-            blk.hxd_duty[t] = b.usc_mp.fs.hxd.heat_duty[0]
-            blk.hxd_salt_Tin[t] = b.usc_mp.fs.hxd.shell_inlet.temperature[0]
-            blk.hxd_salt_Tout[t] = b.usc_mp.fs.hxd.shell_outlet.temperature[0]
-            blk.hxd_steam_Tout[t] = b.usc_mp.fs.hxd.cold_side.properties_out[0].temperature
-            blk.hxd_steam_vfrac[t] = b.usc_mp.fs.hxd.cold_side.properties_out[0].vapor_frac
+            blk.plant_power[t] = b.fs.plant_power_out[0]
+            blk.hxc_salt[t] = b.fs.hxc.tube_inlet.flow_mass[0]
+            blk.hxc_duty[t] = b.fs.hxc.heat_duty[0]
+            blk.hxc_salt_Tin[t] = b.fs.hxc.tube_inlet.temperature[0]
+            blk.hxc_salt_Tout[t] = b.fs.hxc.tube_outlet.temperature[0]
+            blk.hxd_salt[t] = b.fs.hxd.shell_inlet.flow_mass[0]
+            blk.hxd_duty[t] = b.fs.hxd.heat_duty[0]
+            blk.hxd_salt_Tin[t] = b.fs.hxd.shell_inlet.temperature[0]
+            blk.hxd_salt_Tout[t] = b.fs.hxd.shell_outlet.temperature[0]
+            blk.hxd_steam_Tout[t] = b.fs.hxd.cold_side.properties_out[0].temperature
+            blk.hxd_steam_vfrac[t] = b.fs.hxd.cold_side.properties_out[0].vapor_frac
 
 
         self.multiperiod_usc = multiperiod_usc
         return
 
-    def update_model(self, b, implemented_power_output, realized_soc):
+    @staticmethod
+    def update_model(b, implemented_power_output, realized_soc):
 
         """
         Update `blk` variables using the actual implemented power output.
@@ -131,8 +168,7 @@ def update_model(self, b, implemented_power_output, realized_soc):
          Returns:
              None
         """
-        blk = b
-        multiperiod_usc = blk.usc_mp
+        multiperiod_usc = b.usc_mp
         active_blks = multiperiod_usc.get_active_process_blocks()
 
         implemented_power = round(implemented_power_output[-1])
@@ -140,8 +176,8 @@ def update_model(self, b, implemented_power_output, realized_soc):
         print("Implemented Power (MPC)", implemented_power)
         print("Realized SOC (MPC)", realized_soc)
 
-        active_blks[0].usc_mp.previous_power.fix(implemented_power)
-        active_blks[0].usc_mp.previous_salt_inventory_hot.fix(realized_soc)
+        active_blks[0].fs.previous_power.fix(implemented_power)
+        active_blks[0].fs.previous_salt_inventory_hot.fix(realized_soc)
 
         return
 
@@ -157,10 +193,10 @@ def get_last_delivered_power(b, last_implemented_time_step):
                                         step
 
         Returns:
-            Float64: Value of power output in the last time step
+            Float64: Value of power output in last time step
         """
-        blk = b
-        return pyo.value(blk.P_T[last_implemented_time_step])
+        # blk = b
+        return pyo.value(b.P_T[last_implemented_time_step])
 
     @staticmethod
     def get_implemented_profile(b, last_implemented_time_step):
@@ -175,18 +211,18 @@ def get_implemented_profile(b, last_implemented_time_step):
          Returns:
              profile: the intended profile, {unit: [...]}
         """
-        blk = b
-        multiperiod_usc = blk.usc_mp
+        # blk = b
+        multiperiod_usc = b.usc_mp
         active_blks = multiperiod_usc.get_active_process_blocks()
         implemented_power_output = deque(
             [
-                pyo.value(active_blks[t].usc_mp.fs.net_power)
+                pyo.value(active_blks[t].fs.net_power)
                 for t in range(last_implemented_time_step + 1)
             ]
         )
         realized_soc = deque(
             [
-                pyo.value(active_blks[t].usc_mp.salt_inventory_hot)
+                pyo.value(active_blks[t].fs.salt_inventory_hot)
                 for t in range(last_implemented_time_step + 1)
             ]
         )
@@ -196,7 +232,7 @@ def get_implemented_profile(b, last_implemented_time_step):
             "realized_soc": realized_soc,
         }
 
-    def record_results(self, b, date=None, hour=None, **kwargs):
+    def record_results(self, blk, date=None, hour=None, **kwargs):
 
         """
         Record the operations stats for the model.
@@ -210,7 +246,7 @@ def record_results(self, b, date=None, hour=None, **kwargs):
             None
 
         """
-        blk = b
+        # blk = b
         df_list = []
         df_listimp = []
         for t in blk.HOUR:
@@ -220,6 +256,7 @@ def record_results(self, b, date=None, hour=None, **kwargs):
             result_dict["Date"] = date
             result_dict["Hour"] = hour
 
+
             # simulation inputs
             result_dict["Horizon [hr]"] = int(t)
 
diff --git a/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_integrated_storage_usc.py b/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_integrated_storage_usc.py
index 83537e7a..80d4084f 100644
--- a/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_integrated_storage_usc.py
+++ b/dispatches/case_studies/fossil_case/ultra_supercritical_plant/storage/multiperiod_integrated_storage_usc.py
@@ -1,7 +1,7 @@
 #################################################################################
-# DISPATCHES was produced under the DOE Design Integration and Synthesis
-# Platform to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES),
-# and is copyright (c) 2022 by the software owners: The Regents of the University
+# DISPATCHES was produced under the DOE Design Integration and Synthesis Platform
+# to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES), and is
+# copyright (c) 2020-2023 by the software owners: The Regents of the University
 # of California, through Lawrence Berkeley National Laboratory, National
 # Technology & Engineering Solutions of Sandia, LLC, Alliance for Sustainable
 # Energy, LLC, Battelle Energy Alliance, LLC, University of Notre Dame du Lac, et
@@ -10,7 +10,6 @@
 # Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
 # information, respectively. Both files are also available online at the URL:
 # "https://github.com/gmlc-dispatches/dispatches".
-#
 #################################################################################
 
 """
@@ -25,21 +24,23 @@
 
 __author__ = "Naresh Susarla and Soraya Rawlings"
 
-try:
-    from importlib import resources  # Python 3.8+
-except ImportError:
-    import importlib_resources as resources  # Python 3.7
-
+from importlib import resources  # Python 3.8+
+from idaes.core.util import from_json, to_json
+import idaes.logger as idaeslog
 
-from pyomo.environ import (NonNegativeReals, ConcreteModel, Constraint, Var)
+from pyomo.environ import (NonNegativeReals, Constraint, Var, ConcreteModel)
 from idaes.apps.grid_integration.multiperiod.multiperiod import (
     MultiPeriodModel)
+from dispatches.case_studies.fossil_case.ultra_supercritical_plant import (
+    ultra_supercritical_powerplant as usc)
 from dispatches.case_studies.fossil_case.ultra_supercritical_plant.storage import (
-    integrated_storage_with_ultrasupercritical_power_plant as usc)
-from dispatches.case_studies.fossil_case.ultra_supercritical_plant import storage
+    integrated_storage_with_ultrasupercritical_power_plant as usc_w_tes)
+
 
+def create_usc_model(m=None, pmin=None, pmax=None):
 
-def create_usc_model(pmin, pmax):
+    if m is None:
+        m = ConcreteModel(name="Integrated Model")
 
     # Set bounds for plant power
     min_storage_heat_duty = 10  # in MW
@@ -47,70 +48,206 @@ def create_usc_model(pmin, pmax):
 
     max_power = 436  # in MW
     min_power = int(0.65 * max_power)  # 283 in MW
+    if pmin is None:
+        pmin = int(0.65 * 436) + 1
+    if pmax is None:
+        pmax = 436 + 30
+
+    m = usc.build_plant_model(m)
+
+    # Create a flowsheet, add properties, unit models, and arcs
+    m = usc_w_tes.create_integrated_model(m, max_power=max_power)
+
+    # Give all the required inputs to the model
+    usc_w_tes.set_model_input(m)
+
+    # Add scaling factor
+    usc_w_tes.set_scaling_factors(m)
 
-    m = ConcreteModel()
+    # Initialize the model with a sequential initialization and custom
+    # Add cost correlations
+    m = usc_w_tes.build_costing(m)
 
-    with resources.path(storage, "initialized_integrated_storage_usc.json") as data_file_path:
-        assert data_file_path.is_file()
-        m.usc_mp = usc.main(max_power=max_power,
-                            load_from_file=str(data_file_path))
+    # Add bounds
+    usc_w_tes.add_bounds(m)
+    
 
-    m.usc_mp.fs.plant_min_power_eq = Constraint(
-        expr=m.usc_mp.fs.plant_power_out[0] >= min_power
+    m.fs.plant_min_power_eq = Constraint(
+        expr=m.fs.plant_power_out[0] >= min_power
     )
-    m.usc_mp.fs.plant_max_power_eq = Constraint(
-        expr=m.usc_mp.fs.plant_power_out[0] <= max_power
+    m.fs.plant_max_power_eq = Constraint(
+        expr=m.fs.plant_power_out[0] <= max_power
     )
 
-    m.usc_mp.fs.hxc.heat_duty.setlb(min_storage_heat_duty * 1e6)
-    m.usc_mp.fs.hxd.heat_duty.setlb(min_storage_heat_duty * 1e6)
+    m.fs.hxc.heat_duty.setlb(min_storage_heat_duty * 1e6)
+    m.fs.hxd.heat_duty.setlb(min_storage_heat_duty * 1e6)
+
+    m.fs.hxc.heat_duty.setub(max_storage_heat_duty * 1e6)
+    m.fs.hxd.heat_duty.setub(max_storage_heat_duty * 1e6)
+
+    # Add coupling variables
+    m.fs.previous_power = Var(
+        domain=NonNegativeReals,
+        initialize=300,
+        bounds=(pmin, pmax),
+        doc="Previous period power (MW)"
+        )
+
+    inventory_max = 1e7
+    inventory_min = 75000
+    tank_max = 6739292           # Units in kg
 
-    m.usc_mp.fs.hxc.heat_duty.setub(max_storage_heat_duty * 1e6)
-    m.usc_mp.fs.hxd.heat_duty.setub(max_storage_heat_duty * 1e6)
+    m.fs.previous_salt_inventory_hot = Var(
+        domain=NonNegativeReals,
+        initialize=inventory_min,
+        bounds=(0, inventory_max),
+        doc="Hot salt at the beginning of the hour (or time period), kg"
+        )
+    m.fs.salt_inventory_hot = Var(
+        domain=NonNegativeReals,
+        initialize=inventory_min,
+        bounds=(0, inventory_max),
+        doc="Hot salt inventory at the end of the hour (or time period), kg"
+        )
+    m.fs.previous_salt_inventory_cold = Var(
+        domain=NonNegativeReals,
+        initialize=tank_max-inventory_min,
+        bounds=(0, inventory_max),
+        doc="Cold salt at the beginning of the hour (or time period), kg"
+        )
+    m.fs.salt_inventory_cold = Var(
+        domain=NonNegativeReals,
+        initialize=tank_max-inventory_min,
+        bounds=(0, inventory_max),
+        doc="Cold salt inventory at the end of the hour (or time period), kg"
+        )
 
+    @m.fs.Constraint(doc="Plant ramping down constraint")
+    def constraint_ramp_down(b):
+        return (
+            m.fs.previous_power - 60 <=
+            m.fs.plant_power_out[0])
+
+    @m.fs.Constraint(doc="Plant ramping up constraint")
+    def constraint_ramp_up(b):
+        return (
+            m.fs.previous_power + 60 >=
+            m.fs.plant_power_out[0])
+
+    @m.fs.Constraint(doc="Inventory balance at the end of the time period")
+    def constraint_salt_inventory_hot(b):
+        return (
+            m.fs.salt_inventory_hot ==
+            m.fs.previous_salt_inventory_hot
+            + (3600*m.fs.hxc.tube_inlet.flow_mass[0]
+               - 3600*m.fs.hxd.shell_inlet.flow_mass[0])
+        )
+
+    @m.fs.Constraint(doc="Max salt flow to hxd based on available hot salt")
+    def constraint_salt_maxflow_hot(b):
+        return (
+            3600*m.fs.hxd.shell_inlet.flow_mass[0] <=
+            m.fs.previous_salt_inventory_hot
+        )
+
+    @m.fs.Constraint(doc="Max salt flow to hxc based on available cold salt")
+    def constraint_salt_maxflow_cold(b):
+        return (
+            3600*m.fs.hxc.tube_inlet.flow_mass[0] <=
+            m.fs.previous_salt_inventory_cold
+        )
+
+    @m.fs.Constraint(doc="Maximum salt inventory at any time")
+    def constraint_salt_inventory(b):
+        return (
+            m.fs.salt_inventory_hot +
+            m.fs.salt_inventory_cold == m.fs.salt_amount)
+
+    return m
+
+
+def usc_unfix_dof(m):
     # Unfix data
-    m.usc_mp.fs.boiler.inlet.flow_mol[0].unfix()
+    m.fs.boiler.inlet.flow_mol[0].unfix()
 
     # Unfix storage system data
-    m.usc_mp.fs.ess_hp_split.split_fraction[0, "to_hxc"].unfix()
-    m.usc_mp.fs.ess_bfp_split.split_fraction[0, "to_hxd"].unfix()
-    for salt_hxc in [m.usc_mp.fs.hxc]:
+    m.fs.ess_hp_split.split_fraction[0, "to_hxc"].unfix()
+    m.fs.ess_bfp_split.split_fraction[0, "to_hxd"].unfix()
+    for salt_hxc in [m.fs.hxc]:
         salt_hxc.shell_inlet.unfix()
         salt_hxc.tube_inlet.flow_mass.unfix()  # kg/s, 1 DOF
         salt_hxc.area.unfix()  # 1 DOF
 
-    for salt_hxd in [m.usc_mp.fs.hxd]:
+    for salt_hxd in [m.fs.hxd]:
         salt_hxd.tube_inlet.unfix()
         salt_hxd.shell_inlet.flow_mass.unfix()  # kg/s, 1 DOF
         salt_hxd.area.unfix()  # 1 DOF
 
-    for unit in [m.usc_mp.fs.cooler]:
+    for unit in [m.fs.cooler]:
         unit.inlet.unfix()
-    m.usc_mp.fs.cooler.outlet.enth_mol[0].unfix()  # 1 DOF
+    m.fs.cooler.outlet.enth_mol[0].unfix()  # 1 DOF
 
     # Fix storage heat exchangers area and salt temperatures
-    m.usc_mp.fs.hxc.area.fix(1904)
-    m.usc_mp.fs.hxd.area.fix(2830)
-    m.usc_mp.fs.hxc.tube_outlet.temperature[0].fix(831)
-    m.usc_mp.fs.hxd.shell_inlet.temperature[0].fix(831)
-    m.usc_mp.fs.hxd.shell_outlet.temperature[0].fix(513.15)
+    m.fs.hxc.area.fix(1904)
+    m.fs.hxd.area.fix(2830)
+    m.fs.hxc.tube_outlet.temperature[0].fix(831)
+    m.fs.hxd.shell_inlet.temperature[0].fix(831)
+    m.fs.hxd.shell_outlet.temperature[0].fix(513.15)
+    
 
-    return m
 
+def usc_custom_init(m):
+
+    blk = usc.build_plant_model()
+    usc.initialize(blk)
+
+    # Create a flowsheet, add properties, unit models, and arcs
+    max_power = 436  # in MW
+    min_power = int(0.65 * 436)
+    pmin = int(0.65 * 436) + 1
+    pmax = 436 + 30
+    min_storage_heat_duty = 10  # in MW
+    max_storage_heat_duty = 200  # in MW
+
+
+    blk = usc_w_tes.create_integrated_model(blk, max_power=max_power)
+
+    # Give all the required inputs to the model
+    usc_w_tes.set_model_input(blk)
+
+    # Add scaling factor
+    usc_w_tes.set_scaling_factors(blk)
 
-def create_usc_mp_block(pmin=None, pmax=None):
-    print('>>> Creating USC model and initialization for each time period')
+    # Initialize the model with a sequential initialization and custom
+    # routines
+    usc_w_tes.initialize(blk)
+
+    # Add cost correlations
+    blk = usc_w_tes.build_costing(blk)
+
+    # Initialize with bounds
+    usc_w_tes.initialize_with_costing(blk)
+
+    blk.fs.plant_min_power_eq = Constraint(
+        expr=blk.fs.plant_power_out[0] >= min_power
+    )
+    blk.fs.plant_max_power_eq = Constraint(
+        expr=blk.fs.plant_power_out[0] <= max_power
+    )
+
+    blk.fs.hxc.heat_duty.setlb(min_storage_heat_duty * 1e6)
+    blk.fs.hxd.heat_duty.setlb(min_storage_heat_duty * 1e6)
+
+    blk.fs.hxc.heat_duty.setub(max_storage_heat_duty * 1e6)
+    blk.fs.hxd.heat_duty.setub(max_storage_heat_duty * 1e6)
 
     if pmin is None:
         pmin = int(0.65 * 436) + 1
     if pmax is None:
         pmax = 436 + 30
 
-    m = create_usc_model(pmin, pmax)
-    b1 = m.usc_mp
-
     # Add coupling variables
-    b1.previous_power = Var(
+    blk.fs.previous_power = Var(
         domain=NonNegativeReals,
         initialize=300,
         bounds=(pmin, pmax),
@@ -121,73 +258,76 @@ def create_usc_mp_block(pmin=None, pmax=None):
     inventory_min = 75000
     tank_max = 6739292           # Units in kg
 
-    b1.previous_salt_inventory_hot = Var(
+    blk.fs.previous_salt_inventory_hot = Var(
         domain=NonNegativeReals,
         initialize=inventory_min,
         bounds=(0, inventory_max),
         doc="Hot salt at the beginning of the hour (or time period), kg"
         )
-    b1.salt_inventory_hot = Var(
+    blk.fs.salt_inventory_hot = Var(
         domain=NonNegativeReals,
         initialize=inventory_min,
         bounds=(0, inventory_max),
         doc="Hot salt inventory at the end of the hour (or time period), kg"
         )
-    b1.previous_salt_inventory_cold = Var(
+    blk.fs.previous_salt_inventory_cold = Var(
         domain=NonNegativeReals,
         initialize=tank_max-inventory_min,
         bounds=(0, inventory_max),
         doc="Cold salt at the beginning of the hour (or time period), kg"
         )
-    b1.salt_inventory_cold = Var(
+    blk.fs.salt_inventory_cold = Var(
         domain=NonNegativeReals,
         initialize=tank_max-inventory_min,
         bounds=(0, inventory_max),
         doc="Cold salt inventory at the end of the hour (or time period), kg"
         )
 
-    @b1.fs.Constraint(doc="Plant ramping down constraint")
+    @blk.fs.Constraint(doc="Plant ramping down constraint")
     def constraint_ramp_down(b):
         return (
-            b1.previous_power - 60 <=
-            b1.fs.plant_power_out[0])
+            blk.fs.previous_power - 60 <=
+            blk.fs.plant_power_out[0])
 
-    @b1.fs.Constraint(doc="Plant ramping up constraint")
+    @blk.fs.Constraint(doc="Plant ramping up constraint")
     def constraint_ramp_up(b):
         return (
-            b1.previous_power + 60 >=
-            b1.fs.plant_power_out[0])
+            blk.fs.previous_power + 60 >=
+            blk.fs.plant_power_out[0])
 
-    @b1.fs.Constraint(doc="Inventory balance at the end of the time period")
+    @blk.fs.Constraint(doc="Inventory balance at the end of the time period")
     def constraint_salt_inventory_hot(b):
         return (
-            b1.salt_inventory_hot ==
-            b1.previous_salt_inventory_hot
-            + (3600*b1.fs.hxc.tube_inlet.flow_mass[0]
-               - 3600*b1.fs.hxd.shell_inlet.flow_mass[0])
+            blk.fs.salt_inventory_hot ==
+            blk.fs.previous_salt_inventory_hot
+            + (3600*blk.fs.hxc.tube_inlet.flow_mass[0]
+               - 3600*blk.fs.hxd.shell_inlet.flow_mass[0])
         )
 
-    @b1.fs.Constraint(doc="Max salt flow to hxd based on available hot salt")
+    @blk.fs.Constraint(doc="Max salt flow to hxd based on available hot salt")
     def constraint_salt_maxflow_hot(b):
         return (
-            3600*b1.fs.hxd.shell_inlet.flow_mass[0] <=
-            b1.previous_salt_inventory_hot
+            3600*blk.fs.hxd.shell_inlet.flow_mass[0] <=
+            blk.fs.previous_salt_inventory_hot
         )
 
-    @b1.fs.Constraint(doc="Max salt flow to hxc based on available cold salt")
+    @blk.fs.Constraint(doc="Max salt flow to hxc based on available cold salt")
     def constraint_salt_maxflow_cold(b):
         return (
-            3600*b1.fs.hxc.tube_inlet.flow_mass[0] <=
-            b1.previous_salt_inventory_cold
+            3600*blk.fs.hxc.tube_inlet.flow_mass[0] <=
+            blk.fs.previous_salt_inventory_cold
         )
 
-    @b1.fs.Constraint(doc="Maximum salt inventory at any time")
+    @blk.fs.Constraint(doc="Maximum salt inventory at any time")
     def constraint_salt_inventory(b):
         return (
-            b1.salt_inventory_hot +
-            b1.salt_inventory_cold == b1.fs.salt_amount)
+            blk.fs.salt_inventory_hot +
+            blk.fs.salt_inventory_cold == blk.fs.salt_amount)
 
-    return m
+    init_model = to_json(blk, return_dict=True)
+    from_json(m, sd=init_model)
+
+    return
 
 
 # The tank level and power output are linked between the contiguous time periods
@@ -196,10 +336,10 @@ def get_usc_link_variable_pairs(b1, b2):
         b1: current time block
         b2: next time block
     """
-    return [(b1.usc_mp.salt_inventory_hot,
-             b2.usc_mp.previous_salt_inventory_hot),
-            (b1.usc_mp.fs.plant_power_out[0],
-             b2.usc_mp.previous_power)]
+    return [(b1.fs.salt_inventory_hot,
+             b2.fs.previous_salt_inventory_hot),
+            (b1.fs.plant_power_out[0],
+             b2.fs.previous_power)]
 
 
 # The tank level at the end of the last time period must be the same as at the
@@ -210,8 +350,8 @@ def get_usc_periodic_variable_pairs(b1, b2):
         b2: first time block
     """
     # return
-    return [(b1.usc_mp.salt_inventory_hot,
-             b2.usc_mp.previous_salt_inventory_hot)]
+    return [(b1.fs.salt_inventory_hot,
+             b2.fs.previous_salt_inventory_hot)]
 
 # Create the multiperiod model object. You can pass arguments to your
 # "process_model_func" for each time period using a dict of dicts as
@@ -227,13 +367,15 @@ def create_multiperiod_usc_model(n_time_points=4, pmin=None, pmax=None):
     n_time_points: Number of time blocks to create
     """
     multiperiod_usc = MultiPeriodModel(
-        n_time_points,
-        lambda: create_usc_mp_block(pmin=None, pmax=None),
-        get_usc_link_variable_pairs,
-        get_usc_periodic_variable_pairs
+        n_time_points=n_time_points,
+        process_model_func=create_usc_model,
+        initialization_func=usc_custom_init,
+        unfix_dof_func=usc_unfix_dof,
+        linking_variable_func=get_usc_link_variable_pairs,
+        flowsheet_options={"pmin": pmin,
+                           "pmax": pmax},
+        use_stochastic_build=True,
+        outlvl=idaeslog.INFO,
         )
 
-    # If you have no arguments, you don't actually need to pass in
-    # anything. NOTE: building the model will initialize each time block
-    multiperiod_usc.build_multi_period_model()
     return multiperiod_usc