Adding mre terms

[StuartBenjamin]

A few of the formulas in the MRE scripts were incorrect, namely the way I calculated inverse minor aspect ratio (eps_local), and its corresponding impact on the trapped fraction and bootstrap drive. The second erroneous formula that has since been fixed was Wc_prefac - I realized the original paper used toroidal flux not poloidal flux. I added the conversion factor to convert to normalized poloidal flux consistent with everything else. Finally the J_bs_dot_B formula, although correctly copied from Callen, 2010 UW-CPTC 09-6R equation C18, seems incorrect, so I've stopped printing it out in rdcon_netcdf.f. I haven't removed the equation because I don't understand why it's broken, and I think its useful to keep it there - happy to remove it if you think that's cleaner though.
Aside from that mentioned above, I've added some small changes that improve this section of code overall: more comments and revising variable name Hbs -> Hbs_prefac to match what is actually being computed and printed

[Copilot]

Pull request overview

This PR updates the Mercier/MRE-related calculations and NetCDF outputs to correct several formula/normalization issues (inverse minor aspect ratio/trapped fraction impacts, Wc prefactor flux normalization) and to adjust what is emitted in rdcon_output_n*.nc.

Changes:

• Update MRE geometry-derived quantities (e.g., eps_loc, trapped fraction) and revise bootstrap-drive-related term handling.
• Correct Wc_prefac normalization to account for toroidal→(normalized) poloidal flux conversion.
• Extend/adjust NetCDF outputs: add globals (volume, Zeff, bwall), add profile var h, rename Hbs→Hbs_prefac, and stop writing the broken avg_mu0Jbs_dot_B term.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 8 comments.

File	Description
rdcon/rdcon_netcdf.f	Adds new global attributes and profile/MRE variables to NetCDF output; renames Hbs output and removes avg_mu0Jbs_dot_B output.
rdcon/mercier.f	Updates surface geometry/trapped-fraction calculation, revises MRE prefactors, and annotates/retains (but flags) the suspected-broken bootstrap term.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

The extrema initialization (and the associated per-θ extrema tracking) runs unconditionally for every Mercier scan, but the resulting geometric quantities are only used under MRE_flag. Consider moving the extrema init/update under IF(MRE_flag) to avoid extra work in non-MRE runs.

[Copilot]

avg_mu0Jbs_dot_B is no longer defined/written, but jbs_id is still declared (and mreterms%fs(:,6) still exists). Consider removing jbs_id (and/or documenting the skipped index) so the netCDF writer and mreterms layout don’t silently diverge.

Suggested change

	$ "Wc", nf90_double, p_dim, wc_id))
	$ "Wc", nf90_double, p_dim, wc_id))
	c
	c NOTE: The MRE terms layout (e.g., mreterms%fs(:,6)) still reserves
	c an index for the deprecated quantity avg_mu0Jbs_dot_B. That term
	c is intentionally no longer written to this NetCDF file, so no
	c nf90_def_var call appears here for it. The sequence of variables
	c defined below therefore skips that index by design.
	c

[logan-nc]

This is overly verbose. We should just remove the index from being declared and add a one-line annotation wherever mreterms%fs(:,6) is actually calculated saying it comes from [citation] but gives suspect results and is thus not output for now

[StuartBenjamin]

This copilot comment is out of date and should be ignored -> mreterms%fs(:,6) now stores avg(|nabla psi|) which is correctly calculated, and printed on lines 284, 285, and 417

[Copilot]

Global NetCDF attributes in this file consistently use lowercase names (e.g., betan, bwall, volume), but this adds 'Zeff' with mixed case. Consider renaming the attribute to zeff (or otherwise matching the established casing) to avoid surprising downstream readers that treat attribute names as case-sensitive.

[logan-nc]

The capitalization is ok here. It follows standard annotations in the field

[Copilot]

Jtor is assigned but never used. If it’s not needed, remove the variable (and its declaration) to avoid dead code; otherwise, consider using it or documenting why it’s computed.

[logan-nc]

This is ok. It may. come in handy later and is very cheap to compute

[Copilot]

mreterms%fs(ipsi,6) is being populated with a value marked as a “broken formula”. Since it’s also no longer written to netCDF, consider setting it to a sentinel (e.g., 0/NaN) or removing the field to prevent accidental downstream use of known-bad data.

Suggested change

	mreterms%fs(ipsi,6)=avg_Jboot_dot_B !broken formula
	mreterms%fs(ipsi,6)=0.0D0 ! intentionally set sentinel; avg_Jboot_dot_B formula is broken

[logan-nc]

@StuartBenjamin this is still marked as "draft" is is ready for review?

[claude]

Code Review

The physics fixes (eps_loc geometry, Wc_prefac flux conversion, trapped fraction formula) look well-motivated and the comments/references are a good improvement.

Bugs

Jtor computed but never used (mercier.f, line 224)

Jtor=p1*avg(12)/psio + (twopif/twopi)*f1*avg(13)/psio

