diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults
index 399bd83d1..eb00f74a0 100644
--- a/star/defaults/controls.defaults
+++ b/star/defaults/controls.defaults
@@ -542,7 +542,7 @@
       !   Format originally defined for the GONG solar model project.  A
       !   definition was given in 2005 for the `CoRoT/ESTA project`_ and
       !   `GONG itself <https://phys.au.dk/~jcd/solar_models/file-format.pdf>`_.
-      !   MESA's implementation largely follows this subsequent 
+      !   MESA's implementation largely follows this subsequent
       !   `2015 definition <https://phys.au.dk/~jcd/ASTEC_results/file-format.pdf>`_.
 
       ! ``'OSC'``
@@ -673,7 +673,7 @@
       ! format_for_OSC_data
       ! ~~~~~~~~~~~~~~~~~~~
 
-      ! float format for ``'OSC'`` data format 
+      ! float format for ``'OSC'`` data format
 
       ! ::
 
@@ -1338,7 +1338,7 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
       !     envelope_fraction_left = (star_mass - he_core_mass)/(initial_mass - he_core_mass)
-      
+
       ! Stop when ``envelope_fraction_left`` < this limit.
 
       ! ::
@@ -1350,7 +1350,7 @@
       ! ~~~~~~~~~~~~~~~~
 
       !     xmstar = mstar - M_center
-      
+
       ! stop when xmstar in grams is < this.  <= 0 means no limit.
 
       ! ::
@@ -1362,7 +1362,7 @@
       ! ~~~~~~~~~~~~~~~~
 
       !     xmstar = mstar - M_center
-      
+
       ! stop when xmstar in grams is > this.  <= 0 means no limit.
 
       ! ::
@@ -1802,7 +1802,7 @@
       ! fe_core_infall_mass
       ! ~~~~~~~~~~~~~~~~~~~
 
-      ! Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``fe_core_infall`` but there will be a greater chance of a 
+      ! Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``fe_core_infall`` but there will be a greater chance of a
       ! transistent velocity spike causing the model to prematurely exit.  In solar masses
 
 
@@ -1821,7 +1821,7 @@
       ! non_fe_core_infall_mass
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``non_fe_core_infall`` but there will be a greater chance of a 
+      ! Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``non_fe_core_infall`` but there will be a greater chance of a
       ! transistent velocity spike causing the model to prematurely exit. In solar masses
 
 
@@ -1916,8 +1916,8 @@
       ! ::
 
     stop_near_zams = .false.
-    
-    
+
+
       ! stop_at_phase_PreMS
       ! ~~~~~~~~~~~~~~~~~~~
       ! stop_at_phase_ZAMS
@@ -1952,7 +1952,7 @@
       ! ::
 
       ! ::
-    
+
     stop_at_phase_PreMS = .false.
     stop_at_phase_ZAMS = .false.
     stop_at_phase_IAMS = .false.
@@ -1966,7 +1966,7 @@
     stop_at_phase_O_Burn = .false.
     stop_at_phase_Si_Burn = .false.
     stop_at_phase_WDCS = .false.
-    
+
 
       ! Lnuc_div_L_upper_limit
       ! ~~~~~~~~~~~~~~~~~~~~~~
@@ -2256,7 +2256,7 @@
 
       ! If false, then stick to the usual definition -- P/(g*rho).
       ! If true, use min of the usual and sound speed * hydro time scale, sqrt(P/G)/rho.
-      ! Note that the 'TDC' ``MLT_option`` does not respect the ``alt_scale_height`` option, and continues to use ``h = P / rho g``   
+      ! Note that the 'TDC' ``MLT_option`` does not respect the ``alt_scale_height`` option, and continues to use ``h = P / rho g``
       ! even if that flag is set.
 
       ! ::
@@ -2372,7 +2372,6 @@
       ! ::
 
     mlt_make_surface_no_mixing = .false.
-    
 
 
       ! T_mix_limit
@@ -3112,7 +3111,7 @@
     RSP_nz_outer = 40
     RSP_T_anchor = 11d3
     RSP_T_inner = 2d6
