diff --git a/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/genx_settings.yml b/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/genx_settings.yml
index 044864cf44..19d3dd63d7 100644
--- a/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/genx_settings.yml
+++ b/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/genx_settings.yml
@@ -1,8 +1,8 @@
 OverwriteResults: 0 # Overwrite existing results in output folder or create a new one; 0 = create new folder; 1 = overwrite existing results
 PrintModel: 0 # Write the model formulation as an output; 0 = active; 1 = not active
-NetworkExpansion: 0 # Transmission network expansionl; 0 = not active; 1 = active systemwide
+NetworkExpansion: 1 # Transmission network expansionl; 0 = not active; 1 = active systemwide
 Trans_Loss_Segments: 1 # Number of segments used in piecewise linear approximation of transmission losses; 1 = linear, >2 = piecewise quadratic
-Reserves: 0 # Regulation (primary) and operating (secondary) reserves; 0 = not active, 1 = active systemwide
+Reserves: 1 # Regulation (primary) and operating (secondary) reserves; 0 = not active, 1 = active systemwide
 EnergyShareRequirement: 1 # Minimum qualifying renewables penetration; 0 = not active; 1 = active systemwide
 CapacityReserveMargin: 1 # Number of capacity reserve margin constraints; 0 = not active; 1 = active systemwide
 CO2Cap: 1 # CO2 emissions cap; 0 = not active (no CO2 emission limit); 1 = mass-based emission limit constraint; 2 = demand + rate-based emission limit constraint; 3 = generation + rate-based emission limit constraint
@@ -18,3 +18,4 @@ ModelingToGenerateAlternatives: 0 # Modeling to generate alternatives; 0 = not a
 ModelingtoGenerateAlternativeSlack: 0.1 # Slack value as a fraction of least-cost objective in budget constraint used for evaluating alternative model solutions; positive float value
 ModelingToGenerateAlternativeIterations: 3 # Number of MGA iterations with maximization and minimization objective
 MethodofMorris: 0 #Flag for turning on the Method of Morris analysis
+WriteOutputs: "annual"
diff --git a/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/output_settings.yml b/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/output_settings.yml
new file mode 100644
index 0000000000..283e141ca7
--- /dev/null
+++ b/Example_Systems/RealSystemExample/ISONE_Singlezone/Settings/output_settings.yml
@@ -0,0 +1,27 @@
+WriteCosts: false
+WriteCapacity: false
+WriteCapacityValue: false
+WriteCapacityFactor: false
+WriteCharge: false
+WriteChargingCost: false
+WriteCO2: false
+WriteCO2Cap: false
+WriteCommit: false
+WriteCurtailment: false
+WriteEmissions: false
+WriteEnergyRevenue: false
+WriteESRPrices: false
+WriteESRRevenue: false
+WriteFuelConsumption: false
+WriteHourlyMatchingPrices: false
+WriteHydrogenPrices: false
+WriteMaintenance: false
+WriteMaxCapReq: false
+WriteMinCapReq: false
+WriteNetRevenue: false
+WriteNSE: false
+WriteNWExpansion: false
+WriteTransmissionFlows: false
+WriteTransmissionLosses: false
+WriteVirtualDischarge: false
+WriteVREStor: false
\ No newline at end of file
diff --git a/docs/src/data_documentation.md b/docs/src/data_documentation.md
index 77b5ff5051..487085e16c 100644
--- a/docs/src/data_documentation.md
+++ b/docs/src/data_documentation.md
@@ -70,6 +70,9 @@ Note that all settings parameters are case sensitive.
 |PrintModel | Flag for printing the model equations as .lp file.|
 ||1 = including the model equation as an output|
 ||0 = the model equation won't be included as an output|
+| WriteOutputs | Flag for writing the model outputs with hourly resolution or just the annual sum.|
+|| "full" = write the model outputs with hourly resolution.|
+|| "annual" = write only the annual sum of the model outputs.|
 
 Additionally, Solver related settings parameters are specified in the appropriate .yml file (e.g. `gurobi_settings.yml` or `cplex_settings.yml`),
 which should be located in the current working directory.
@@ -77,6 +80,8 @@ Note that GenX supplies default settings for most solver settings in the various
 To overwrite default settings, you can specify the below Solver specific settings.
 Settings are specific to each solver.
 
+(Optional) The user can also select the output files that they want to export using the `output_settings.yml` file. This file containes a list of `yes/no` options for each output file, and should be located in the `Settings` folder. By default, if `output_settings.yml` is not included, GenX will export all output files. 
+
 ###### Table 1b: Summary of the Solver settings parameters
 ---
 |**Settings Parameter** | **Description**|
diff --git a/src/case_runners/case_runner.jl b/src/case_runners/case_runner.jl
index 9ac5431f82..f0dc105591 100644
--- a/src/case_runners/case_runner.jl
+++ b/src/case_runners/case_runner.jl
@@ -14,8 +14,9 @@ end
 case - folder for the case
 """
 function run_genx_case!(case::AbstractString, optimizer::Any=HiGHS.Optimizer)
-    genx_settings = get_settings_path(case, "genx_settings.yml") #Settings YAML file path
-    mysetup = configure_settings(genx_settings) # mysetup dictionary stores settings and GenX-specific parameters
+    genx_settings = get_settings_path(case, "genx_settings.yml") # Settings YAML file path
+    writeoutput_settings = get_settings_path(case, "output_settings.yml") # Write-output settings YAML file path
+    mysetup = configure_settings(genx_settings, writeoutput_settings) # mysetup dictionary stores settings and GenX-specific parameters
 
     if mysetup["MultiStage"] == 0
         run_genx_case_simple!(case, mysetup, optimizer)
diff --git a/src/configure_settings/configure_settings.jl b/src/configure_settings/configure_settings.jl
index 92074cada0..f93709ce95 100644
--- a/src/configure_settings/configure_settings.jl
+++ b/src/configure_settings/configure_settings.jl
@@ -24,19 +24,23 @@ function default_settings()
         "IncludeLossesInESR" => 0,
         "HydrogenHourlyMatching" => 0,
         "EnableJuMPStringNames" => false,
+        "HydrogenHourlyMatching" => 0,
+        "WriteOutputs" => "full",
         "ComputeConflicts" => 0,
         "ResourcePath" => "Resources",
     )
 end
 
-function configure_settings(settings_path::String)
+function configure_settings(settings_path::String, output_settings_path::String)
     println("Configuring Settings")
     model_settings = YAML.load(open(settings_path))
 
     settings = default_settings()
-
     merge!(settings, model_settings)
 
+    output_settings = configure_writeoutput(output_settings_path, settings)
+    settings["WriteOutputsSettingsDict"] = output_settings
+
     validate_settings!(settings)
     return settings
 end
@@ -44,6 +48,10 @@ end
 function validate_settings!(settings::Dict{Any,Any})
     # Check for any settings combinations that are not allowed.
     # If we find any then make a response and issue a note to the user.
+    
+    # make WriteOutputs setting lowercase and check for valid value
+    settings["WriteOutputs"] = lowercase(settings["WriteOutputs"])
+    @assert settings["WriteOutputs"] ∈ ["annual", "full"]
 
     if "OperationWrapping" in keys(settings)
         @warn """The behavior of the TimeDomainReduction and OperationWrapping
@@ -57,3 +65,79 @@ function validate_settings!(settings::Dict{Any,Any})
     end
 
 end