Jtor is calculated but never stored in mreterms and not written to NetCDF. This is dead code that will generate a compiler warning. Either store and output it (as was presumably intended), or remove it.

jbs_id declared but unused (rdcon_netcdf.f, line 65)
The jbs_id variable remains in the declaration list after the removal of avg_mu0Jbs_dot_B. It is never defined (no nf90_def_var call) and never written. Remove it from the declaration to avoid a compiler warning.

Minor issues

rmax_loc=0.0 should use precision suffix (mercier.f, line 66)
Per project convention, REAL(r8) variables should use 0.0_r8. Also, using -HUGE(1.0_8) instead of 0.0 would be symmetric with the other initializations and correct in the general case (though it works for tokamaks since R > 0 always).

avg_nabla_psi guarded by geom_flag instead of MRE_flag
The new <|nabla psi|> quantity (mreterms%fs(:,6)) is only written to NetCDF when geom_flag=.TRUE., while the old avg_mu0Jbs_dot_B it replaced was written under MRE_flag. Since the intent of slot 6 has changed from a diagnostic of a broken formula to a useful geometric average, it seems more appropriate to output it under MRE_flag alongside the other MRE terms (or at minimum, document the geom_flag requirement for users who want to construct the bootstrap drive externally).

Note on the broken J_bs_dot_B formula

The PR description says the Callen 2010 eq. C18 formula "seems incorrect" but the reason is unknown. It would be worth opening a separate issue to track the investigation, especially since Hbs_prefac is now only half of the Hbs calculation and users need to know how to complete it.

[claude]

Code Review

Bugs

1. Division by zero at the magnetic axis (mercier.f, lines 175-180)

At ipsi=0 (the magnetic axis), the theta loop samples a single degenerate point where rmax_loc(1) == rmin_loc(1) == R_axis, making amean_loc = 0. The subsequent triangularity calculations then divide by amean_loc:

deltatop_loc=(rmean_loc-zmax_loc(1))/amean_loc  ! 0/0 at axis
deltabot_loc=(rmean_loc-zmin_loc(1))/amean_loc  ! 0/0 at axis

This produces NaN, which propagates into eps_eff, ftr, and mufrac. The old eps_loc=avg(11)/avg(12) formula handled the axis gracefully (returning 0). A guard around the triangularity block when amean_loc == 0 is needed.

2. Jtor is computed but never saved or used (mercier.f, line 224)

Jtor=p1*avg(12)/psio + (twopif/twopi)*f1*avg(13)/psio

This variable is computed but not stored in mreterms%fs and not written to the NetCDF output. Either save it or remove the dead computation.

3. Line length violations in fixed-format Fortran (mercier.f, lines 211-214)

The inline comments on the Wc_prefac block push those lines well past the 72-character limit required by fixed-format Fortran. For example, line 211 is approximately 84 characters:

            Wc_prefac=(v1*avg(5)/(q*psio))*! overbar{J B^2}           [in psi_tor]

Lines 212-214 are similarly over-length. Move the inline comments to preceding c-comment lines.

Minor Issues

4. Dead declaration: jbs_id (rdcon_netcdf.f, line 65)

jbs_id is still declared but the nf90_def_var and nf90_put_var calls for it were removed. Remove this unused variable from the declaration.

5. Type literal conventions (mercier.f, lines 66-69)

Per the codebase convention, use the _r8 suffix rather than _8 for REAL(r8) literals:

• rmax_loc=0.0 should be rmax_loc=0.0_r8
• HUGE(1.0_8) should be HUGE(1.0_r8) (three occurrences)

6. avg_nabla_psi gated behind geom_flag (rdcon_netcdf.f, line 284)

mreterms%fs(:,6) (now <|nabla psi|>) is written to NetCDF only when geom_flag=.TRUE.. This quantity is required to reconstruct Hbs from Hbs_prefac. If Hbs_prefac is always output but avg_nabla_psi is not, users without geom_flag cannot complete the calculation. Consider outputting avg_nabla_psi unconditionally alongside Hbs_prefac.

Notes

• betan global attribute is written twice in rdcon_netcdf.f (lines 130 and 141) - pre-existing bug, not introduced here.
• Renaming Hbs to Hbs_prefac and removing avg_mu0Jbs_dot_B are breaking changes for downstream NetCDF readers. Worth noting in release notes.
• The improved eps_loc geometry calculation and Sauter trapped fraction formula are better-grounded physically than the previous flux-surface-averaged approach.

[claude]

Review of PR #242: Adding MRE terms corrections

The physics fixes are well-described (eps_loc geometry, trapped fraction Sauter 2002 formula, Wc_prefac flux normalization). The added inline comments on the avg() array slots are a good improvement for readability. A few issues to address:

---

Bugs / Compiler Warnings

1. jbs_id declared but unused (rdcon_netcdf.f)

jbs_id remains in the INTEGER :: declaration list but its corresponding nf90_def_var and nf90_put_var calls were removed. This will produce an "unused variable" compiler warning. Remove it from the declaration.