-    RSP_testing = .false. 
+    RSP_testing = .false.
 
       ! ::
 
@@ -3359,7 +3358,7 @@
       ! Mixing coefficients are multiplied by this factor.
       ! The ``mix_factor`` is applied in subroutine ``get_convection_sigmas`` in ``star/private/mix_info.f90`` --
       ! the lagrangian diffusion coefficient sigma(k) at cell boundary k is set to
-      ! ``mix_factor*D*(4*pi*r(k)^2*rho_face(k))^2``. 
+      ! ``mix_factor*D*(4*pi*r(k)^2*rho_face(k))^2``.
       ! Note that the value of D is not changed -- it is just used as a term in calculating sigma.
 
       ! ::
@@ -3995,7 +3994,7 @@
       ! use_other_j_for_adjust_J_lost option to specify a specific angular momentum
       ! of removed material different from j_rot_surf
 
-      ! In order to prevent the algorithm from digging to deep to adjust J, 
+      ! In order to prevent the algorithm from digging to deep to adjust J,
       ! there is a timestep limit adjust_J_q_limit
 
       ! ::
@@ -4060,7 +4059,7 @@
 
       ! When w_div_wc_flag is true, rather than a hard limit on w_div_wcrit
       ! we use w_div_wcrit_max2<w_div_wcrit_max to provide a smooth transition.
-      ! In the limit of j_rot->infinity, the resulting w_div_wc will match 
+      ! In the limit of j_rot->infinity, the resulting w_div_wc will match
       ! w_div_wcrit_max, while nothing is done when w_div_wcrit_max<w_div_wcrit_max2
 
       ! ::
@@ -4170,13 +4169,13 @@
       ! setting up outer boundary conditions. We caution that the use of ``'fixed_'`` atmosphere
       ! options might conflict with mlt_option = ``TDC``.
 
-      ! + ``'T_tau'``: 
+      ! + ``'T_tau'``:
       !   set Tsurf and Psurf by solving for the atmosphere structure given a T(tau) relation.
       !   The choice of relation is set by the ``atm_T_tau_relation`` control. See also
       !   the ``atm_T_tau_opacity``, ``atm_T_tau_errtol``, ``atm_T_tau_max_tries`` and
       !   ``atm_T_tau_max_steps`` controls.
       !
-      ! + ``'table'``: 
+      ! + ``'table'``:
       !   set Tsurf and Psurf by interpolating in pre-calculated tables based on model
       !   atmospheres. The choice of table is set by the ``atm_table`` control.
       !   Requires tau_factor = 1, as surface of the model must always attach at the base
@@ -4187,14 +4186,14 @@
       !   Teff and log_g. This is consistent with the assumptions used for table construction:
       !   geometrically thin atmospheres with constant flux.
       !
-      ! + ``'irradiated_grey'``: 
+      ! + ``'irradiated_grey'``:
       !   set Tsurf by solving for the atmosphere structure given the irradiated-grey
       !   T(tau) relation of Guillot, T, and Havel, M., A&A 527, A20 (2011).
       !   See also the ``atm_irradiated_opacity``, ``atm_irradiated_errtol``, ``atm_irradiated_T_eq``,
       !   ``atm_irradiated_T_eq``, ``atm_irradiated_kap_v``, ``atm_irradiated_kap_v_div_kap_th``,
       !   ``atm_irradiated_P_surf`` and ``atm_irradiated_max_tries`` controls.
       !
-      ! + ``'fixed_Teff'`` : 
+      ! + ``'fixed_Teff'`` :
       !    set Tsurf from Eddington T(tau) relation
       !    for current surface tau and Teff = ``atm_fixed_Teff``.
       !    set Psurf = Radiation_Pressure(Tsurf)
@@ -4212,7 +4211,7 @@
       ! + ``'fixed_Psurf_and_Tsurf'`` :
       !    get value of Psurf from control parameter ``atm_fixed_Psurf``.
       !    get value of Tsurf from control parameter ``atm_fixed_Tsurf``.
-      !    see the ``conductive_flame`` test_suite for an example of 
+      !    see the ``conductive_flame`` test_suite for an example of
       !    this boundary condition implemented via the ``other_surface_PT`` hook.
 
       ! ::
