diff --git a/src/initialization/init_device.jl b/src/initialization/init_device.jl
index e22c51b8a..31144f771 100644
--- a/src/initialization/init_device.jl
+++ b/src/initialization/init_device.jl
@@ -571,6 +571,7 @@ function initialize_dynamic_device!(
     end
 
     set_P_ref(dynamic_wrapper, Pref)
+    PSY.set_P_ref!(dynamic_device, Pref)
     set_Q_ref(dynamic_wrapper, Qref)
     set_V_ref(dynamic_wrapper, Vref)
     set_ω_ref(dynamic_wrapper, Freq_ref)
diff --git a/src/models/device.jl b/src/models/device.jl
index 289fd6960..09b764e48 100644
--- a/src/models/device.jl
+++ b/src/models/device.jl
@@ -1015,6 +1015,7 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     sys_ω = global_vars[GLOBAL_VAR_SYS_FREQ_INDEX]
+    Sbase = get_system_base_power(dynamic_device)
     Vt = sqrt(voltage_r^2 + voltage_i^2)
 
     #Obtain References (from wrapper and device)
@@ -1050,6 +1051,8 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     rrpwr = PSY.get_rrpwr(get_device(dynamic_device))
     Tv = PSY.get_Tv(get_device(dynamic_device))
     Iq_lim = PSY.get_Iq_lim(get_device(dynamic_device))
+    basepower = PSY.get_base_power(get_device(dynamic_device))
+    base_power_ratio = basepower / Sbase
 
     #STATE Vmeas
     _, dVmeas_dt = low_pass_mass_matrix(Vt, Vmeas, 1.0, T_rv)
@@ -1115,8 +1118,8 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     Iq_neg = -Iq
     I_r = real(complex(Ip_limited, Iq_neg) * exp(im * θ))
     I_i = imag(complex(Ip_limited, Iq_neg) * exp(im * θ))
-    current_r[1] = I_r
-    current_i[1] = I_i
+    current_r[1] = I_r * base_power_ratio
+    current_i[1] = I_i * base_power_ratio
 end
 
 #Freq_Flag = 1
@@ -1134,6 +1137,7 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     t,
 ) where {T <: ACCEPTED_REAL_TYPES}
     sys_ω = global_vars[GLOBAL_VAR_SYS_FREQ_INDEX]
+    Sbase = get_system_base_power(dynamic_device)
     Vt = sqrt(voltage_r^2 + voltage_i^2)
 
     #Obtain References (from wrapper and device)
@@ -1182,6 +1186,8 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     rrpwr = PSY.get_rrpwr(get_device(dynamic_device))
     Tv = PSY.get_Tv(get_device(dynamic_device))
     Iq_lim = PSY.get_Iq_lim(get_device(dynamic_device))
+    basepower = PSY.get_base_power(get_device(dynamic_device))
+    base_power_ratio = basepower / Sbase
 
     #STATE Vmeas
     _, dVmeas_dt = low_pass_mass_matrix(Vt, Vmeas, 1.0, T_rv)
