diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 98a6e91..e172a71 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -614,8 +614,8 @@ contains
             results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, &
             jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt)
 
-          call postproc_profiles(equil, pabs_beam, icd_beam, rhot_tab, dpdv_beam,   &
-              jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,      &
+          call postproc_profiles(equil, rhot_tab, dpdv_beam,                    &
+              jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,  &
               rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx)              !  *compute profiles width for current beam
 
           if (results%tables%summary%active) then
diff --git a/src/main.f90 b/src/main.f90
index 8c9a484..ef34ad6 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -279,8 +279,8 @@ contains
     end if
 
     ! Open combined output files
-    open(newunit=prof, file=output_dir//'/'//'sum-ec-profiles.txt', action='write')
-    open(newunit=summ, file=output_dir//'/'//'sum-summary.txt', action='write')
+    open(newunit=prof, file=output_dir//'/sum-ec-profiles.txt', action='write')
+    open(newunit=summ, file=output_dir//'/sum-summary.txt', action='write')
 
     ! Move to the filelist directory
     ! (the filepaths are assumed relative to it)
@@ -358,9 +358,9 @@ contains
 
       ! Recompute the summary values
       call postproc_profiles(equil, &
-        results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi,  &
-        rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,     &
-        rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx)
+        rhot_tab, results%dPdV, jphi, rhotpav, drhotpav, rhotjava, &
+        drhotjava, dpdvp, jphip, rhotp, drhotp, rhotj, drhotj,     &
+        dpdvmx, jphimx, ratjamx, ratjbmx)
 
       write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12
       write(summ, trim(fmtstr)) &
diff --git a/src/pec.f90 b/src/pec.f90
index 6b035b0..b191826 100644
--- a/src/pec.f90
+++ b/src/pec.f90
@@ -244,8 +244,8 @@ contains
   end subroutine pec_tab
 
 
-  subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, &
-        rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip,                 &
+  subroutine postproc_profiles(equil, rhot_tab, dpdv, ajphiv, &
+        rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip,    &
         rhotp,  drhotp,  rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx)
 
     ! Radial average values over power and current density profile
@@ -253,7 +253,6 @@ contains
     use gray_equil,           only : abstract_equil
 
     class(abstract_equil), intent(in)  :: equil
-    real(wp_),             intent(in)  :: pabs,currt
     real(wp_),             intent(in)  :: rhot_tab(:)
     real(wp_),             intent(in)  :: dpdv(:), ajphiv(:)
     real(wp_),             intent(out) :: rhotpav, drhotpav, dpdvp
@@ -262,39 +261,44 @@ contains
     real(wp_),             intent(out) :: rhotjfi, drhotjfi, ajmxfi
     real(wp_),             intent(out) :: ratjamx, ratjbmx
 
-    real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav
-    real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava
+    real(wp_) :: P_tot, I_tot, ratjplmx, rhopjava, rhoppav
+    real(wp_) :: rhotjav, dvdrhotav, dadrhotava
 
-    rhotpav   = 0
-    rhot2pav  = 0
-    rhotjav   = 0
-    rhotjava  = 0
-    rhot2java = 0
-
-    if (pabs > 0) then
-      rhotpav  = sum(rhot_tab   *dpdv*dvol)/pabs
-      rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs
+    ! ⟨ρ_t⟩, Δρ_t using dP/dV⋅dV as a PDF
+    ! Note: the factor √8 is to match the FWHM of a Gaussian
+    P_tot = sum(dpdv*dvol)
+    if (P_tot > 0) then
+      rhotpav  = sum(rhot_tab * dpdv*dvol)/P_tot
+      drhotpav = sqrt(8 * sum((rhot_tab - rhotpav)**2 * dpdv*dvol)/P_tot)
+    else
+      rhotpav  = 0
+      drhotpav = 0
     end if
 
-    if (abs(currt) > 0) then
-      rhotjav = sum(rhot_tab*ajphiv*darea)/currt
-    end if
-    sccsa = sum(abs(ajphiv)*darea)
-    if (sccsa > 0) then
-      rhotjava  = sum(rhot_tab   *abs(ajphiv)*darea)/sccsa
-      rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa
+    ! ⟨ρ_t⟩ using J_φ⋅dA as a PDF
+    I_tot = sum(ajphiv*darea)
+    if (abs(I_tot) > 0) then
+      rhotjav = sum(rhot_tab*ajphiv*darea)/I_tot
+    else
+      rhotjav = 0
     end if
 
-    ! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
-    drhotpav  = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps)))
-    drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps)))
+    ! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF
+    I_tot = sum(abs(ajphiv)*darea)
+    if (I_tot > 0) then
+      rhotjava  = sum(rhot_tab * abs(ajphiv)*darea)/I_tot
+      drhotjava = sqrt(8* sum((rhot_tab - rhotjava)**2 * abs(ajphiv)*darea)/I_tot)
+    else
+      rhotjava  = 0
+      drhotjava = 0
+    end if
 
     rhoppav  = equil%tor2pol(rhotpav)
     rhopjava = equil%tor2pol(rhotjava)
 
-    if (pabs > 0) then
+    if (P_tot > 0) then
       call equil%flux_average(rhoppav, dVdrho_t=dvdrhotav)
-      dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav)
+      dpdvp = P_tot * 2/(sqrt(pi)*drhotpav*dvdrhotav)
       call profwidth(rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
     else
       dpdvp  = 0
@@ -303,10 +307,10 @@ contains
       drhotp = 0
     end if
 
-    if (sccsa > 0) then
+    if (I_tot > 0) then
       call equil%flux_average(rhopjava, dAdrho_t=dadrhotava, ratio_astra_tor=ratjamx, &
                               ratio_jintrac_tor=ratjbmx, ratio_paral_tor=ratjplmx)
-      ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava)
+      ajphip = I_tot * 2/(sqrt(pi)*drhotjava*dadrhotava)
       call profwidth(rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
     else
       ajphip   = 0