@@ -4323,11 +4322,11 @@
       ! + ``'fixed'``:
       !   use a uniform opacity, fixed to the opacity of the outermost cell of the
       !   interior model
-      
+
       ! + ``'iterated'``:
       !   use a uniform opacity, iterated to be consistent with the final Tsurf and
       !   Psurf at the base of the atmosphere.
-      
+
       ! + ``'varying'``:
       !   use a varying opacity consistent with the local T and P throughout the
       !   atmosphere. Involves numerical integration of the hydrostatic equilibrium
@@ -4387,7 +4386,7 @@
       !   model atmospheres (Allard et al. 2001) for Teff < 3000 K.
       !   where no published results are available, the table has been
       !   filled in using integration of the Eddington T(tau) relation
-      
+
       ! + ``'photosphere'``:
       !   use model atmosphere tables for photosphere; range of Z's,
       !   as described in MESA Paper I (Paxton et al. 2011), Sec. 5.3.
@@ -4399,7 +4398,7 @@
       !   and the models by Castelli & Kurucz (2003).  where
       !   neither is available, an entry is generated by
       !   integrating the Eddington T(tau) relation
-      
+
       ! + ``'WD_tau_25'``:
       !   hydrogen atmosphere tables for cool white dwarfs
       !   giving Pgas and T at log10(tau) = 1.4 (tau = 25.11886)
@@ -4408,7 +4407,7 @@
       !   R.D. Rohrmann, L.G. Althaus, and S.O. Kepler,
       !   Lyman α wing absorption in cool white dwarf stars,
       !   Mon. Not. R. Astron. Soc. 411, 781–791 (2011)
-      
+
       ! + ``'DB_WD_tau_25'``:
       !   helium dominated (log(H/He)=-5.0) atmosphere tables for DB
       !   white dwarfs, provided by Odette Toloza and Detlev Koester.
@@ -4631,30 +4630,35 @@
 
       ! If ``.true.``, modify the radiative gradient so that the
       ! temperature profile of the optically thin layers follow the
-      ! T(τ) relation chosen by ``atm_T_tau_relation``.  
+      ! T(τ) relation chosen by ``atm_T_tau_relation``.
       ! ::
 
     use_T_tau_gradr_factor = .false.
 
-! star spot controls
-! =================
 
-      ! The impedence of the surface flux due to magnetic pressure from starspots is parameterized 
-      ! in the style of an atmospheric boundary modification. As first described by 
-      ! Somers et al. (2015; ApJ) https://ui.adsabs.harvard.edu/abs/2015ApJ...807..174S, 
-      ! the degree of "spottiness" on the stellar surface is characterized using two parameters:
+! starspots
+! =========
 
-      ! * SPOTF (hereafter fspot): a coverage fraction, or "spot filling factor" (in the notation of the YREC documentation); and
 
-      ! * SPOTX (hereafter xspot): the temperature contrast between the spotted and unspotted regions: xspot = T_spot/T_photosphere.
+      ! fspot
+      ! ~~~~~~~~~~~
 