@@ -1273,6 +1279,6 @@ function _mdl_ode_AggregateDistributedGenerationA!(
     Iq_neg = -Iq
     I_r = real(complex(Ip_limited, Iq_neg) * exp(im * θ))
     I_i = imag(complex(Ip_limited, Iq_neg) * exp(im * θ))
-    current_r[1] = I_r
-    current_i[1] = I_i
+    current_r[1] = I_r * base_power_ratio
+    current_i[1] = I_i * base_power_ratio
 end
diff --git a/src/post_processing/post_proc_generator.jl b/src/post_processing/post_proc_generator.jl
index ecf3adecf..e28fc73f8 100644
--- a/src/post_processing/post_proc_generator.jl
+++ b/src/post_processing/post_proc_generator.jl
@@ -1,3 +1,22 @@
+"""
+Default function to compute output current. 
+Returns an error
+
+"""
+function compute_output_current(
+    ::SimulationResults,
+    dynamic_device::I,
+    ::Vector{Float64},
+    ::Vector{Float64},
+    ::Union{Nothing, Float64},
+) where {I <: PSY.DynamicInjection}
+
+    #Return error
+    error(
+        "Output current for device type $(typeof(dynamic_device)) is not implemented yet.",
+    )
+end
+
 """
 Function to obtain the output current time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
 the dynamic device and bus voltage. It is dispatched for device type to compute the specific current.
@@ -30,6 +49,99 @@ function compute_output_current(
     )
 end
 
+"""
+Function to obtain the output current time series of a AggregateDistributedGenerationA (DERA) model out of the DAE Solution.
+It receives the simulation inputs, the dynamic device and bus voltage.
+
+"""
+function compute_output_current(
+    res::SimulationResults,
+    dynamic_device::PSY.AggregateDistributedGenerationA,
+    V_R::Vector{Float64},
+    V_I::Vector{Float64},
+    dt::Union{Nothing, Float64},
+)
+
+    #Obtain Data
+    sys = get_system(res)
+    Sbase = PSY.get_base_power(sys)
+    basepower = PSY.get_base_power(dynamic_device)
+    base_power_ratio = basepower / Sbase
+    Freq_Flag = PSY.get_Freq_Flag(dynamic_device)
+    name = PSY.get_name(dynamic_device)
+    if Freq_Flag == 1
+        _, Pord = post_proc_state_series(res, (name, :Pord), dt)
+        _, dPord = post_proc_state_series(res, (name, :dPord), dt)
+        Tpord = PSY.get_Tpord(dynamic_device)
+        P_lim = PSY.get_P_lim(dynamic_device)
+    end
+
+    # Get states
+    θ = atan.(V_I ./ V_R)
+    ts, Ip = post_proc_state_series(res, (name, :Ip), dt)
+    ts, Iq = post_proc_state_series(res, (name, :Iq), dt)
+    _, Mult = post_proc_state_series(res, (name, :Mult), dt)
+    _, Vmeas = post_proc_state_series(res, (name, :Vmeas), dt)
+    Iq_neg = -Iq
+    Ip_cmd = Ip
+    Iq_cmd = Iq
+
+    # Get Params
+    Tg = PSY.get_Tg(dynamic_device)
+    rrpwr = PSY.get_rrpwr(dynamic_device)
+    P_ref = PSY.get_P_ref(dynamic_device)
+
+    I_R = similar(Ip)
+    I_I = similar(Iq)
+    for (ix, Ip_cmd_val) in enumerate(Ip_cmd)
+        Ip_min, Ip_max, _, _ = current_limit_logic(dynamic_device, Ip_cmd_val, Iq_cmd[ix])
+        if Ip[ix] >= 0
+            Rup = abs(rrpwr)
+            Rdown = -Inf
+        else
+            Rdown = -abs(rrpwr)
+            Rup = Inf
+        end
+        if Freq_Flag == 0
+            Ip_input = clamp(P_ref / max(Vmeas[ix], 0.01), Ip_min, Ip_max) * Mult[ix]
+            Ip_limited, _ = low_pass_nonwindup_ramp_limits(
+                Ip_input,
+                Ip[ix],
+                1.0,
+                Tg,
+                -Inf,
+                Inf,
+                Rdown,
+                Rup,
+            )
+        else
+            Pord_limited, _ = low_pass_nonwindup_mass_matrix(
+                dPord[ix],
+                Pord[ix],
+                1.0,
+                Tpord,
+                P_lim[:min],
+                P_lim[:max],
+            )
+            Ip_input = clamp(Pord_limited / max(Vmeas[ix], 0.01), Ip_min, Ip_max) * Mult[ix]
+            Ip_limited, _ = low_pass_nonwindup_ramp_limits(
+                Ip_input,
+                Ip[ix],
+                1.0,
+                Tg,
+                -Inf,
+                Inf,
+                Rdown,
+                Rup,
+            )
+        end
+        I_R[ix] = real(complex(Ip_limited, Iq_neg[ix]) * exp(im * θ[ix])) * base_power_ratio
+        I_I[ix] = imag(complex(Ip_limited, Iq_neg[ix]) * exp(im * θ[ix])) * base_power_ratio
+    end
+
+    return ts, I_R, I_I
+end
+
 """
 Function to obtain the field current time series of a Dynamic Generator model out of the DAE Solution. It receives the simulation inputs,
 the dynamic device and bus voltage. It is dispatched for device type to compute the specific current.
diff --git a/test/test_case42_dera.jl b/test/test_case42_dera.jl
index 0f35588da..3d1f4de88 100644
--- a/test/test_case42_dera.jl
+++ b/test/test_case42_dera.jl
@@ -79,6 +79,7 @@ function test_dera_residual(freqflag_value, csv_file, init_cond, eigs_value)
         t_psid = get_voltage_angle_series(results, 102)[1]
         θ_psid = get_voltage_angle_series(results, 102)[2]
         V_psid = get_voltage_magnitude_series(results, 102)[2]
+        power = get_activepower_series(results, "generator-102-1")
 
         M = get_csv_data(csv_file)
         t_psse, V_psse = clean_extra_timestep!(M[:, 1], M[:, 2])
@@ -143,6 +144,7 @@ function test_dera_massmatrix(freqflag_value, csv_file, init_cond, eigs_value)
         t_psid = get_voltage_angle_series(results, 102)[1]
         θ_psid = get_voltage_angle_series(results, 102)[2]
         V_psid = get_voltage_magnitude_series(results, 102)[2]
+        power = get_activepower_series(results, "generator-102-1")
 
         M = get_csv_data(csv_file)
         t_psse, V_psse = clean_extra_timestep!(M[:, 1], M[:, 2])