+
+function default_writeoutput()
+    Dict{String,Bool}(
+        "WriteCosts" => true,
+        "WriteCapacity" => true,
+        "WriteCapacityValue" => true,
+        "WriteCapacityFactor" => true,
+        "WriteCharge" => true,
+        "WriteChargingCost" => true,
+        "WriteCO2" => true,
+        "WriteCO2Cap" => true,
+        "WriteCommit" => true,
+        "WriteCurtailment" => true,
+        "WriteEmissions" => true,
+        "WriteEnergyRevenue" => true,
+        "WriteESRPrices" => true,
+        "WriteESRRevenue" => true,
+        "WriteFuelConsumption" => true,
+        "WriteHourlyMatchingPrices" => true,
+        "WriteHydrogenPrices" => true,
+        "WriteMaintenance" => true,
+        "WriteMaxCapReq" => true,
+        "WriteMinCapReq" => true,
+        "WriteNetRevenue" => true,
+        "WriteNSE" => true,
+        "WriteNWExpansion" => true,
+        "WriteOpWrapLDSdStor" => true,
+        "WriteOpWrapLDSStorInit" => true,
+        "WritePower" => true,
+        "WritePowerBalance" => true,
+        "WritePrice" => true,
+        "WriteReg" => true,
+        "WriteReliability" => true,
+        "WriteReserveMargin" => true,
+        "WriteReserveMarginRevenue" => true,
+        "WriteReserveMarginSlack" => true,
+        "WriteReserveMarginWithWeights" => true,
+        "WriteRsv" => true,
+        "WriteShutdown" => true,
+        "WriteStart" => true,
+        "WriteStatus" => true,
+        "WriteStorage" => true,
+        "WriteStorageDual" => true,
+        "WriteSubsidyRevenue" => true,
+        "WriteTimeWeights" => true,
+        "WriteTransmissionFlows" => true,
+        "WriteTransmissionLosses" => true,
+        "WriteVirtualDischarge" => true,
+        "WriteVREStor" => true
+    )
+end
+
+function configure_writeoutput(output_settings_path::String, settings::Dict)
+    
+    writeoutput = default_writeoutput()
+
+    # don't write files with hourly data if settings["WriteOutputs"] == "annual"
+    if settings["WriteOutputs"] == "annual"
+        writeoutput["WritePrice"] = false
+        writeoutput["WriteReliability"] = false
+        writeoutput["WriteStorage"] = false
+        writeoutput["WriteStorageDual"] = false
+        writeoutput["WriteTimeWeights"] = false
+        writeoutput["WriteCommit"] = false
+        writeoutput["WriteCapacityValue"] = false
+        writeoutput["WriteReserveMargin"] = false
+        writeoutput["WriteReserveMarginWithWeights"] = false
+    end
+
+    # read in YAML file if provided
+    if isfile(output_settings_path)
+        model_writeoutput = YAML.load(open(output_settings_path))
+        merge!(writeoutput, model_writeoutput)
+    end
+    return writeoutput
+end
\ No newline at end of file
diff --git a/src/multi_stage/dual_dynamic_programming.jl b/src/multi_stage/dual_dynamic_programming.jl
index b3e2c23c14..b9ed9565f0 100644
--- a/src/multi_stage/dual_dynamic_programming.jl
+++ b/src/multi_stage/dual_dynamic_programming.jl
@@ -336,7 +336,7 @@ function write_multi_stage_outputs(stats_d::Dict, outpath::String, settings_d::D
     	write_multi_stage_network_expansion(outpath, multi_stage_settings_d)
     end
     write_multi_stage_costs(outpath, multi_stage_settings_d, inputs_dict)
-    write_multi_stage_stats(outpath, stats_d)
+    multi_stage_settings_d["Myopic"] == 0 && write_multi_stage_stats(outpath, stats_d)
     write_multi_stage_settings(outpath, settings_d)
 
 end
diff --git a/src/write_outputs/capacity_reserve_margin/write_reserve_margin.jl b/src/write_outputs/capacity_reserve_margin/write_reserve_margin.jl
index e889425c89..6d2f8c2e80 100644
--- a/src/write_outputs/capacity_reserve_margin/write_reserve_margin.jl
+++ b/src/write_outputs/capacity_reserve_margin/write_reserve_margin.jl
@@ -5,5 +5,5 @@ function write_reserve_margin(path::AbstractString, setup::Dict, EP::Model)
 	end
 	dfResMar = DataFrame(temp_ResMar, :auto)
 	CSV.write(joinpath(path, "ReserveMargin.csv"), dfResMar)
-	return dfResMar
+	return nothing
 end
diff --git a/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl b/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
index 221241c278..99b0e9e0f6 100644
--- a/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
+++ b/src/write_outputs/capacity_reserve_margin/write_reserve_margin_slack.jl
@@ -4,14 +4,22 @@ function write_reserve_margin_slack(path::AbstractString, inputs::Dict, setup::D
     dfResMar_slack = DataFrame(CRM_Constraint = [Symbol("CapRes_$res") for res = 1:NCRM], 
                                 AnnualSum = value.(EP[:eCapResSlack_Year]),
                                 Penalty = value.(EP[:eCCapResSlack]))
-    temp_ResMar_slack = value.(EP[:vCapResSlack])
+   
     if setup["ParameterScale"] == 1
         dfResMar_slack.AnnualSum .*= ModelScalingFactor # Convert GW to MW
         dfResMar_slack.Penalty .*= ModelScalingFactor^2 # Convert Million $ to $
-        temp_ResMar_slack .*= ModelScalingFactor # Convert GW to MW
     end
-    dfResMar_slack = hcat(dfResMar_slack, DataFrame(temp_ResMar_slack, [Symbol("t$t") for t in 1:T]))
-    CSV.write(joinpath(path, "ReserveMargin_prices_and_penalties.csv"), dftranspose(dfResMar_slack, false), header=false)
-    return dfResMar_slack
+
+    if setup["WriteOutputs"] == "annual"
+        CSV.write(joinpath(path, "ReserveMargin_prices_and_penalties.csv"), dfResMar_slack)
+    else     # setup["WriteOutputs"] == "full"
+        temp_ResMar_slack = value.(EP[:vCapResSlack])
+        if setup["ParameterScale"] == 1
+            temp_ResMar_slack .*= ModelScalingFactor # Convert GW to MW
+        end
+        dfResMar_slack = hcat(dfResMar_slack, DataFrame(temp_ResMar_slack, [Symbol("t$t") for t in 1:T]))
+        CSV.write(joinpath(path, "ReserveMargin_prices_and_penalties.csv"), dftranspose(dfResMar_slack, false), writeheader=false)
+    end
+    return nothing
 end
 
diff --git a/src/write_outputs/capacity_reserve_margin/write_virtual_discharge.jl b/src/write_outputs/capacity_reserve_margin/write_virtual_discharge.jl
index b093e984d9..9a4d8308a6 100644
--- a/src/write_outputs/capacity_reserve_margin/write_virtual_discharge.jl
+++ b/src/write_outputs/capacity_reserve_margin/write_virtual_discharge.jl
@@ -5,30 +5,25 @@ Function for writing the "virtual" discharge of each storage technology. Virtual
 	allow storage resources to contribute to the capacity reserve margin without actually discharging.
 """
 function write_virtual_discharge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
 
 	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
 	T = inputs["T"]     # Number of time steps (hours)
 	STOR_ALL = inputs["STOR_ALL"]
 
-	dfVirtualDischarge = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones, AnnualSum = Array{Union{Missing,Float64}}(undef, G))
-	virtual_discharge = zeros(G,T)
-
 	scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-	if !isempty(STOR_ALL)
-	    virtual_discharge[STOR_ALL, :] = (value.(EP[:vCAPRES_discharge][STOR_ALL, :]).data - value.(EP[:vCAPRES_charge][STOR_ALL, :]).data) * scale_factor
-	end
+	
+	resources = inputs["RESOURCE_NAMES"][STOR_ALL]
+	zones = inputs["R_ZONES"][STOR_ALL]
+	virtual_discharge = (value.(EP[:vCAPRES_discharge][STOR_ALL, :].data) - value.(EP[:vCAPRES_charge][STOR_ALL, :].data)) * scale_factor
 
+	dfVirtualDischarge = DataFrame(Resource = resources, Zone = zones)
 	dfVirtualDischarge.AnnualSum .= virtual_discharge * inputs["omega"]
-	dfVirtualDischarge = hcat(dfVirtualDischarge, DataFrame(virtual_discharge, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfVirtualDischarge,auxNew_Names)
-	total = DataFrame(["Total" 0 sum(dfVirtualDischarge[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
 
-	total[:, 4:T+3] .= sum(virtual_discharge, dims = 1)
-	rename!(total,auxNew_Names)
-	dfVirtualDischarge = vcat(dfVirtualDischarge, total)
-	CSV.write(joinpath(path, "virtual_discharge.csv"), dftranspose(dfVirtualDischarge, false), header=false)
-	return dfVirtualDischarge
-end
+	filepath = joinpath(path, "virtual_discharge.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfVirtualDischarge)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filepath, virtual_discharge, dfVirtualDischarge)
+	end
+	return nothing
+end 
diff --git a/src/write_outputs/co2_cap/write_co2_cap.jl b/src/write_outputs/co2_cap/write_co2_cap.jl
index fbe5e28d4f..fa8e479ec8 100644
--- a/src/write_outputs/co2_cap/write_co2_cap.jl
+++ b/src/write_outputs/co2_cap/write_co2_cap.jl
@@ -21,5 +21,5 @@ function write_co2_cap(path::AbstractString, inputs::Dict, setup::Dict, EP::Mode
 
     CSV.write(joinpath(path, "CO2_prices_and_penalties.csv"), dfCO2Price)
 
-    return dfCO2Price
+    return nothing
 end
diff --git a/src/write_outputs/hydrogen/write_hourly_matching_prices.jl b/src/write_outputs/hydrogen/write_hourly_matching_prices.jl
index c804d810c4..393544e4e4 100644
--- a/src/write_outputs/hydrogen/write_hourly_matching_prices.jl
+++ b/src/write_outputs/hydrogen/write_hourly_matching_prices.jl
@@ -13,5 +13,5 @@ function write_hourly_matching_prices(path::AbstractString, inputs::Dict, setup:
 
 	CSV.write(joinpath(path, "hourly_matching_prices.csv"), dftranspose(dfHourlyMatchPrices, false), header=false)
 
-	return dfHourlyMatchPrices
+	return nothing
 end
diff --git a/src/write_outputs/hydrogen/write_hydrogen_prices.jl b/src/write_outputs/hydrogen/write_hydrogen_prices.jl
index 7c6948844b..1d7d905491 100644
--- a/src/write_outputs/hydrogen/write_hydrogen_prices.jl
+++ b/src/write_outputs/hydrogen/write_hydrogen_prices.jl
@@ -3,5 +3,5 @@ function write_hydrogen_prices(path::AbstractString, inputs::Dict, setup::Dict,
 	dfHydrogenPrice = DataFrame(Hydrogen_Price_Per_Tonne = convert(Array{Float64}, dual.(EP[:cHydrogenMin])*scale_factor))
 
 	CSV.write(joinpath(path, "hydrogen_prices.csv"), dfHydrogenPrice)
-	return dfHydrogenPrice
+	return nothing
 end
diff --git a/src/write_outputs/reserves/write_reg.jl b/src/write_outputs/reserves/write_reg.jl
index 8d8d51f217..4b984fcc14 100644
--- a/src/write_outputs/reserves/write_reg.jl
+++ b/src/write_outputs/reserves/write_reg.jl
@@ -1,23 +1,20 @@
 function write_reg(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
-
-	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-	T = inputs["T"]     # Number of time steps (hours)
 	REG = inputs["REG"]
 	scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
+
+	resources = inputs["RESOURCE_NAMES"][REG]
+	zones = inputs["R_ZONES"][REG]
 	# Regulation contributions for each resource in each time step
-	dfReg = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
-	reg = zeros(G,T)
-	reg[REG, :] = value.(EP[:vREG][REG, :])
-	dfReg.AnnualSum = (reg*scale_factor) * inputs["omega"]
-	dfReg = hcat(dfReg, DataFrame(reg*scale_factor, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfReg,auxNew_Names)
+	reg = value.(EP[:vREG][REG, :].data) * scale_factor
+
+	dfReg = DataFrame(Resource = resources, Zone = zones)
+	dfReg.AnnualSum = reg * inputs["omega"]
 
-	total = DataFrame(["Total" 0 sum(dfReg.AnnualSum) fill(0.0, (1,T))], :auto)
-	total[!, 4:T+3] .= sum(reg, dims = 1)
-	rename!(total,auxNew_Names)
-	dfReg = vcat(dfReg, total)
-	CSV.write(joinpath(path, "reg.csv"), dftranspose(dfReg, false), header=false)
+	filepath = joinpath(path, "reg.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfReg)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filepath, reg, dfReg)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/reserves/write_rsv.jl b/src/write_outputs/reserves/write_rsv.jl
index 7db756e394..5eb853db8c 100644
--- a/src/write_outputs/reserves/write_rsv.jl
+++ b/src/write_outputs/reserves/write_rsv.jl
@@ -1,28 +1,31 @@
 function write_rsv(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
-
-	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-	T = inputs["T"]     # Number of time steps (hours)
 	RSV = inputs["RSV"]
 	scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
-	dfRsv = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
-	rsv = zeros(G,T)
-	unmet_vec = zeros(T)
-	rsv[RSV, :] = value.(EP[:vRSV][RSV, :]) * scale_factor
-	unmet_vec = value.(EP[:vUNMET_RSV]) * scale_factor
-	total_unmet = sum(unmet_vec)
+
+	resources = inputs["RESOURCE_NAMES"][RSV]
+	zones = inputs["R_ZONES"][RSV]
+	rsv = value.(EP[:vRSV][RSV, :].data) * scale_factor
+
+	dfRsv = DataFrame(Resource = resources, Zone = zones)
+
 	dfRsv.AnnualSum = rsv * inputs["omega"]
-	dfRsv = hcat(dfRsv, DataFrame(rsv, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfRsv,auxNew_Names)
 
-	total = DataFrame(["Total" 0 sum(dfRsv.AnnualSum) zeros(1, T)], :auto)
-	unmet = DataFrame(["unmet" 0 total_unmet zeros(1, T)], :auto)
-	total[!, 4:T+3] .= sum(rsv, dims = 1)
-	unmet[!, 4:T+3] .= transpose(unmet_vec)
-	rename!(total,auxNew_Names)
-	rename!(unmet,auxNew_Names)
-	dfRsv = vcat(dfRsv, unmet, total)
-	CSV.write(joinpath(path, "reg_dn.csv"), dftranspose(dfRsv, false), header=false)
+	if setup["WriteOutputs"] == "annual"
+		write_annual(joinpath(path, "reg_dn.csv"), dfRsv)
+	else 	# setup["WriteOutputs"] == "full"
+		unmet_vec = value.(EP[:vUNMET_RSV]) * scale_factor
+		total_unmet = sum(unmet_vec)
+		dfRsv = hcat(dfRsv, DataFrame(rsv, :auto))
+		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+		rename!(dfRsv,auxNew_Names)
+		
+		total = DataFrame(["Total" 0 sum(dfRsv.AnnualSum) zeros(1, T)], :auto)
+		unmet = DataFrame(["unmet" 0 total_unmet zeros(1, T)], :auto)
+		total[!, 4:T+3] .= sum(rsv, dims = 1)
+		unmet[!, 4:T+3] .= transpose(unmet_vec)
+		rename!(total,auxNew_Names)
+		rename!(unmet,auxNew_Names)
+		dfRsv = vcat(dfRsv, unmet, total)
+		CSV.write(joinpath(path, "reg_dn.csv"), dftranspose(dfRsv, false), writeheader=false)
+	end
 end
diff --git a/src/write_outputs/transmission/write_transmission_flows.jl b/src/write_outputs/transmission/write_transmission_flows.jl
index 28a6fddffa..008c3ac263 100644
--- a/src/write_outputs/transmission/write_transmission_flows.jl
+++ b/src/write_outputs/transmission/write_transmission_flows.jl
@@ -8,13 +8,22 @@ function write_transmission_flows(path::AbstractString, inputs::Dict, setup::Dic
 	if setup["ParameterScale"] == 1
 	    flow *= ModelScalingFactor
 	end
+
 	dfFlow.AnnualSum = flow * inputs["omega"]
-	dfFlow = hcat(dfFlow, DataFrame(flow, :auto))
-	auxNew_Names=[Symbol("Line");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfFlow,auxNew_Names)
-	total = DataFrame(["Total" sum(dfFlow.AnnualSum) fill(0.0, (1,T))], :auto)
-	total[:, 3:T+2] .= sum(flow, dims = 1)
-	rename!(total,auxNew_Names)
-	dfFlow = vcat(dfFlow, total)
-	CSV.write(joinpath(path, "flow.csv"), dftranspose(dfFlow, false), header=false)
+
+	filepath = joinpath(path, "flow.csv")	
+	if setup["WriteOutputs"] == "annual"
+		total = DataFrame(["Total" sum(dfFlow.AnnualSum)], [:Line, :AnnualSum])
+		dfFlow = vcat(dfFlow, total)
+		CSV.write(filepath, dfFlow)
+	else 	# setup["WriteOutputs"] == "full" 
+		dfFlow = hcat(dfFlow, DataFrame(flow, :auto))
+		auxNew_Names=[Symbol("Line");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+		rename!(dfFlow,auxNew_Names)
+		total = DataFrame(["Total" sum(dfFlow.AnnualSum) fill(0.0, (1,T))], auxNew_Names)
+		total[:, 3:T+2] .= sum(flow, dims = 1)
+		dfFlow = vcat(dfFlow, total)
+		CSV.write(filepath, dftranspose(dfFlow, false), writeheader=false)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/transmission/write_transmission_losses.jl b/src/write_outputs/transmission/write_transmission_losses.jl
index 583e00dc9f..8f5bb51977 100644
--- a/src/write_outputs/transmission/write_transmission_losses.jl
+++ b/src/write_outputs/transmission/write_transmission_losses.jl
@@ -9,14 +9,21 @@ function write_transmission_losses(path::AbstractString, inputs::Dict, setup::Di
 	if setup["ParameterScale"] == 1
 	    tlosses[LOSS_LINES, :] *= ModelScalingFactor
 	end
-	dfTLosses.AnnualSum = tlosses * inputs["omega"]
-	dfTLosses = hcat(dfTLosses, DataFrame(tlosses, :auto))
-	auxNew_Names=[Symbol("Line");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfTLosses,auxNew_Names)
 
-	total = DataFrame(["Total" sum(dfTLosses.AnnualSum) fill(0.0, (1,T))], :auto)
-	total[:, 3:T+2] .= sum(tlosses, dims = 1)
-	rename!(total,auxNew_Names)
-	dfTLosses = vcat(dfTLosses, total)
-	CSV.write(joinpath(path, "tlosses.csv"), dftranspose(dfTLosses, false), header=false)
+	dfTLosses.AnnualSum = tlosses * inputs["omega"]
+	
+	if setup["WriteOutputs"] == "annual"
+		total = DataFrame(["Total" sum(dfTLosses.AnnualSum)], [:Line, :AnnualSum])
+		dfTLosses = vcat(dfTLosses, total)
+		CSV.write(joinpath(path, "tlosses.csv"), dfTLosses)
+	else
+		dfTLosses = hcat(dfTLosses, DataFrame(tlosses, :auto))
+		auxNew_Names=[Symbol("Line");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+		rename!(dfTLosses,auxNew_Names)
+		total = DataFrame(["Total" sum(dfTLosses.AnnualSum) fill(0.0, (1,T))], auxNew_Names)
+		total[:, 3:T+2] .= sum(tlosses, dims = 1)
+		dfTLosses = vcat(dfTLosses, total)
+		CSV.write(joinpath(path, "tlosses.csv"), dftranspose(dfTLosses, false), writeheader=false)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/ucommit/write_commit.jl b/src/write_outputs/ucommit/write_commit.jl
index 3e4eff5387..685ad53e0a 100644
--- a/src/write_outputs/ucommit/write_commit.jl
+++ b/src/write_outputs/ucommit/write_commit.jl
@@ -1,14 +1,13 @@
 function write_commit(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
-
-	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-	T = inputs["T"]     # Number of time steps (hours)
+	
 	COMMIT = inputs["COMMIT"]
+	T = inputs["T"]
+
 	# Commitment state for each resource in each time step
-	commit = zeros(G,T)
-	commit[COMMIT, :] = value.(EP[:vCOMMIT][COMMIT, :])
-	dfCommit = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
+	resources = inputs["RESOURCE_NAMES"][COMMIT]
+	zones = inputs["R_ZONES"][COMMIT]
+  commit = value.(EP[:vCOMMIT][COMMIT, :].data)
+	dfCommit = DataFrame(Resource = resources, Zone = zones)
 	dfCommit = hcat(dfCommit, DataFrame(commit, :auto))
 	auxNew_Names=[Symbol("Resource");Symbol("Zone");[Symbol("t$t") for t in 1:T]]
 	rename!(dfCommit,auxNew_Names)
diff --git a/src/write_outputs/ucommit/write_shutdown.jl b/src/write_outputs/ucommit/write_shutdown.jl
index f524fb4daf..56325b25f6 100644
--- a/src/write_outputs/ucommit/write_shutdown.jl
+++ b/src/write_outputs/ucommit/write_shutdown.jl
@@ -1,22 +1,19 @@
 function write_shutdown(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
-
-	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-	T = inputs["T"]     # Number of time steps (hours)
 	# Operational decision variable states
-	COMMIT = inputs["COMMIT"]
+  COMMIT = inputs["COMMIT"]
+  zones = inputs["R_ZONES"][COMMIT]
 	# Shutdown state for each resource in each time step
-	dfShutdown = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
-	shut = zeros(G,T)
-	shut[COMMIT, :] = value.(EP[:vSHUT][COMMIT, :])
+	shut = value.(EP[:vSHUT][COMMIT, :].data)
+	resources = inputs["RESOURCE_NAMES"][COMMIT]
+
+	dfShutdown = DataFrame(Resource = resources, Zone = zones)
 	dfShutdown.AnnualSum = shut * inputs["omega"]
-	dfShutdown = hcat(dfShutdown, DataFrame(shut, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfShutdown,auxNew_Names)
-	total=DataFrame(["Total" 0 sum(dfShutdown.AnnualSum) fill(0.0, (1,T))], :auto)
-	total[:, 4:T+3] .= sum(shut, dims = 1)
-	rename!(total,auxNew_Names)
-	dfShutdown = vcat(dfShutdown, total)
-	CSV.write(joinpath(path, "shutdown.csv"), dftranspose(dfShutdown, false), header=false)
+
+	filepath = joinpath(path, "shutdown.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfShutdown)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filepath, shut, dfShutdown)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/ucommit/write_start.jl b/src/write_outputs/ucommit/write_start.jl
index 71aa7b723e..461d522a17 100644
--- a/src/write_outputs/ucommit/write_start.jl
+++ b/src/write_outputs/ucommit/write_start.jl
@@ -1,22 +1,19 @@
 function write_start(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
-	gen = inputs["RESOURCES"]
-	zones = zone_id.(gen)
-	
-	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-	T = inputs["T"]     # Number of time steps (hours)
+
 	COMMIT = inputs["COMMIT"]
 	# Startup state for each resource in each time step
-	dfStart = DataFrame(Resource = inputs["RESOURCE_NAMES"], Zone = zones)
-	start = zeros(G,T)
-	start[COMMIT, :] = value.(EP[:vSTART][COMMIT, :])
+	resources = inputs["RESOURCE_NAMES"][COMMIT]
+	zones = inputs["R_ZONES"][COMMIT]
+
+	dfStart = DataFrame(Resource = resources, Zone = zones)
+  start = value.(EP[:vSTART][COMMIT, :].data)
 	dfStart.AnnualSum = start * inputs["omega"]
-	dfStart = hcat(dfStart, DataFrame(start, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfStart,auxNew_Names)
 
-	total = DataFrame(["Total" 0 sum(dfStart.AnnualSum) fill(0.0, (1,T))], :auto)
-	total[:, 4:T+3] .= sum(start, dims = 1)
-	rename!(total,auxNew_Names)
-	dfStart = vcat(dfStart, total)
-	CSV.write(joinpath(path, "start.csv"), dftranspose(dfStart, false), header=false)
+	filepath = joinpath(path, "start.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfStart)
+	else 	# setup["WriteOutputs"] == "full"	
+		write_fulltimeseries(filepath, start, dfStart)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/write_capacityfactor.jl b/src/write_outputs/write_capacityfactor.jl
index 959f4366e2..03c2a50e4b 100644
--- a/src/write_outputs/write_capacityfactor.jl
+++ b/src/write_outputs/write_capacityfactor.jl
@@ -52,5 +52,5 @@ function write_capacityfactor(path::AbstractString, inputs::Dict, setup::Dict, E
     end
 
     CSV.write(joinpath(path, "capacityfactor.csv"), dfCapacityfactor)
-    return dfCapacityfactor
+    return nothing
 end
diff --git a/src/write_outputs/write_charge.jl b/src/write_outputs/write_charge.jl
index 7913799c6f..74d00ad65a 100644
--- a/src/write_outputs/write_charge.jl
+++ b/src/write_outputs/write_charge.jl
@@ -34,14 +34,12 @@ function write_charge(path::AbstractString, inputs::Dict, setup::Dict, EP::Model
 	end
 
 	dfCharge.AnnualSum .= charge * inputs["omega"]
-	dfCharge = hcat(dfCharge, DataFrame(charge, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfCharge,auxNew_Names)
-	total = DataFrame(["Total" 0 sum(dfCharge[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
 
-	total[:, 4:T+3] .= sum(charge, dims = 1)
-	rename!(total,auxNew_Names)
-	dfCharge = vcat(dfCharge, total)
-	CSV.write(joinpath(path, "charge.csv"), dftranspose(dfCharge, false), header=false)
-	return dfCharge
+	filepath = joinpath(path, "charge.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfCharge)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filepath, charge, dfCharge)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/write_co2.jl b/src/write_outputs/write_co2.jl
index 6fa2edf574..77faf3a1ee 100644
--- a/src/write_outputs/write_co2.jl
+++ b/src/write_outputs/write_co2.jl
@@ -13,25 +13,23 @@ end
 function write_co2_emissions_plant(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
     gen = inputs["RESOURCES"]
     G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
-    T = inputs["T"]     # Number of time steps (hours)
-    Z = inputs["Z"]     # Number of zones
 
     # CO2 emissions by plant
     dfEmissions_plant = DataFrame(Resource=inputs["RESOURCE_NAMES"], Zone=zone_id.(gen), AnnualSum=zeros(G))
     emissions_plant = value.(EP[:eEmissionsByPlant])
+
     if setup["ParameterScale"] == 1
         emissions_plant *= ModelScalingFactor
     end
     dfEmissions_plant.AnnualSum .= emissions_plant * inputs["omega"]
-    dfEmissions_plant = hcat(dfEmissions_plant, DataFrame(emissions_plant, :auto))
 
-    auxNew_Names = [Symbol("Resource"); Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t = 1:T]]
-    rename!(dfEmissions_plant, auxNew_Names)
-
-    total = DataFrame(["Total" 0 sum(dfEmissions_plant[!, :AnnualSum]) fill(0.0, (1, T))], auxNew_Names)
-    total[!, 4:T+3] .= sum(emissions_plant, dims=1)
-    dfEmissions_plant = vcat(dfEmissions_plant, total)
-    CSV.write(joinpath(path, "emissions_plant.csv"), dftranspose(dfEmissions_plant, false), header=false)
+    filepath = joinpath(path, "emissions_plant.csv")
+    if setup["WriteOutputs"] == "annual"
+        write_annual(filepath, dfEmissions_plant)
+    else 	# setup["WriteOutputs"] == "full"
+        write_fulltimeseries(filepath, emissions_plant, dfEmissions_plant)
+    end
+    return nothing
 end
 
 function write_co2_capture_plant(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
@@ -44,21 +42,19 @@ function write_co2_capture_plant(path::AbstractString, inputs::Dict, setup::Dict
     dfCapturedEmissions_plant = DataFrame(Resource=inputs["RESOURCE_NAMES"][CCS], Zone=zone_id.(gen[CCS]), AnnualSum=zeros(length(CCS)))
     if !isempty(CCS)
         # Captured CO2 emissions by plant
-        emissions_captured_plant = zeros(length(CCS), T)
-        emissions_captured_plant = (value.(EP[:eEmissionsCaptureByPlant]).data)
+        emissions_captured_plant = (value.(EP[:eEmissionsCaptureByPlant]))
+
         if setup["ParameterScale"] == 1
             emissions_captured_plant *= ModelScalingFactor
         end
         dfCapturedEmissions_plant.AnnualSum .= emissions_captured_plant * inputs["omega"]
-        dfCapturedEmissions_plant = hcat(dfCapturedEmissions_plant, DataFrame(emissions_captured_plant, :auto))
-
-        auxNew_Names = [Symbol("Resource"); Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t = 1:T]]
-        rename!(dfCapturedEmissions_plant, auxNew_Names)
-
-        total = DataFrame(["Total" 0 sum(dfCapturedEmissions_plant[!, :AnnualSum]) fill(0.0, (1, T))], auxNew_Names)
-        total[!, 4:T+3] .= sum(emissions_captured_plant, dims=1)
-        dfCapturedEmissions_plant = vcat(dfCapturedEmissions_plant, total)
 
-        CSV.write(joinpath(path, "captured_emissions_plant.csv"), dftranspose(dfCapturedEmissions_plant, false), header=false)
+        filepath = joinpath(path, "captured_emissions_plant.csv")
+        if setup["WriteOutputs"] == "annual"
+            write_annual(filepath, dfCapturedEmissions_plant)
+        else     # setup["WriteOutputs"] == "full"
+            write_fulltimeseries(filepath, emissions_captured_plant, dfCapturedEmissions_plant)
+        end
+        return nothing
     end
-end
+end
\ No newline at end of file
diff --git a/src/write_outputs/write_curtailment.jl b/src/write_outputs/write_curtailment.jl
index 76061269a5..6cb151f448 100644
--- a/src/write_outputs/write_curtailment.jl
+++ b/src/write_outputs/write_curtailment.jl
@@ -33,15 +33,14 @@ function write_curtailment(path::AbstractString, inputs::Dict, setup::Dict, EP::
 				+ (value.(EP[:eTotalCap_WIND][SOLAR_WIND]).data .* inputs["pP_Max_Wind"][SOLAR_WIND, :] .- value.(EP[:vP_WIND][SOLAR_WIND, :]).data))
 		end
 	end
+
 	dfCurtailment.AnnualSum = curtailment * inputs["omega"]
-	dfCurtailment = hcat(dfCurtailment, DataFrame(curtailment, :auto))
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfCurtailment,auxNew_Names)
 
-	total = DataFrame(["Total" 0 sum(dfCurtailment[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-	total[:, 4:T+3] .= sum(curtailment, dims = 1)
-	rename!(total,auxNew_Names)
-	dfCurtailment = vcat(dfCurtailment, total)
-	CSV.write(joinpath(path, "curtail.csv"), dftranspose(dfCurtailment, false), header=false)
-	return dfCurtailment
+	filename = joinpath(path, "curtail.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filename, dfCurtailment)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filename, curtailment, dfCurtailment)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/write_emissions.jl b/src/write_outputs/write_emissions.jl
index 53129e81bc..036ab35667 100644
--- a/src/write_outputs/write_emissions.jl
+++ b/src/write_outputs/write_emissions.jl
@@ -11,7 +11,6 @@ function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Mo
 
 	scale_factor = setup["ParameterScale"] == 1 ? ModelScalingFactor : 1
 
-
 	if (setup["WriteShadowPrices"]==1 || setup["UCommit"]==0 || (setup["UCommit"]==2 && (setup["Reserves"]==0 || (setup["Reserves"]>0 && inputs["pDynamic_Contingency"]==0)))) # fully linear model
 		# CO2 emissions by zone
 
@@ -34,48 +33,64 @@ function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Mo
 			dfEmissions = DataFrame(Zone = 1:Z, AnnualSum = Array{Float64}(undef, Z))
 		end
 
+		emissions_by_zone = value.(EP[:eEmissionsByZone])
 		for i in 1:Z
-			dfEmissions[i,:AnnualSum] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone][i,:]))*scale_factor
+			dfEmissions[i,:AnnualSum] = sum(inputs["omega"] .* emissions_by_zone[i,:]) * scale_factor
 		end
 
-		dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone])*scale_factor, :auto))
-
-
-		if setup["CO2Cap"]>=1
-			auxNew_Names=[Symbol("Zone");[Symbol("CO2_Price_$cap") for cap in 1:inputs["NCO2Cap"]];Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-			rename!(dfEmissions,auxNew_Names)
-			total = DataFrame(["Total" zeros(1,inputs["NCO2Cap"]) sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-			for t in 1:T
-				total[:,t+inputs["NCO2Cap"]+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+		if setup["WriteOutputs"] == "annual"
+			total = DataFrame(["Total" sum(dfEmissions.AnnualSum)], [:Zone;:AnnualSum])
+			if setup["CO2Cap"]>=1
+				total = DataFrame(["Total" zeros(1,inputs["NCO2Cap"]) sum(dfEmissions.AnnualSum)], [:Zone;[Symbol("CO2_Price_$cap") for cap in 1:inputs["NCO2Cap"]];:AnnualSum])
 			end
-		else
-			auxNew_Names=[Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t in 1:T]]
-			rename!(dfEmissions,auxNew_Names)
-			total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-			for t in 1:T
-				total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+			dfEmissions = vcat(dfEmissions, total)
+			CSV.write(joinpath(path, "emissions.csv"), dfEmissions)
+		else	# setup["WriteOutputs"] == "full"
+			dfEmissions = hcat(dfEmissions, DataFrame(emissions_by_zone * scale_factor, :auto))
+			if setup["CO2Cap"]>=1
+				auxNew_Names=[Symbol("Zone");[Symbol("CO2_Price_$cap") for cap in 1:inputs["NCO2Cap"]];Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+				rename!(dfEmissions,auxNew_Names)
+				total = DataFrame(["Total" zeros(1,inputs["NCO2Cap"]) sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
+				for t in 1:T
+					total[:,t+inputs["NCO2Cap"]+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+				end
+			else
+				auxNew_Names=[Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t in 1:T]]
+				rename!(dfEmissions,auxNew_Names)
+				total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
+				for t in 1:T
+					total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+				end
 			end
+			rename!(total,auxNew_Names)
+			dfEmissions = vcat(dfEmissions, total)
+			CSV.write(joinpath(path, "emissions.csv"), dftranspose(dfEmissions, false), writeheader=false)
 		end
-        rename!(total,auxNew_Names)
-        dfEmissions = vcat(dfEmissions, total)
-
-
 ## Aaron - Combined elseif setup["Dual_MIP"]==1 block with the first block since they were identical. Why do we have this third case? What is different about it?
 	else
 		# CO2 emissions by zone
+		emissions_by_zone = value.(EP[:eEmissionsByZone])
 		dfEmissions = hcat(DataFrame(Zone = 1:Z), DataFrame(AnnualSum = Array{Float64}(undef, Z)))
 		for i in 1:Z
-			dfEmissions[i,:AnnualSum] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone][i,:])) * scale_factor
+			dfEmissions[i,:AnnualSum] = sum(inputs["omega"] .* emissions_by_zone[i,:]) * scale_factor
 		end
-		dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone]) * scale_factor, :auto))
-		auxNew_Names=[Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfEmissions,auxNew_Names)
-		total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		for t in 1:T
-			total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+
+		if setup["WriteOutputs"] == "annual"
+			total = DataFrame(["Total" sum(dfEmissions.AnnualSum)], [:Zone;:AnnualSum])
+			dfEmissions = vcat(dfEmissions, total)
+			CSV.write(joinpath(path, "emissions.csv"), dfEmissions)
+		else	# setup["WriteOutputs"] == "full"
+			dfEmissions = hcat(dfEmissions, DataFrame(emissions_by_zone * scale_factor, :auto))
+			auxNew_Names=[Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+			rename!(dfEmissions,auxNew_Names)
+			total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
+			for t in 1:T
+				total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
+			end
+			rename!(total,auxNew_Names)
+			dfEmissions = vcat(dfEmissions, total)
+			CSV.write(joinpath(path, "emissions.csv"), dftranspose(dfEmissions, false), writeheader=false)
 		end
-		rename!(total,auxNew_Names)
-		dfEmissions = vcat(dfEmissions, total)
 	end
-	CSV.write(joinpath(path, "emissions.csv"), dftranspose(dfEmissions, false), header=false)
+	return nothing
 end
diff --git a/src/write_outputs/write_fuel_consumption.jl b/src/write_outputs/write_fuel_consumption.jl
index 427e318a4c..7a661b9386 100644
--- a/src/write_outputs/write_fuel_consumption.jl
+++ b/src/write_outputs/write_fuel_consumption.jl
@@ -6,7 +6,9 @@ Write fuel consumption of each power plant.
 function write_fuel_consumption(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
 
 	write_fuel_consumption_plant(path::AbstractString,inputs::Dict, setup::Dict, EP::Model)
-	write_fuel_consumption_ts(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+	if setup["WriteOutputs"] != "annual"
+		write_fuel_consumption_ts(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
+	end
 	write_fuel_consumption_tot(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
 end
 
diff --git a/src/write_outputs/write_nse.jl b/src/write_outputs/write_nse.jl
index ae896fe416..5d30dcc987 100644
--- a/src/write_outputs/write_nse.jl
+++ b/src/write_outputs/write_nse.jl
@@ -15,15 +15,22 @@ function write_nse(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
 		nse[((z-1)*SEG+1):z*SEG, :] = value.(EP[:vNSE])[:, :, z] * scale_factor
 	end
 	dfNse.AnnualSum .= nse * inputs["omega"]
-	dfNse = hcat(dfNse, DataFrame(nse, :auto))
-	auxNew_Names=[Symbol("Segment");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfNse,auxNew_Names)
 
-	total = DataFrame(["Total" 0 sum(dfNse[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-	total[:, 4:T+3] .= sum(nse, dims = 1)
-	rename!(total,auxNew_Names)
-	dfNse = vcat(dfNse, total)
+	if setup["WriteOutputs"] == "annual"
+		total = DataFrame(["Total" 0 sum(dfNse[!,:AnnualSum])], [:Segment, :Zone, :AnnualSum])
+		dfNse = vcat(dfNse, total)
+		CSV.write(joinpath(path, "nse.csv"), dfNse)
+	else 	# setup["WriteOutputs"] == "full"
+		dfNse = hcat(dfNse, DataFrame(nse, :auto))
+		auxNew_Names=[Symbol("Segment");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+		rename!(dfNse,auxNew_Names)
 
-	CSV.write(joinpath(path, "nse.csv"),  dftranspose(dfNse, false), header=false)
-	return dfNse
+		total = DataFrame(["Total" 0 sum(dfNse[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
+		total[:, 4:T+3] .= sum(nse, dims = 1)
+		rename!(total,auxNew_Names)
+		dfNse = vcat(dfNse, total)
+
+		CSV.write(joinpath(path, "nse.csv"),  dftranspose(dfNse, false), writeheader=false)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/write_outputs.jl b/src/write_outputs/write_outputs.jl
index 0cef6db417..113100120e 100644
--- a/src/write_outputs/write_outputs.jl
+++ b/src/write_outputs/write_outputs.jl
@@ -41,44 +41,99 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
 		end
 	end
 
-	write_status(path, inputs, setup, EP)
-	elapsed_time_costs = @elapsed write_costs(path, inputs, setup, EP)
-	println("Time elapsed for writing costs is")
-	println(elapsed_time_costs)
-	dfCap = write_capacity(path, inputs, setup, EP)
-	dfPower = write_power(path, inputs, setup, EP)
-	dfCharge = write_charge(path, inputs, setup, EP)
-	dfCapacityfactor = write_capacityfactor(path, inputs, setup, EP)
-	elapsed_time_storage = @elapsed write_storage(path, inputs, setup, EP)
-	println("Time elapsed for writing storage is")
-	println(elapsed_time_storage)
-	dfCurtailment = write_curtailment(path, inputs, setup, EP)
-	elapsed_time_nse = @elapsed write_nse(path, inputs, setup, EP)
-	println("Time elapsed for writing nse is")
-	println(elapsed_time_nse)
-	elapsed_time_power_balance = @elapsed write_power_balance(path, inputs, setup, EP)
-	println("Time elapsed for writing power balance is")
-	println(elapsed_time_power_balance)
+	# Dict containing the list of outputs to write
+	output_settings_d = setup["WriteOutputsSettingsDict"]
+	write_settings_file(path, setup)
+	
+	output_settings_d["WriteStatus"] && write_status(path, inputs, setup, EP)
+
+	if output_settings_d["WriteCosts"]
+		elapsed_time_costs = @elapsed write_costs(path, inputs, setup, EP)
+		println("Time elapsed for writing costs is")
+		println(elapsed_time_costs)
+	end
+	
+	if output_settings_d["WriteCapacity"] || output_settings_d["WriteNetRevenue"]
+		elapsed_time_capacity = @elapsed dfCap = write_capacity(path, inputs, setup, EP)
+		println("Time elapsed for writing capacity is")
+		println(elapsed_time_capacity)
+	end
+
+	if output_settings_d["WritePower"] || output_settings_d["WriteNetRevenue"]
+		elapsed_time_power = @elapsed dfPower = write_power(path, inputs, setup, EP)
+		println("Time elapsed for writing power is")
+		println(elapsed_time_power)
+	end
+
+	if output_settings_d["WriteCharge"]
+		elapsed_time_charge = @elapsed write_charge(path, inputs, setup, EP)
+		println("Time elapsed for writing charge is")
+		println(elapsed_time_charge)
+	end
+
+	if output_settings_d["WriteCapacityFactor"]
+		elapsed_time_capacityfactor = @elapsed write_capacityfactor(path, inputs, setup, EP)
+		println("Time elapsed for writing capacity factor is")
+		println(elapsed_time_capacityfactor)
+	end
+
+	if output_settings_d["WriteStorage"]
+		elapsed_time_storage = @elapsed write_storage(path, inputs, setup, EP)
+		println("Time elapsed for writing storage is")
+		println(elapsed_time_storage)
+	end
+
+	if output_settings_d["WriteCurtailment"]
+		elapsed_time_curtailment = @elapsed write_curtailment(path, inputs, setup, EP)
+		println("Time elapsed for writing curtailment is")
+		println(elapsed_time_curtailment)
+	end
+
+	if output_settings_d["WriteNSE"]
+		elapsed_time_nse = @elapsed write_nse(path, inputs, setup, EP)
+		println("Time elapsed for writing nse is")
+		println(elapsed_time_nse)
+	end
+
+	if output_settings_d["WritePowerBalance"]
+		elapsed_time_power_balance = @elapsed write_power_balance(path, inputs, setup, EP)
+		println("Time elapsed for writing power balance is")
+		println(elapsed_time_power_balance)
+	end
+
 	if inputs["Z"] > 1
-		elapsed_time_flows = @elapsed write_transmission_flows(path, inputs, setup, EP)
-		println("Time elapsed for writing transmission flows is")
-		println(elapsed_time_flows)
-		elapsed_time_losses = @elapsed write_transmission_losses(path, inputs, setup, EP)
-		println("Time elapsed for writing transmission losses is")
-		println(elapsed_time_losses)
-		if setup["NetworkExpansion"] == 1
+		if output_settings_d["WriteTransmissionFlows"]
+			elapsed_time_flows = @elapsed write_transmission_flows(path, inputs, setup, EP)
+			println("Time elapsed for writing transmission flows is")
+			println(elapsed_time_flows)
+		end
+
+		if output_settings_d["WriteTransmissionLosses"]
+			elapsed_time_losses = @elapsed write_transmission_losses(path, inputs, setup, EP)
+			println("Time elapsed for writing transmission losses is")
+			println(elapsed_time_losses)
+		end
+
+		if setup["NetworkExpansion"] == 1 && output_settings_d["WriteNWExpansion"]
 			elapsed_time_expansion = @elapsed write_nw_expansion(path, inputs, setup, EP)
 			println("Time elapsed for writing network expansion is")
 			println(elapsed_time_expansion)
 		end
 	end
-	elapsed_time_emissions = @elapsed write_emissions(path, inputs, setup, EP)
-	println("Time elapsed for writing emissions is")
-	println(elapsed_time_emissions)
+
+	if output_settings_d["WriteEmissions"]
+		elapsed_time_emissions = @elapsed write_emissions(path, inputs, setup, EP)
+		println("Time elapsed for writing emissions is")
+		println(elapsed_time_emissions)
+	end
 
 	dfVreStor = DataFrame()
 	if !isempty(inputs["VRE_STOR"])
-		dfVreStor = write_vre_stor(path, inputs, setup, EP)
+		if output_settings_d["WriteVREStor"] || output_settings_d["WriteNetRevenue"]
+			elapsed_time_vrestor = @elapsed dfVreStor = write_vre_stor(path, inputs, setup, EP)
+			println("Time elapsed for writing vre stor is")
+			println(elapsed_time_vrestor)
+		end
 		VS_LDS = inputs["VS_LDS"]
 		VS_STOR = inputs["VS_STOR"]
 	else
@@ -87,97 +142,178 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
 	end
 
 	if has_duals(EP) == 1
-		elapsed_time_reliability = @elapsed write_reliability(path, inputs, setup, EP)
-		println("Time elapsed for writing reliability is")
-		println(elapsed_time_reliability)
+		if output_settings_d["WriteReliability"]
+			elapsed_time_reliability = @elapsed write_reliability(path, inputs, setup, EP)
+			println("Time elapsed for writing reliability is")
+			println(elapsed_time_reliability)
+		end
 		if !isempty(inputs["STOR_ALL"]) || !isempty(VS_STOR)
-			elapsed_time_stordual = @elapsed write_storagedual(path, inputs, setup, EP)
-			println("Time elapsed for writing storage duals is")
-			println(elapsed_time_stordual)
+			if output_settings_d["WriteStorageDual"]
+				elapsed_time_stordual = @elapsed write_storagedual(path, inputs, setup, EP)
+				println("Time elapsed for writing storage duals is")
+				println(elapsed_time_stordual)
+			end
 		end
 	end
 
 	if setup["UCommit"] >= 1
-		elapsed_time_commit = @elapsed write_commit(path, inputs, setup, EP)
-		println("Time elapsed for writing commitment is")
-		println(elapsed_time_commit)
-		elapsed_time_start = @elapsed write_start(path, inputs, setup, EP)
-		println("Time elapsed for writing startup is")
-		println(elapsed_time_start)
-		elapsed_time_shutdown = @elapsed write_shutdown(path, inputs, setup, EP)
-		println("Time elapsed for writing shutdown is")
-		println(elapsed_time_shutdown)
+		if output_settings_d["WriteCommit"]
+			elapsed_time_commit = @elapsed write_commit(path, inputs, setup, EP)
+			println("Time elapsed for writing commitment is")
+			println(elapsed_time_commit)
+		end
+
+		if output_settings_d["WriteStart"]
+			elapsed_time_start = @elapsed write_start(path, inputs, setup, EP)
+			println("Time elapsed for writing startup is")
+			println(elapsed_time_start)
+		end
+
+		if output_settings_d["WriteShutdown"]
+			elapsed_time_shutdown = @elapsed write_shutdown(path, inputs, setup, EP)
+			println("Time elapsed for writing shutdown is")
+			println(elapsed_time_shutdown)
+		end
 		if setup["Reserves"] == 1
-			elapsed_time_reg = @elapsed write_reg(path, inputs, setup, EP)
-			println("Time elapsed for writing regulation is")
-			println(elapsed_time_reg)
-			elapsed_time_rsv = @elapsed write_rsv(path, inputs, setup, EP)
-			println("Time elapsed for writing reserves is")
-			println(elapsed_time_rsv)
+			if output_settings_d["WriteReg"]
+				elapsed_time_reg = @elapsed write_reg(path, inputs, setup, EP)
+				println("Time elapsed for writing regulation is")
+				println(elapsed_time_reg)
+			end
+
+			if output_settings_d["WriteRsv"]
+				elapsed_time_rsv = @elapsed write_rsv(path, inputs, setup, EP)
+				println("Time elapsed for writing reserves is")
+				println(elapsed_time_rsv)
+			end
 		end
 	end
 
 	# Output additional variables related inter-period energy transfer via storage
 	representative_periods = inputs["REP_PERIOD"]
 	if representative_periods > 1 && (!isempty(inputs["STOR_LONG_DURATION"]) || !isempty(VS_LDS))
-		elapsed_time_lds_init = @elapsed write_opwrap_lds_stor_init(path, inputs, setup, EP)
-		println("Time elapsed for writing lds init is")
-		println(elapsed_time_lds_init)
-		elapsed_time_lds_dstor = @elapsed write_opwrap_lds_dstor(path, inputs, setup, EP)
-		println("Time elapsed for writing lds dstor is")
-		println(elapsed_time_lds_dstor)
+		if output_settings_d["WriteOpWrapLDSStorInit"]
+			elapsed_time_lds_init = @elapsed write_opwrap_lds_stor_init(path, inputs, setup, EP)
+			println("Time elapsed for writing lds init is")
+			println(elapsed_time_lds_init)
+		end
+
+		if output_settings_d["WriteOpWrapLDSdStor"]
+			elapsed_time_lds_dstor = @elapsed write_opwrap_lds_dstor(path, inputs, setup, EP)
+			println("Time elapsed for writing lds dstor is")
+			println(elapsed_time_lds_dstor)
+		end
 	end
 
-	elapsed_time_fuel_consumption = @elapsed write_fuel_consumption(path, inputs, setup, EP)
-	println("Time elapsed for writing fuel consumption is")
-	println(elapsed_time_fuel_consumption)
+	if output_settings_d["WriteFuelConsumption"]
+		elapsed_time_fuel_consumption = @elapsed write_fuel_consumption(path, inputs, setup, EP)
+		println("Time elapsed for writing fuel consumption is")
+		println(elapsed_time_fuel_consumption)
+	end
 
-	elapsed_time_emissions = @elapsed write_co2(path, inputs, setup, EP)
-	println("Time elapsed for writing co2 is")
-	println(elapsed_time_emissions)
+	if output_settings_d["WriteCO2"]
+		elapsed_time_emissions = @elapsed write_co2(path, inputs, setup, EP)
+		println("Time elapsed for writing co2 is")
+		println(elapsed_time_emissions)
+	end
 
-    if has_maintenance(inputs)
+    if has_maintenance(inputs) && output_settings_d["WriteMaintenance"]
         write_maintenance(path, inputs, EP)
     end
 
 	# Temporary! Suppress these outputs until we know that they are compatable with multi-stage modeling
 	if setup["MultiStage"] == 0
-		dfPrice = DataFrame()
 		dfEnergyRevenue = DataFrame()
 		dfChargingcost = DataFrame()
 		dfSubRevenue = DataFrame()
 		dfRegSubRevenue = DataFrame()
 		if has_duals(EP) == 1
-			dfPrice = write_price(path, inputs, setup, EP)
-			dfEnergyRevenue = write_energy_revenue(path, inputs, setup, EP)
-			dfChargingcost = write_charging_cost(path, inputs, setup, EP)
-			dfSubRevenue, dfRegSubRevenue = write_subsidy_revenue(path, inputs, setup, EP)
+			if output_settings_d["WritePrice"]
+				elapsed_time_price = @elapsed write_price(path, inputs, setup, EP)
+				println("Time elapsed for writing price is")
+				println(elapsed_time_price)
+			end
+
+			if output_settings_d["WriteEnergyRevenue"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_energy_rev = @elapsed dfEnergyRevenue = write_energy_revenue(path, inputs, setup, EP)
+				println("Time elapsed for writing energy revenue is")
+				println(elapsed_time_energy_rev)
+			end
+
+			if output_settings_d["WriteChargingCost"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_charging_cost = @elapsed dfChargingcost = write_charging_cost(path, inputs, setup, EP)
+				println("Time elapsed for writing charging cost is")
+				println(elapsed_time_charging_cost)
+			end
+
+			if output_settings_d["WriteSubsidyRevenue"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_subsidy = @elapsed dfSubRevenue, dfRegSubRevenue = write_subsidy_revenue(path, inputs, setup, EP)
+				println("Time elapsed for writing subsidy is")
+				println(elapsed_time_subsidy)
+			end
+		end
+
+		if output_settings_d["WriteTimeWeights"] 
+			elapsed_time_time_weights = @elapsed write_time_weights(path, inputs)
+			println("Time elapsed for writing time weights is")
+			println(elapsed_time_time_weights)
 		end
 
-		elapsed_time_time_weights = @elapsed write_time_weights(path, inputs)
-	    println("Time elapsed for writing time weights is")
-	    println(elapsed_time_time_weights)
-		dfESR = DataFrame()
 		dfESRRev = DataFrame()
-		if setup["EnergyShareRequirement"]==1 && has_duals(EP)
-			dfESR = write_esr_prices(path, inputs, setup, EP)
-			dfESRRev = write_esr_revenue(path, inputs, setup, dfPower, dfESR, EP)
+		if setup["EnergyShareRequirement"] == 1 && has_duals(EP)
+			dfESR = DataFrame()
+			if output_settings_d["WriteESRPrices"] || output_settings_d["WriteESRRevenue"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_esr_prices = @elapsed dfESR = write_esr_prices(path, inputs, setup, EP)
+				println("Time elapsed for writing esr prices is")
+				println(elapsed_time_esr_prices)
+			end 
+
+			if output_settings_d["WriteESRRevenue"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_esr_revenue = @elapsed dfESRRev = write_esr_revenue(path, inputs, setup, dfPower, dfESR, EP)
+				println("Time elapsed for writing esr revenue is")
+				println(elapsed_time_esr_revenue)
+			end
+
 		end
-		dfResMar = DataFrame()
+
 		dfResRevenue = DataFrame()
 		if setup["CapacityReserveMargin"]==1 && has_duals(EP)
-			dfResMar = write_reserve_margin(path, setup, EP)
-			elapsed_time_rsv_margin = @elapsed write_reserve_margin_w(path, inputs, setup, EP)
-			dfVirtualDischarge = write_virtual_discharge(path, inputs, setup, EP)
-		  println("Time elapsed for writing reserve margin is")
-		  println(elapsed_time_rsv_margin)
-			dfResRevenue = write_reserve_margin_revenue(path, inputs, setup, EP)
-			elapsed_time_cap_value = @elapsed write_capacity_value(path, inputs, setup, EP)
-		  println("Time elapsed for writing capacity value is")
-		  println(elapsed_time_cap_value)
-			if haskey(inputs, "dfCapRes_slack")
-				dfResMar_slack = write_reserve_margin_slack(path, inputs, setup, EP)
-			end		  
+			if output_settings_d["WriteReserveMargin"]
+				elapsed_time_reserve_margin = @elapsed write_reserve_margin(path, setup, EP)
+				println("Time elapsed for writing reserve margin is")
+				println(elapsed_time_reserve_margin)
+			end
+
+			if output_settings_d["WriteReserveMarginWithWeights"]
+				elapsed_time_rsv_margin_w = @elapsed write_reserve_margin_w(path, inputs, setup, EP)
+				println("Time elapsed for writing reserve margin with weights is")
+				println(elapsed_time_rsv_margin_w)
+			end
+
+			if output_settings_d["WriteVirtualDischarge"]
+				elapsed_time_virtual_discharge = @elapsed write_virtual_discharge(path, inputs, setup, EP)
+				println("Time elapsed for writing virtual discharge is")
+				println(elapsed_time_virtual_discharge)
+			end
+			
+			if output_settings_d["WriteReserveMarginRevenue"] || output_settings_d["WriteNetRevenue"]
+				elapsed_time_res_rev = @elapsed dfResRevenue = write_reserve_margin_revenue(path, inputs, setup, EP)
+				println("Time elapsed for writing reserve revenue is")
+				println(elapsed_time_res_rev)
+			end
+
+			if haskey(inputs, "dfCapRes_slack") && output_settings_d["WriteReserveMarginSlack"]
+				elapsed_time_rsv_slack = @elapsed write_reserve_margin_slack(path, inputs, setup, EP)
+				println("Time elapsed for writing reserve margin slack is")
+				println(elapsed_time_rsv_slack)
+			end	
+
+			if output_settings_d["WriteCapacityValue"]
+				elapsed_time_cap_value = @elapsed write_capacity_value(path, inputs, setup, EP)
+				println("Time elapsed for writing capacity value is")
+				println(elapsed_time_cap_value)
+			end
+			
 		end
 
 		dfOpRegRevenue = DataFrame()
@@ -187,33 +323,76 @@ function write_outputs(EP::Model, path::AbstractString, setup::Dict, inputs::Dic
 		  	println("Time elapsed for writing oerating reserve and regulation revenue is")
 		  	println(elapsed_time_op_res_rev)
 		end
-
-		if setup["CO2Cap"]>0 && has_duals(EP)
-			dfCO2Cap = write_co2_cap(path, inputs, setup, EP)
+		
+		if setup["CO2Cap"]>0 && has_duals(EP) == 1 && output_settings_d["WriteCO2Cap"]
+			elapsed_time_co2_cap = @elapsed	write_co2_cap(path, inputs, setup, EP)
+			println("Time elapsed for writing co2 cap is")
+			println(elapsed_time_co2_cap)
 		end
-		if setup["MinCapReq"] == 1 && has_duals(EP)
-			dfMinCapReq = write_minimum_capacity_requirement(path, inputs, setup, EP)
+		if setup["MinCapReq"] == 1 && has_duals(EP) == 1 && output_settings_d["WriteMinCapReq"]
+			elapsed_time_min_cap_req = @elapsed write_minimum_capacity_requirement(path, inputs, setup, EP)
+			println("Time elapsed for writing minimum capacity requirement is")
+			println(elapsed_time_min_cap_req)
 		end
 
-		if setup["MaxCapReq"] == 1 && has_duals(EP)
-			dfMaxCapReq = write_maximum_capacity_requirement(path, inputs, setup, EP)
+		if setup["MaxCapReq"] == 1 && has_duals(EP) == 1 && output_settings_d["WriteMaxCapReq"]
+			elapsed_time_max_cap_req = @elapsed write_maximum_capacity_requirement(path, inputs, setup, EP)
+			println("Time elapsed for writing maximum capacity requirement is")
+			println(elapsed_time_max_cap_req)
 		end
 
 		if !isempty(inputs["ELECTROLYZER"]) && has_duals(EP)
-			dfHydrogenPrice = write_hydrogen_prices(path, inputs, setup, EP)
-			if setup["HydrogenHourlyMatching"] == 1
-				dfHourlyMatchingPrices = write_hourly_matching_prices(path, inputs, setup, EP)
+			if output_settings_d["WriteHydrogenPrices"]
+				elapsed_time_hydrogen_prices = @elapsed	write_hydrogen_prices(path, inputs, setup, EP)
+				println("Time elapsed for writing hydrogen prices is")
+				println(elapsed_time_hydrogen_prices)
+			end
+			if setup["HydrogenHourlyMatching"] == 1 && output_settings_d["WriteHourlyMatchingPrices"]
+				elapsed_time_hourly_matching_prices = @elapsed	write_hourly_matching_prices(path, inputs, setup, EP)
+				println("Time elapsed for writing hourly matching prices is")
+				println(elapsed_time_hourly_matching_prices)
 			end
 		end
 
-
-		elapsed_time_net_rev = @elapsed write_net_revenue(path, inputs, setup, EP, dfCap, dfESRRev, dfResRevenue, dfChargingcost, dfPower, dfEnergyRevenue, dfSubRevenue, dfRegSubRevenue, dfVreStor, dfOpRegRevenue, dfOpRsvRevenue)
-	  	println("Time elapsed for writing net revenue is")
-	  	println(elapsed_time_net_rev)
+		if output_settings_d["WriteNetRevenue"]
+			elapsed_time_net_rev = @elapsed write_net_revenue(path, inputs, setup, EP, dfCap, dfESRRev, dfResRevenue, dfChargingcost, dfPower, dfEnergyRevenue, dfSubRevenue, dfRegSubRevenue, dfVreStor, dfOpRegRevenue, dfOpRsvRevenue)
+		  	println("Time elapsed for writing net revenue is")
+	  		println(elapsed_time_net_rev)
+		end
 	end
 	## Print confirmation
 	println("Wrote outputs to $path")
 
 	return path
-
 end # END output()
+
+
+"""
+	write_annual(fullpath::AbstractString, dfOut::DataFrame)
+
+Internal function for writing annual outputs. 
+"""
+function write_annual(fullpath::AbstractString, dfOut::DataFrame)
+	push!(dfOut, ["Total" 0 sum(dfOut[!, :AnnualSum])])
+	CSV.write(fullpath, dfOut)
+	return nothing
+end
+
+"""
+	write_fulltimeseries(fullpath::AbstractString, dataOut::Matrix{Float64}, dfOut::DataFrame)
+
+Internal function for writing full time series outputs. This function wraps the instructions for creating the full time series output files. 
+"""
+function write_fulltimeseries(fullpath::AbstractString, dataOut::Matrix{Float64}, dfOut::DataFrame)
+	T = size(dataOut, 2)
+	dfOut = hcat(dfOut, DataFrame(dataOut, :auto))
+	auxNew_Names = [Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
+	rename!(dfOut, auxNew_Names)
+	total = DataFrame(["Total" 0 sum(dfOut[!, :AnnualSum]) fill(0.0, (1, T))], auxNew_Names)
+	total[!, 4:T+3] .= sum(dataOut, dims=1)
+	dfOut = vcat(dfOut, total)
+	CSV.write(fullpath, dftranspose(dfOut, false), writeheader=false)
+	return nothing
+end
+
+write_settings_file(path, setup) = YAML.write_file(joinpath(path, "run_settings.yml"), setup)
diff --git a/src/write_outputs/write_power.jl b/src/write_outputs/write_power.jl
index 7e1cc986da..3be5e83bf3 100644
--- a/src/write_outputs/write_power.jl
+++ b/src/write_outputs/write_power.jl
@@ -17,16 +17,13 @@ function write_power(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
 		power *= ModelScalingFactor
 	end
 	dfPower.AnnualSum .= power * inputs["omega"]
-	dfPower = hcat(dfPower, DataFrame(power, :auto))
 
-	auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-	rename!(dfPower,auxNew_Names)
-
-	total = DataFrame(["Total" 0 sum(dfPower[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-	total[:, 4:T+3] .= sum(power, dims = 1)
+	filepath = joinpath(path, "power.csv")
+	if setup["WriteOutputs"] == "annual"
+		write_annual(filepath, dfPower)
+	else 	# setup["WriteOutputs"] == "full"
+		write_fulltimeseries(filepath, power, dfPower)
+	end
 
-	rename!(total,auxNew_Names)
-	dfPower = vcat(dfPower, total)
-	CSV.write(joinpath(path, "power.csv"), dftranspose(dfPower, false), header=false)
-	return dfPower
+	return dfPower #Shouldn't this be return nothing
 end
diff --git a/src/write_outputs/write_power_balance.jl b/src/write_outputs/write_power_balance.jl
index 82bd6e2b27..627f3a4821 100644
--- a/src/write_outputs/write_power_balance.jl
+++ b/src/write_outputs/write_power_balance.jl
@@ -58,8 +58,14 @@ function write_power_balance(path::AbstractString, inputs::Dict, setup::Dict, EP
 		powerbalance *= ModelScalingFactor
 	end
 	dfPowerBalance.AnnualSum .= powerbalance * inputs["omega"]
-	dfPowerBalance = hcat(dfPowerBalance, DataFrame(powerbalance, :auto))
-	auxNew_Names = [Symbol("BalanceComponent"); Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t in 1:T]]
-	rename!(dfPowerBalance,auxNew_Names)
-	CSV.write(joinpath(path, "power_balance.csv"), dftranspose(dfPowerBalance, false), header=false)
+
+	if setup["WriteOutputs"] == "annual"
+		CSV.write(joinpath(path, "power_balance.csv"), dfPowerBalance)
+	else 	# setup["WriteOutputs"] == "full"	
+		dfPowerBalance = hcat(dfPowerBalance, DataFrame(powerbalance, :auto))
+		auxNew_Names = [Symbol("BalanceComponent"); Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t in 1:T]]
+		rename!(dfPowerBalance,auxNew_Names)
+		CSV.write(joinpath(path, "power_balance.csv"), dftranspose(dfPowerBalance, false), writeheader=false)
+	end
+	return nothing
 end
diff --git a/src/write_outputs/write_price.jl b/src/write_outputs/write_price.jl
index 3bb03c130b..7d240a0f02 100644
--- a/src/write_outputs/write_price.jl
+++ b/src/write_outputs/write_price.jl
@@ -18,8 +18,8 @@ function write_price(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
 	rename!(dfPrice,auxNew_Names)
 
 	## Linear configuration final output
-	CSV.write(joinpath(path, "prices.csv"), dftranspose(dfPrice, false), header=false)
-	return dfPrice
+	CSV.write(joinpath(path, "prices.csv"), dftranspose(dfPrice, false), writeheader=false)
+	return nothing
 end
 
 @doc raw"""
diff --git a/src/write_outputs/write_vre_stor.jl b/src/write_outputs/write_vre_stor.jl
index 818cb60c58..a7af4bf5ea 100644
--- a/src/write_outputs/write_vre_stor.jl
+++ b/src/write_outputs/write_vre_stor.jl
@@ -300,14 +300,14 @@ function write_vre_stor_charge(path::AbstractString, inputs::Dict, setup::Dict,
 		charge_dc = zeros(size(DC_CHARGE)[1], T)
 		charge_dc = value.(EP[:vP_DC_CHARGE]).data ./ etainverter.(gen_VRE_STOR[(gen_VRE_STOR.stor_dc_discharge.!=0)]) * (setup["ParameterScale"]==1 ? ModelScalingFactor : 1)
 		dfCharge_DC.AnnualSum .= charge_dc * inputs["omega"]
-		dfCharge_DC = hcat(dfCharge_DC, DataFrame(charge_dc, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfCharge_DC,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfCharge_DC[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(charge_dc, dims = 1)
-		rename!(total,auxNew_Names)
-		dfCharge_DC = vcat(dfCharge_DC, total)
-		CSV.write(joinpath(path,"vre_stor_dc_charge.csv"), dftranspose(dfCharge_DC, false), header=false)
+
+
+		filepath = joinpath(path,"vre_stor_dc_charge.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfCharge_DC)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, charge_dc, dfCharge_DC)
+		end
 	end
 
 	# AC charging of battery dataframe
@@ -316,15 +316,15 @@ function write_vre_stor_charge(path::AbstractString, inputs::Dict, setup::Dict,
 		charge_ac = zeros(size(AC_CHARGE)[1], T)
 		charge_ac = value.(EP[:vP_AC_CHARGE]).data * (setup["ParameterScale"]==1 ? ModelScalingFactor : 1)
 		dfCharge_AC.AnnualSum .= charge_ac * inputs["omega"]
-		dfCharge_AC = hcat(dfCharge_AC, DataFrame(charge_ac, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfCharge_AC,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfCharge_AC[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(charge_ac, dims = 1)
-		rename!(total,auxNew_Names)
-		dfCharge_AC = vcat(dfCharge_AC, total)
-		CSV.write(joinpath(path,"vre_stor_ac_charge.csv"), dftranspose(dfCharge_AC, false), header=false)
+
+		filepath = joinpath(path,"vre_stor_ac_charge.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfCharge_AC)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, charge_ac, dfCharge_AC)
+		end
 	end
+	return nothing
 end
 
 @doc raw"""
@@ -349,14 +349,13 @@ function write_vre_stor_discharge(path::AbstractString, inputs::Dict, setup::Dic
 			power_vre_stor *= ModelScalingFactor
 		end
 		dfDischarge_DC.AnnualSum .= power_vre_stor * inputs["omega"]
-		dfDischarge_DC = hcat(dfDischarge_DC, DataFrame(power_vre_stor, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfDischarge_DC,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfDischarge_DC[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(power_vre_stor, dims = 1)
-		rename!(total,auxNew_Names)
-		dfDischarge_DC = vcat(dfDischarge_DC, total)
-		CSV.write(joinpath(path, "vre_stor_dc_discharge.csv"), dftranspose(dfDischarge_DC, false), header=false)
+
+		filepath = joinpath(path,"vre_stor_dc_discharge.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfDischarge_DC)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, power_vre_stor, dfDischarge_DC)
+		end
 	end
 
 	# AC discharging of battery dataframe
@@ -367,14 +366,13 @@ function write_vre_stor_discharge(path::AbstractString, inputs::Dict, setup::Dic
 			power_vre_stor *= ModelScalingFactor
 		end
 		dfDischarge_AC.AnnualSum .= power_vre_stor * inputs["omega"]
-		dfDischarge_AC = hcat(dfDischarge_AC, DataFrame(power_vre_stor, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfDischarge_AC,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfDischarge_AC[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(power_vre_stor, dims = 1)
-		rename!(total,auxNew_Names)
-		dfDischarge_AC = vcat(dfDischarge_AC, total)
-		CSV.write(joinpath(path, "vre_stor_ac_discharge.csv"), dftranspose(dfDischarge_AC, false), header=false)
+
+		filepath = joinpath(path,"vre_stor_ac_discharge.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfDischarge_AC)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, power_vre_stor, dfDischarge_AC)
+		end
 	end
 
 	# Wind generation of co-located resource dataframe
@@ -385,14 +383,13 @@ function write_vre_stor_discharge(path::AbstractString, inputs::Dict, setup::Dic
 			vre_vre_stor *= ModelScalingFactor
 		end
 		dfVP_VRE_STOR.AnnualSum .= vre_vre_stor * inputs["omega"]
-		dfVP_VRE_STOR = hcat(dfVP_VRE_STOR, DataFrame(vre_vre_stor, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfVP_VRE_STOR,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfVP_VRE_STOR[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(vre_vre_stor, dims = 1)
-		rename!(total,auxNew_Names)
-		dfVP_VRE_STOR = vcat(dfVP_VRE_STOR, total)
-		CSV.write(joinpath(path,"vre_stor_wind_power.csv"), dftranspose(dfVP_VRE_STOR, false), header=false)
+
+		filepath = joinpath(path,"vre_stor_wind_power.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfVP_VRE_STOR)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, vre_vre_stor, dfVP_VRE_STOR)
+		end
 	end
 
 	# Solar generation of co-located resource dataframe
@@ -403,13 +400,13 @@ function write_vre_stor_discharge(path::AbstractString, inputs::Dict, setup::Dic
 			vre_vre_stor *= ModelScalingFactor
 		end
 		dfVP_VRE_STOR.AnnualSum .= vre_vre_stor * inputs["omega"]
-		dfVP_VRE_STOR = hcat(dfVP_VRE_STOR, DataFrame(vre_vre_stor, :auto))
-		auxNew_Names=[Symbol("Resource");Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
-		rename!(dfVP_VRE_STOR,auxNew_Names)
-		total = DataFrame(["Total" 0 sum(dfVP_VRE_STOR[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
-		total[:, 4:T+3] .= sum(vre_vre_stor, dims = 1)
-		rename!(total,auxNew_Names)
-		dfVP_VRE_STOR = vcat(dfVP_VRE_STOR, total)
-		CSV.write(joinpath(path,"vre_stor_solar_power.csv"), dftranspose(dfVP_VRE_STOR, false), header=false)
+
+		filepath = joinpath(path,"vre_stor_solar_power.csv")
+		if setup["WriteOutputs"] == "annual"
+			write_annual(filepath, dfVP_VRE_STOR)
+		else 	# setup["WriteOutputs"] == "full"
+			write_fulltimeseries(filepath, vre_vre_stor, dfVP_VRE_STOR)
+		end
 	end
+	return nothing
 end