-      ! The coverage fraction is set to fspot = 0.34 
-      ! (for consistency with observations of low-mass stars: 
-      ! Cao et al., 2022 https://ui.adsabs.harvard.edu/abs/2022ApJ...924...84C
-      ! and the temperature contrast is set to xspot = 0.85 (also from fits to observations).
+      ! Filling Factor of starspots (in [0,1])
+
+      ! ::
+
+    fspot = 0d0
+
+
+      ! xspot
+      ! ~~~~~~~~~~~
+
+      ! Temperature contrast between the spotted and unspotted regions (in [0,1])
+
+      ! ::
+
+    xspot = 1d0
 
-    fspot = 1.0 
-    xspot = 0.0 
 
 ! mass gain or loss
 ! =================
@@ -4671,10 +4675,9 @@
     mass_change = 0d0
 
 
-
       ! mdot_omega_power
       ! ~~~~~~~~~~~~~~~~
-      
+
       ! Enhanced mass loss due to rotation
       ! as in Heger, Langer, and Woosley, 2000, ApJ, 528:368-396.
 
@@ -7107,8 +7110,6 @@
    dq_D_mix_zero_at_H_C_crossover = -1d0
 
 
-
-
       ! diffusion_min_T_for_radaccel
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! diffusion_max_T_for_radaccel
@@ -7319,12 +7320,12 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
       ! Choice of appropriate option for the WD core mixture:
-      
+
       ! + ``'CO'`` : carbon-oxygen phase separation using the two-component
       !   phase diagram of `Blouin & Daligault (2021a) <https://ui.adsabs.harvard.edu/abs/2021PhRvE.103d3204B/abstract>`_.
       ! + ``'ONe'`` : oxygen-neon phase separation using the two-component
       !   phase diagram of `Blouin & Daligault (2021b) <https://ui.adsabs.harvard.edu/abs/2021ApJ...919...87B/abstract>`_.
-      
+
       ! ::
 
     phase_separation_option = 'CO'
@@ -7355,7 +7356,7 @@
       ! ::
 
     phase_separation_mixing_use_brunt = .true.
-    
+
 
 ! eos controls
 ! ============
@@ -7510,7 +7511,7 @@
       ! or use the faster method by allowing for a small tolerance on the mixture used for the computations of these variables (``'mombarg'``).
 
       ! ::
-      
+
       op_mono_method = 'hu'
 
       ! emesh_data_for_op_mono_path
@@ -7526,7 +7527,7 @@
       ! and then unpack it in place with ``unxz -v OP_mono_master_grid_MESA_emesh.txt.xz``
 
       ! ::
-      
+
     emesh_data_for_op_mono_path = ''      ! '' defaults to $MESA_OP_MONO_MASTER_GRID
 
 
@@ -7733,14 +7734,10 @@
     min_magnitude_brunt_B = -1d99
 
 
-
-
-
 ! structure equations
 ! ===================
 
 
-
       ! energy_eqn_option
       ! ~~~~~~~~~~~~~~~~~
       ! Available options are ``'dedt'`` or ``'eps_grav'``.
@@ -7908,7 +7905,7 @@
       ! ::
 
     velocity_q_upper_bound = 1d99
-    
+
       ! velocity_tau_lower_bound
       ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -8221,9 +8218,9 @@
       ! ::
 
     use_dPrad_dm_form_of_T_gradient_eqn = .false.
-    use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.    
-    
-    
+    use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.
+
+
       ! drag_coefficient
       ! ~~~~~~~~~~~~~~~~
       ! min_q_for_drag
@@ -8236,7 +8233,7 @@
       ! burning in massive stars and AGB stars.
 
       ! Under certain circumstances we will not have drag in the surface k=1 zone
-      ! To force the drag term to be on in the outer zone you must enable one of the 
+      ! To force the drag term to be on in the outer zone you must enable one of the
       ! following surface boundary conditions:
       ! ``use_momentum_outer_BC``, ``use_zero_Pgas_outer_BC``, or ``use_fixed_Psurf_outer_BC``
 
@@ -8254,7 +8251,7 @@
     use_drag_energy = .true.
     drag_coefficient = 0d0
     min_q_for_drag = 0d0
-    
+
 
       ! for hydro comparison tests (e.g., Sedov)
 
@@ -8324,7 +8321,7 @@
     RTI_m_full_boost = 4d0
     RTI_m_no_boost = 5d0
 
-    
+
       ! retry_for_v_above_clight
       ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -8335,7 +8332,7 @@
       ! ::
 
     retry_for_v_above_clight = .true.
-    
+
 
 ! solver controls
 ! ===============
@@ -8542,7 +8539,6 @@
     include_w_in_correction_limits = .true.
 
 
-
       ! max_X_for_conv_timescale
       ! ~~~~~~~~~~~~~~~~~~~~~~~~
       ! min_X_for_conv_timescale
@@ -8557,7 +8553,7 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
       ! ::
-    
+
     max_X_for_conv_timescale = 1d0
     min_X_for_conv_timescale = 0d0
     max_q_for_conv_timescale = 1d0
@@ -8706,7 +8702,7 @@
 
       ! do_normalize_dqs_as_part_of_set_qs
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
+
       ! normalize_dqs destroys bit-for-bit read as inverse of write for models.
       ! ok for create pre ms etc., but not for read model
       ! create_pre_ms calls normalize_dqs even if this flag is false.
@@ -8714,12 +8710,12 @@
       ! ::
 
    do_normalize_dqs_as_part_of_set_qs = .false.
-   
+
 
       ! use_DGESVX_in_bcyclic
       ! use_equilibration_in_DGESVX
       ! report_min_rcond_from_DGESXV
-      
+
       ! FOR DEBUGGING ONLY.   NOT FOR GENERAL USE.
 
       ! ::
@@ -9207,7 +9203,7 @@
       !          'Lnuc_cat'             delta_lgL_nuc_cat_limit
       !          'Lnuc_H'               delta_lgL_H_limit
       !          'Lnuc_He'              delta_lgL_He_limit
-      !          'lgL_power_phot'       delta_lgL_power_photo_limit 
+      !          'lgL_power_phot'       delta_lgL_power_photo_limit
       !          'Lnuc_z'               delta_lgL_z_limit
       !          'bad_X_sum'            (solver found bad mass sum)
       !          'dL/L'                 dL_div_L_limit
@@ -9269,12 +9265,12 @@
       ! time_delta_coeff - smaller forces smaller timesteps giving better time resolution.
       ! multiplier for all real number timestep limits and hard limits.
       ! does not apply to integer valued limits such as
-      
+
       ! + solver_iters_timestep_limit
       ! + burn_steps_limit
       ! + diffusion_steps_limit
       ! + diffusion_iters_limit
-      
+
       ! does not apply to varcontrol_target.
       ! analogous to mesh_delta_coeff for better spatial resolution.
 
@@ -10345,7 +10341,7 @@
     min_lgT_for_lgL_power_photo_limit = 9d0
     lgL_power_photo_burn_min = 10d0
     lgL_power_photo_drop_factor = 10
-    
+
 
       ! limits based on changes at photosphere
 
@@ -10442,7 +10438,7 @@
       ! ``abs_du_div_cs`` > ``min_abs_du_div_cs_for_dt_div_min_dr_div_cs_limit`` and
       ! and
       ! ``abs_u_div_cs`` > ``min_abs_u_div_cs_for_dt_div_min_dr_div_cs_limit``
-      ! 
+      !
       ! allow focus on regions near shock face.
 
       ! ::
@@ -10879,7 +10875,6 @@
     delta_lg_XSi_cntr_hard_limit = -1
 
 
-
       ! delta_XH_cntr_limit
       ! ~~~~~~~~~~~~~~~~~~~
 
@@ -11038,9 +11033,9 @@
       ! ~~~~~~~~~~~~~~~~~~
       ! delta_XSi_drop_only
       ! ~~~~~~~~~~~~~~~~~~~
-         
 
-      ! If true, then only limit drops in abundance. 
+
+      ! If true, then only limit drops in abundance.
 
       ! ::
 
@@ -11121,7 +11116,7 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! hard_limit_for_rel_error_in_energy_conservation
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    
+
       ! ::
 
       !     rel_error_in_energy_conservation = abs(error_in_energy_conservation/total_energy)
@@ -11252,7 +11247,7 @@
     solver_test_partials_write_eos_call_info = .false.
     solver_epsder_struct = 1d-5
     solver_epsder_chem = 1d-5
-    
+
     energy_conservation_dump_model_number = -1
 
 
@@ -11489,7 +11484,7 @@
     use_other_set_pgstar_controls = .false.
     use_other_accreting_state = .false.
     use_other_eval_i_rot = .false.
-    
+
       ! ::
 
     use_other_export_pulse_data = .false.
diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90
index 216092ae8..251379874 100644
--- a/star/private/ctrls_io.f90
+++ b/star/private/ctrls_io.f90
@@ -175,7 +175,7 @@ module ctrls_io
     RSP_relax_initial_model, RSP_trace_RSP_build_model, &
    RSP_GREKM_avg_abs_limit, RSP_GREKM_avg_abs_frac_new, RSP_kap_density_factor, RSP_map_columns_filename, &
    RSP_relax_alfap_before_alfat, RSP_max_outer_dm_tries, RSP_max_inner_scale_tries, RSP_T_anchor_tolerance, &
-    ! mass gain or loss  
+    ! mass gain or loss
     mass_change, mass_change_full_on_dt, mass_change_full_off_dt, trace_dt_control_mass_change, &
     min_wind, max_wind, use_accreted_material_j, accreted_material_j, D_omega_mixing_rate, &
     D_omega_mixing_across_convection_boundary, max_q_for_D_omega_zero_in_convection_region, nu_omega_mixing_rate, &
@@ -451,7 +451,7 @@ module ctrls_io
     delta_dX_div_X_cntr_min, delta_dX_div_X_cntr_max, delta_dX_div_X_cntr_limit, delta_dX_div_X_cntr_hard_limit, &
     delta_dX_div_X_drop_only, delta_lg_XH_drop_only, &
     delta_lg_XHe_drop_only, delta_lg_XC_drop_only, delta_lg_XNe_drop_only, delta_lg_XO_drop_only, delta_lg_XSi_drop_only, &
-    delta_XH_drop_only, delta_XHe_drop_only, delta_XC_drop_only, delta_XNe_drop_only, delta_XO_drop_only, delta_XSi_drop_only, &    
+    delta_XH_drop_only, delta_XHe_drop_only, delta_XC_drop_only, delta_XNe_drop_only, delta_XO_drop_only, delta_XSi_drop_only, &
     delta_lg_XH_cntr_min, delta_lg_XH_cntr_max, delta_lg_XH_cntr_limit, delta_lg_XH_cntr_hard_limit, &
     delta_lg_XHe_cntr_min, delta_lg_XHe_cntr_max, delta_lg_XHe_cntr_limit, delta_lg_XHe_cntr_hard_limit, &
     delta_lg_XC_cntr_min, delta_lg_XC_cntr_max, delta_lg_XC_cntr_limit, delta_lg_XC_cntr_hard_limit, &
@@ -491,7 +491,7 @@ module ctrls_io
 
     use_T_tau_gradr_factor, &
 
-    ! star spots
+    ! starspots
     fspot, xspot, &
     
     ! extra heat near surface to model irradiation
@@ -1288,7 +1288,7 @@ subroutine store_controls(s, ierr)
 
  s% use_T_tau_gradr_factor = use_T_tau_gradr_factor
 
- ! star spots
+ ! starspots
  s% fspot = fspot
  s% xspot = xspot
 
@@ -2973,7 +2973,7 @@ subroutine set_controls_for_writing(s, ierr)
  
  use_T_tau_gradr_factor = s% use_T_tau_gradr_factor
 
- ! star spots
+ ! starspots
  fspot = s% fspot
  xspot = s% xspot
 
@@ -4151,7 +4151,7 @@ subroutine get_control(s, name, val, ierr)
             exit
          end if
          if(is_iostat_end(iostat)) exit
-      end do   
+      end do
 
       if(len_trim(val) == 0 .and. ind==0 ) ierr = -1
 
diff --git a/star/test_suite/starspots/README.rst b/star/test_suite/starspots/README.rst
index 990c0f06a..f7b60e9c1 100644
--- a/star/test_suite/starspots/README.rst
+++ b/star/test_suite/starspots/README.rst
@@ -16,14 +16,14 @@ in the style of an atmospheric boundary modification. As first described by
 `Somers et al. (2015; ApJ) <https://ui.adsabs.harvard.edu/abs/2015ApJ...807..174S>`__, 
 the degree of "spottiness" on the stellar surface is characterized using two parameters:
 
-* SPOTF (hereafter fspot): a coverage fraction, or "spot filling factor" (in the notation of the YREC documentation); and
+* SPOTF (hereafter ``fspot``): a coverage fraction, or "spot filling factor" (in the notation of the YREC documentation); and
 
-* SPOTX (hereafter xspot): the temperature contrast between the spotted and unspotted regions: xspot = T_spot/T_photosphere.
+* SPOTX (hereafter ``xspot``): the temperature contrast between the spotted and unspotted regions: xspot = T_spot/T_photosphere.
 
-The coverage fraction is set to fspot = 0.34 
+The coverage fraction is set to ``fspot = 0.34`` 
 (for consistency with observations of low-mass stars: 
 `Cao et al., 2022 <https://ui.adsabs.harvard.edu/abs/2022ApJ...924...84C>`__)
-and the temperature contrast is set to xspot = 0.85 (also from fits to observations).
+and the temperature contrast is set to ``xspot = 0.85`` (also from fits to observations).
  
 Detailed discussion of this functionality can be found in 
 `MESA Instrument Paper VI: Starspots <https://ui.adsabs.harvard.edu/abs/2023ApJS..265...15J>`__.
diff --git a/star/test_suite/starspots/inlist_starspots b/star/test_suite/starspots/inlist_starspots
index 2fd5fcfe1..ef54ff85d 100644
--- a/star/test_suite/starspots/inlist_starspots
+++ b/star/test_suite/starspots/inlist_starspots
@@ -62,11 +62,11 @@
 
    max_age = 24.63d9 !26d9 !0.8d9
 
-   use_other_mlt_results = .true. !.true.
-   use_other_surface_PT = .true. !.True.
+   use_other_mlt_results = .true.
+   use_other_surface_PT = .true.
 
-   fspot = 0.34d0 ! fspot
-   xspot = 0.85d0 ! xspot
+   fspot = 0.34d0
+   xspot = 0.85d0
 
    atm_option = 'table'
    atm_table = 'photosphere'
@@ -83,7 +83,7 @@
    photo_interval = 100
 
    ! starting specifications
-   initial_mass = 0.9 !1.0 ! in Msun units
+   initial_mass = 0.9 ! in Msun units
    min_timestep_limit=0
    max_model_number = -1 !1500
 
@@ -94,7 +94,7 @@
    Blocker_scaling_factor = 0.2
    max_wind = 1d-3
 
-   mixing_length_alpha = 1.9500 !1.95
+   mixing_length_alpha = 1.9500
    MLT_option = 'Henyey'
 
    do_element_diffusion = .false.
diff --git a/star/test_suite/starspots/src/run_star_extras.f90 b/star/test_suite/starspots/src/run_star_extras.f90
index d448ce3ad..9145d2541 100644
--- a/star/test_suite/starspots/src/run_star_extras.f90
+++ b/star/test_suite/starspots/src/run_star_extras.f90
@@ -27,7 +27,7 @@ module run_star_extras
       use const_def
       use math_lib
       use chem_def !! maybe taking up uncessary space but whatever
-      
+
 
       implicit none
       real(dp) :: fspot, xspot, PB_i, Teff_local
@@ -37,10 +37,10 @@ module run_star_extras
       contains
 
       include "test_suite_extras.inc"
-   
+
 
       ! these routines are called by the standard run_star check_model
-      
+
       subroutine extras_controls(id, ierr)
          integer, intent(in) :: id
          integer, intent(out) :: ierr
@@ -48,7 +48,7 @@ subroutine extras_controls(id, ierr)
          ierr = 0
          call star_ptr(id, s, ierr)
          if (ierr /= 0) return
-         
+
          ! this is the place to set any procedure pointers you want to change
          ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)
 