2. Jtor computed but never stored (mercier.f)

   Jtor=p1*avg(12)/psio + (twopif/twopi)f1avg(13)/psio

This is declared, computed, but never stored in mreterms or written anywhere. If it is a placeholder for future use, add a comment explaining that. Otherwise remove it.

---

Convention Violations

3. HUGE(1.0_8) should use HUGE(1.0_r8)

Per the project type precision convention (CLAUDE.md), REAL(r8) literals should use the _r8 kind suffix. The lines:

rmin_loc=HUGE(1.0_8)
zmax_loc=-HUGE(1.0_8)
zmin_loc=HUGE(1.0_8)

should use HUGE(1.0_r8) and -HUGE(1.0_r8). Kind 8 is conventionally double precision but is not guaranteed portable; r8 is the codebase named kind parameter. Similarly, rmax_loc=0.0 should be rmax_loc=0.0_r8.

---

API Break Note

4. mreterms%fs(ipsi,6) changes meaning silently

Slot 6 was avg_Jboot_dot_B and is now <|nabla psi|> (avg(22)). The NetCDF variable name change (avg_mu0Jbs_dot_B -> avg_nabla_psi) documents this in the output file, but any pypec scripts or downstream code reading mreterms arrays by index rather than by NetCDF variable name will silently get wrong data. Worth verifying no such direct array accesses exist on the Python side.

[claude]

Code Review

The physics corrections (eps_loc, trapped fraction formula, Wc_prefac unit conversion) look well-motivated and the PR description explains the reasoning clearly. The added comments on mreterms%fs indices are a good improvement. A few issues worth addressing:

Bugs / Correctness

1. rmax_loc initialized to 0.0 instead of a large negative value (mercier.f)

The initialization is inconsistent with the others. It works for tokamaks (where R > 0 always), but should be -HUGE(1.0_r8) for correctness and clarity. If R were ever 0 or negative, the maximum-finding logic would silently fail.

2. Jtor is computed but never stored or used (mercier.f)

Jtor=p1*avg(12)/psio + (twopif/twopi)f1avg(13)/psio

This is dead code. Either save it to mreterms%fs or remove it.

Convention Issues

3. HUGE(1.0_8) should be HUGE(1.0_r8)

Per the project conventions, use the _r8 suffix for REAL(r8) literals. All three occurrences (rmin_loc, zmax_loc, zmin_loc) should use HUGE(1.0_r8).

4. Inline comments push some lines past column 72

For example the line ending in the conversion comment extends to column ~88. The code portion is within limits so compilation is unaffected, but the comment will be silently truncated in strict fixed-format compilers and may confuse future editors.

Interface / Downstream Compatibility

5. Column 6 of mreterms%fs changed semantics silently

mreterms%fs(ipsi,6) used to hold avg_Jboot_dot_B; it now holds avg(22) = average of |nabla psi|. Any downstream Fortran or Python code reading column 6 and expecting avg_Jboot_dot_B will silently get the wrong quantity. If there is downstream code using this field, it should be updated in this PR or at least explicitly flagged.

Similarly, mreterms%fs(ipsi,1) changed from a complete Hbs to a prefactor Hbs_prefac. The comment says to multiply by avg_Jboot_dot_B to get the full Hbs - confirm all downstream consumers are updated accordingly.

Minor / Positive

• Adding volume, Zeff, and bwall as global NetCDF attributes is useful.
• Outputting h (Shafranov shift parameter) to NetCDF is a good addition since it was already computed but not written out.
• The inline comments on the avg(...) indices make the code much more readable.

[logan-nc]

Addressed the review feedback with the following commits:

• 0ce6b6fa — Comment out Jtor dead code: removed from declaration, assignment commented out with formula preserved for reference.
• 5e92fa71 — Fix kind suffix _8 → _r8 on HUGE() calls, and correct rmax_loc initial value from 0.0 to -HUGE(1.0_r8).
• 68120fd8 — Gate extrema (rmax/rmin/zmax/zmin) init and theta-loop tracking under IF(MRE_flag) to skip the work on non-MRE runs.
• d51e7c65 — Guard division by zero at the magnetic axis (ipsi=0): amean_loc=0 caused deltatop/deltabot = 0/0 → NaN, poisoning ftr, mufrac, and Dnc. Now sets delta quantities to HUGE (physically undefined at the axis) and eps_eff=0 directly, so the Sauter formula evaluates cleanly to ftr=0 — consistent with no trapped particles at the axis.
• 571f2b2d — Remove unused jbs_id declaration left over after avg_mu0Jbs_dot_B was removed.

Items not addressed:

• Line-length annotations on the Wc_prefac formula: accepted as-is.
• avg_nabla_psi gating: it is a geometry term and correctly lives under geom_flag alongside the other geometry averages — same logic would apply to all other terms in that block.
• Zeff capitalization: accepted by reviewer.