diff --git a/fn_red_tag.m b/fn_red_tag.m
index a9f54e7..dc14709 100644
--- a/fn_red_tag.m
+++ b/fn_red_tag.m
@@ -85,7 +85,7 @@
         
         for dir = 1:3 % Fix assume there are three direction, where direction 3 = nondirectional
             sc_dmg = damage.story{s}.(['qnt_damaged_dir_' num2str(dir)]) .* sc_filt;
-            num_comps = comps.story{s}.(['qty_dir_' num2str(dir)]);
+            num_comps = comps.story(s).(['qty_dir_' num2str(dir)]);
 
             % For each structural system
             structural_systems = unique([damage.comp_ds_table.structural_system; damage.comp_ds_table.structural_system_alt]);
@@ -95,6 +95,7 @@
                 % just ignore system 12 (coupling beams)
                 structural_systems(structural_systems == 12) = [];
             end
+            sys_tag = false(num_reals, length(structural_systems));
             
             for sys = 1:length(structural_systems)
                 ss_filt_ds = damage.comp_ds_table.structural_system' == structural_systems(sys) | damage.comp_ds_table.structural_system_alt' == structural_systems(sys);
@@ -103,10 +104,11 @@
                 % Check damage among each series within this structural system
                 
                 % Initialize so the sub-series is reset for each structural system
-                ser_dmg = [];
-                ser_qty = [];
-                
                 series = unique(damage.comp_ds_table.structural_series_id(ss_filt_ds));
+                
+                ser_dmg = zeros(num_reals, length(series));
+                ser_qty = zeros(num_reals, length(series));
+                
                 for ser = 1:length(series)
                     ser_filt_ds = damage.comp_ds_table.structural_series_id' == series(ser); 
                     ser_filt_comp = comps.comp_table.structural_series_id' == series(ser); 
diff --git a/main_PBEErecovery.m b/main_PBEErecovery.m
index 5be3918..f987a4a 100644
--- a/main_PBEErecovery.m
+++ b/main_PBEErecovery.m
@@ -49,6 +49,24 @@
 import recovery.repair_schedule.main_repair_schedule
 import recovery.functionality.main_functionality
 
+%% Set default red tag options for backward compatibility
+if ~isfield(functionality_options,'red_tag_options') || isempty(functionality_options.red_tag_options)
+    functionality_options.red_tag_options = struct();
+end
+
+if ~isfield(functionality_options.red_tag_options,'ignore_coupling_beam_for_red_tag')
+    functionality_options.red_tag_options.ignore_coupling_beam_for_red_tag = false;
+end
+
+if ~isfield(functionality_options.red_tag_options,'tag_coupling_beams_over_height')
+    functionality_options.red_tag_options.tag_coupling_beams_over_height = false;
+end
+
+%% Calculate Red Tags
+[ damage_consequences.red_tag, damage_consequences.red_tag_impact, damage_consequences.inspection_trigger ] ...
+    = fn_red_tag( functionality_options.calculate_red_tag, damage, building_model.comps,...
+    damage_consequences.simulated_replacement_time, functionality_options.red_tag_options );
+
 %% Combine compoment attributes into recovery filters to expidite recovery assessment
 [damage, tmp_repair_class, damage_consequences] = main_preprocessing(...
     damage.comp_ds_table, damage, repair_time_options, tmp_repair_class, ...