@@ -69,7 +69,7 @@ subroutine extras_controls(id, ierr)
          s% how_many_extra_history_columns => how_many_extra_history_columns
          s% data_for_extra_history_columns => data_for_extra_history_columns
          s% how_many_extra_profile_columns => how_many_extra_profile_columns
-         s% data_for_extra_profile_columns => data_for_extra_profile_columns  
+         s% data_for_extra_profile_columns => data_for_extra_profile_columns
 
          s% how_many_extra_history_header_items => how_many_extra_history_header_items
          s% data_for_extra_history_header_items => data_for_extra_history_header_items
@@ -92,7 +92,7 @@ subroutine extras_startup(id, restart, ierr)
          if (ierr /= 0) return
          call test_suite_startup(s, restart, ierr)
       end subroutine extras_startup
-      
+
 
       integer function extras_start_step(id)
          integer, intent(in) :: id
@@ -115,9 +115,9 @@ integer function extras_start_step(id)
          ! set PB_i here and then other_mlt will know about it, because 
          ! quantities set here are carried over the course of a timestep, 
          ! AKA all newton iterations
-         
+
          mu_ideal_gas = s% mu(1)  !1.00794d0 ! for hydrogen, 1 gram per mole
-         !write(*,*) 'MESA def, my def: ', s% mu(1), mu_ideal_gas  
+         !write(*,*) 'MESA def, my def: ', s% mu(1), mu_ideal_gas
 
          R2 = pow2(s%R(1))
          Teff_local = pow(s%L(1)/(pi4*boltz_sigma*R2), 0.25d0)
