src/pec.f90: don't rely on P_inside, I_cd_inside

When combining multiple independent simulations with `gray -s`, the
statistics in sum-summary.txt are computed by summing the profiles
pointwise and calling `postproc_profiles` with their (supposedly common)
MHD equilibrium.

The subroutine assumes the normalisations of dP/dV⋅dV and J_φ⋅dA are
P_inside(ρ=1), I_cd_inside(ρ=1), however if the individual profiles have
been computed with different versions of Gray, built with different
compilers, or on different machines, the values may be very slightly off.
Further, due to the use of the relation Δρ² = ⟨ρ²⟩ - ⟨ρ⟩², the profile
widths are very sensible to the normalisation, producing wildly
incorrect values.

For example, due to the refactor of `magsurf_data`, a change of about
0.05% in dV resulted in Δρ² being off by 300%.

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,   &
+          call postproc_profiles(equil, rhot_tab, dpdv_beam,                    &
-              jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,      &
+              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

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=prof, file=output_dir//'/sum-ec-profiles.txt', action='write')
-    open(newunit=summ, file=output_dir//'/'//'sum-summary.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,  &
+        rhot_tab, results%dPdV, jphi, rhotpav, drhotpav, rhotjava, &
-        rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip,     &
+        drhotjava, dpdvp, jphip, rhotp, drhotp, rhotj, drhotj,     &
-        rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx)
+        dpdvmx, jphimx, ratjamx, ratjbmx)
 
       write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12
       write(summ, trim(fmtstr)) &

src/pec.f90

@@ -244,8 +244,8 @@ contains
   end subroutine pec_tab
 
 
-  subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, &
+  subroutine postproc_profiles(equil, rhot_tab, dpdv, ajphiv, &
-        rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip,                 &
+        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_) :: P_tot, I_tot, ratjplmx, rhopjava, rhoppav
-    real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava
+    real(wp_) :: rhotjav, dvdrhotav, dadrhotava
 
-    rhotpav   = 0
+    ! ⟨ρ_t⟩, Δρ_t using dP/dV⋅dV as a PDF
-    rhot2pav  = 0
+    ! Note: the factor √8 is to match the FWHM of a Gaussian
-    rhotjav   = 0
+    P_tot = sum(dpdv*dvol)
-    rhotjava  = 0
+    if (P_tot > 0) then
-    rhot2java = 0
+      rhotpav  = sum(rhot_tab * dpdv*dvol)/P_tot
-
+      drhotpav = sqrt(8 * sum((rhot_tab - rhotpav)**2 * dpdv*dvol)/P_tot)
-    if (pabs > 0) then
+    else
-      rhotpav  = sum(rhot_tab   *dpdv*dvol)/pabs
+      rhotpav  = 0
-      rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs
+      drhotpav = 0
     end if
 
-    if (abs(currt) > 0) then
+    ! ⟨ρ_t⟩ using J_φ⋅dA as a PDF
-      rhotjav = sum(rhot_tab*ajphiv*darea)/currt
+    I_tot = sum(ajphiv*darea)
-    end if
+    if (abs(I_tot) > 0) then
-    sccsa = sum(abs(ajphiv)*darea)
+      rhotjav = sum(rhot_tab*ajphiv*darea)/I_tot
-    if (sccsa > 0) then
+    else
-      rhotjava  = sum(rhot_tab   *abs(ajphiv)*darea)/sccsa
+      rhotjav = 0
-      rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa
     end if
 
-    ! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile
+    ! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF
-    drhotpav  = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps)))
+    I_tot = sum(abs(ajphiv)*darea)
-    drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps)))
+    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