diff --git a/desc/vmec.py b/desc/vmec.py
index 8eaffe4d1e..4c6e8a30cb 100644
--- a/desc/vmec.py
+++ b/desc/vmec.py
@@ -271,6 +271,7 @@ def save(cls, eq, path, surfs=128, verbose=1):  # noqa: C901 - FIXME - simplify
             np.array([" " * 100], "S" + str(file.dimensions["dim_00100"].size))
         )  # VMEC input filename: input.[input_extension]
 
+        # TODO: instead of hard-coding for fixed-boundary, allow for free-boundary?
         mgrid_mode = file.createVariable("mgrid_mode", "S1", ("dim_00001",))
         mgrid_mode[:] = stringtochar(
             np.array([""], "S" + str(file.dimensions["dim_00001"].size))
@@ -480,7 +481,7 @@ def save(cls, eq, path, surfs=128, verbose=1):  # noqa: C901 - FIXME - simplify
         chipf.long_name = "d(chi)/ds: poloidal flux derivative"
         chipf[:] = phipf[:] * iotaf[:]
 
-        # scalars on Quadrature grid
+        # scalars computed on a Quadrature grid
         data = eq.compute(
             ["R0/a", "V", "<|B|>_rms", "<beta>_vol", "<beta_pol>_vol", "<beta_tor>_vol"]
         )
@@ -525,7 +526,7 @@ def save(cls, eq, path, surfs=128, verbose=1):  # noqa: C901 - FIXME - simplify
         betator.units = "None"
         betator[:] = data["<beta_tor>_vol"]
 
-        # scalars at last closed flux surface
+        # scalars computed at the last closed flux surface
         grid = LinearGrid(M=eq.M_grid, N=eq.N_grid, rho=[1.0], NFP=NFP)
         data = eq.compute(["G", "I"], grid=grid)
 
@@ -539,7 +540,7 @@ def save(cls, eq, path, surfs=128, verbose=1):  # noqa: C901 - FIXME - simplify
         rbtor.units = "T*m"
         rbtor[:] = data["G"][0]
 
-        # scalars at the magnetic axis
+        # scalars computed at the magnetic axis
         grid = LinearGrid(M=eq.M_grid, N=eq.N_grid, rho=ONAXIS, NFP=NFP)
         data = eq.compute(["G", "p", "R", "<|B|^2>"], grid=grid)