@@ -162,7 +162,7 @@ subroutine YREC_spots_other_mlt_results(id, k, MLT_option, &  ! NOTE: k=0 is a v
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
-         
+
         !------------------------------
          !if (s% star_age >= 10d0) then
          if (.not. s% doing_relax .and. .not. s% doing_first_model_of_run) then
@@ -190,7 +190,7 @@ subroutine starspot_tweak_PT(id, skip_partials, &
             logical, intent(in) :: skip_partials
 
             logical :: need_atm_Psurf, need_atm_Tsurf
-         
+
             real(dp), intent(out) :: &
                   lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
                   lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
@@ -211,7 +211,7 @@ subroutine starspot_tweak_PT(id, skip_partials, &
 
             need_atm_Psurf = .true.
             need_atm_Tsurf = .true. 
-            
+
             alp = 1d0 - fspot + fspot*pow4(xspot)
 
             ! This is the surface-average value for luminosity
@@ -248,7 +248,7 @@ integer function extras_check_model(id)
          ierr = 0
          call star_ptr(id, s, ierr)
          if (ierr /= 0) return
-         extras_check_model = keep_going     
+         extras_check_model = keep_going
 
          if (.false. .and. s% star_mass_h1 < 0.35d0) then
             ! stop when star hydrogen mass drops to specified level
@@ -272,7 +272,7 @@ integer function how_many_extra_history_columns(id)
 
          how_many_extra_history_columns = 0
       end function how_many_extra_history_columns
-      
+
 
       subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
          integer, intent(in) :: id, n
@@ -283,14 +283,14 @@ subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
          ierr = 0
          call star_ptr(id, s, ierr)
          if (ierr /= 0) return
-         
+
          ! note: do NOT add the extras names to history_columns.list
          ! the history_columns.list is only for the built-in history column options.
          ! it must not include the new column names you are adding here.
 
       end subroutine data_for_extra_history_columns
 
-      
+
       integer function how_many_extra_profile_columns(id)
          integer, intent(in) :: id
          integer :: ierr
@@ -300,8 +300,8 @@ integer function how_many_extra_profile_columns(id)
          if (ierr /= 0) return
          how_many_extra_profile_columns = 0
       end function how_many_extra_profile_columns
-      
-      
+
+
       subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr)
          integer, intent(in) :: id, n, nz
          character (len=maxlen_profile_column_name) :: names(n)
@@ -381,8 +381,8 @@ integer function extras_finish_step(id)
          ! by default, indicate where (in the code) MESA terminated
          if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step
       end function extras_finish_step
-      
-      
+
+
       subroutine extras_after_evolve(id, ierr)
          integer, intent(in) :: id
          integer, intent(out) :: ierr
diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc
index cb3381acd..64233afca 100644
--- a/star_data/private/star_controls.inc
+++ b/star_data/private/star_controls.inc
@@ -432,10 +432,8 @@
 
          logical :: use_T_tau_gradr_factor
 
-      ! star spot controls
-
-         real(dp) :: fspot 
-         real(dp) :: xspot
+      ! starspots
+          real(dp) :: fspot, xspot
 
       ! mass gain or loss