main.f

main.f is YREC's main program (it has no PROGRAM statement; the site lists it as (file level)). It reads the run's settings, sets up the physics tables, and then, for each "kind card" (one run of a job: a starting model and what to do with it), loads the starting model and evolves it model by model. Each evolutionary step does, in order: mass loss, mixing and nuclear burning, rezoning, the iterative solution of the structure equations, rotation, light-element burning, the next time step, and output. Around this loop sit optional drivers: automatic calibration of a solar model or of a star of given $L$ and $T_\mathrm{eff}$. This is the reader's starting point.

Running it (from YREC's user documentation, its Quick Start and User Guide, checked against the code).

• The command. yrec control.nml1 physics.nml2: the first file holds the CONTROL namelist (the job's runs, the starting model, the paths of the tables and of the outputs, how often to write), the second the PHYSICS namelist (the physics options and the numerical controls). PARMIN reads them in that order (the Python wrapper accepts either order, but the program does not); without arguments it opens yrec8.nml1 and yrec8.nml2. The names are opened as given, so relative to the working directory, not to the executable's directory as the User Guide says.
• A relaxation code. YREC needs a converged model to start from (FFIRST, read by STARIN: a model from the library under input/models/, or a previous run's .last file). A job is a sequence of up to 50 runs, the kind cards (NUMRUN), each with its type KINDRN(NK): 1 evolves the model (the rescaling options are ignored); 2 rescales it (RSCALE) and relaxes it NMODLS(NK) times with a zero time step (for the zero-age main sequence or horizontal branch, which have static solutions); 3 rescales it and evolves it (for pre-main-sequence models, which have no static solution). LFIRST(NK) says whether a run starts from the starting model or from the previous run's last model (PARMIN forces it true for the first run). A typical job rescales a library model to the wanted mass, composition, mixing length and physics, then evolves it to an end point; the test suite's solar model, for example, rescales and evolves a deuterium-birthline model for two models (KINDRN(1) = 3), then evolves it to the solar age (KINDRN(2) = 1). Each job has one mass and composition; a grid of stars is many jobs, run by scripts or by the Python wrapper's run_parallel.
• Labels, not inputs. XENV0A(NK) and ZENV0A(NK) are copied to XENV0, ZENV0 at the start of each run, and the code uses them only in the output headers (WRTHEAD, WRTOUT, WRTMOD), unless RSCALE overwrites them with the rescaled values; the composition itself comes from the model, as the User Guide says.
• Stop conditions. A run ends after NMODLS(NK) models, or when the age passes ENDAGE(NK) (years) or the central D, $X$ or $Y$ falls below END_DCEN, END_XCEN, END_YCEN (on evolve runs only: see bug 1 below), or when a model fails (step 2 of the evolution step).
• Outputs, each named in the control namelist: .track (FTRACK, one line per model: WRTHEAD, WRTOUT), .short (FSHORT, the numerical details of each model, and the echo of the namelists), .last (FLAST, the last converged model, in the starting model's format: WRTLST, from WRTOUT after each model, and from main.f before the first step of a He-flash run), .store (FSTOR, model snapshots: PUTSTORE), and with LPULSE the pulsation files .pmod, .penv, .patm (WRTMOD). PARMIN creates the directory of FTRACK if it is missing (the Quick Start still says to create output/ by hand); paths may use the placeholders {YREC_INPUT}, {YREC_START}, {YREC_OUTPUT}, which PARMIN replaces from environment variables.

The loop's controls (from the namelist documentation, Control namelist and Physics namelist).

• The kind cards. NUMRUN is the number of runs in the job, and the kind-card loop runs NK = 1 to NUMRUN; KINDRN(NK), LFIRST(NK), NMODLS(NK) and the stop conditions are given for each card in turn (see "A relaxation code" above).

• The iteration limits. NITER1 is the number of Henyey iterations before the surface boundary condition is checked (level 1); NITER2 the iterations after that check, with the start-of-step burning rates (level 2); NITER3 the iterations between the structure and the nuclear burning rates (level 3; the documentation advises 0 when convergence is difficult); NITER4 the iterations between the structure, rotation and rotational mixing (level 4; see the question below). ITDIF1 repeats the block from level 4 on (1 = no repetition; 2 is recommended for rotational mixing); main.f sets it to 1 when LROT is false.

• The first guess. LNEWS extrapolates the structure from the previous step's rates of change (step 7 of the evolution step), which can help convergence. LNEW0 makes CRRECT compute a new envelope triangle (envelopes integrated at three points around the model's $(\log T_\mathrm{eff}, \log L)$, whose values at the fitting point are fitted by a plane to give the outer boundary condition; SURFBC) in the first iteration of every model, ignoring the old one; rotation needs it, since the rotation is held fixed in the solution, and PARMIN turns it on whenever LROT is set.

• Output switches. LRWSH rewinds the .short file before each model, so that it keeps only the last model's details; LSTPCH, with LSTORE, stores the final model of a run in the .store file (label 110, and in WRTOUT).
• Stop conditions: END_DCEN, END_XCEN and END_YCEN are the central deuterium, hydrogen and helium at which a run stops (the documentation's examples: $Y_c$ = 0.97 for the zero-age horizontal branch, $10^{-4}$ for its end).
• Calibration: LCALS calibrates the initial helium and the mixing length $\alpha$ so that the model matches the Sun's $R$ and $L$ at the solar age (with LCALSOLZX, also the surface $Z/X$, CALSOLZX); LCALST, not in the documentation, calibrates a star of given $L$ and $T_\mathrm{eff}$ (SETSCAL).
• Paths in the control namelist may use the placeholders {YREC_INPUT}, {YREC_START} and {YREC_OUTPUT} (above).
• L_ZERODT_PREMS (new at a8a5950; step 3 of the evolution step) is not yet in the public documentation.

• Inputs: everything comes through PARMIN (the run's parameters, per kind card; it also opens files, returns the names of the table files and prints the banner with the version and git hash that SETVERSION sets) and SETUPS (constants, and the opacity, EOS and atmosphere tables). The starting model is read by STARIN from the unit IFIRST.
• Who calls it: nobody; it is the program.
• Outputs: written by the routines it calls: headers by WRTHEAD, each model by WRTOUT (the output driver: short output ISHORT, track ITRACK, model printouts and so on), the last model of a run by PUTSTORE (unit ISTOR, for continuing a sequence), the model before the first step by WRTLST (unit ILAST, for He-flash runs; WRTOUT calls it too, after each model), pulsation output by PDIST, and the solar calibration's summary to ISNU. main.f now owns some values the output routines used to read from COMMON, and passes them on: the mixing length CMIXL (set from CMIXLA(NK) for each kind card), the mass XMSOL (the starting mass, kept for the output), the mixture's name COMPMIX (from STARIN), and NK.
• Modules: main.f uses 20 modules (the COMMON blocks it shared before), mostly to set flags or keep state across the loop: NMLPAR_MOD (the run's parameters, among them per kind card NUMRUN, NMODLS, XENV0A, ZENV0A, CMIXLA, the stop conditions ENDAGE, END_DCEN, END_XCEN, END_YCEN, the iteration limits NITER1–NITER4, ITDIF1, and the switches LCALS, LCALST, LSNU, LNEWS, LEXCOM, LROT, LINSTB, LKUTHE, LOVSTE, ALPHAE, LPTIME, KTTAU, L_ZERODT_PREMS), RUNLOOP_MOD (IRESCA, RESCAL, LENDAG, per kind card), ENTROP_MOD (the rates of change of $\log T$, $\log P$, $\ln L$, $\log r$ from the previous model; shared only with HPOINT, which moves them to the new grid), OLDMOD_MOD (the model at the start of the step, and HCOMPP), PREVCZ_MOD, LIOV_MOD (light-element burning with overshoot), CTOL_MOD (FCORR), DISK_MOD (SAGE), ACDPTH_MOD (outputs at chosen ages), CENT_MOD and FLUXES_MOD (the central values, the convection zone's base and the neutrino fluxes, for the calibration summary), IOUNITS_MOD (the unit numbers), CONSTANTS_MOD, and DIMENSIONS_MOD, LABEL_MOD, CENV_MOD, SOUND_MOD, OVRTRN_MOD, DEUTER_MOD, OP95INTERP_MOD. Former blocks that are now locals or arguments: ORIGSTART (start-of-step $j$ and composition when iterating between structure and rotation; local arrays HJM0, HCOMPP0), CALSUN (DLDX, DRDX, DLDA, DRDA, LSOL, now arguments of SETCAL and CHKCAL), CALSTAR (ALRI, LSTAR, LPASSR, now arguments of SETSCAL and CHKSCAL), ROTPRT (LPRT0, an argument of GETW), BE7 (an output of ENGEB), VERSION, I2O, PULSE (YRECVER, GITHASH, COMPMIX, XMSOL).

The order of the run.

1. Settings and tables. SETVERSION; IOWR = 9 (line 353↗; unused?: PARMIN sets IOWR = 6 at once, parmin.f, line 656↗); PARMIN, then SETUPS (MHP 8/25: the table file names are now passed as arguments instead of through COMMON).
2. Calibration set-up: SETCAL for a solar calibration (LCALS; it changes the kind cards and sets LSOL false, and LPTIME, KTTAU are saved to be restored at each cycle; DLDX, DRDX, DLDA, DRDA, LSOL start at zero and false, so that the final test of LSOL is defined without a calibration), SETSCAL for a star of given $L$, $T_\mathrm{eff}$ (LCALST).

3. Kind-card loop (DO 200 NK = 1, NUMRUN): set the run's $X$, $Z$, mixing length and envelope options; STARIN reads (or keeps) the starting model and rescales it if asked; stop if LROT is set but $\Omega = 0$ (exit status 3, an input error); disable each stop condition the starting model has already passed, by negating its target, and say so; WRTHEAD; for He-flash runs (LKUTHE) store the model with WRTLST; FINDSH (the H-burning shell and convection-zone edges); HTIMER (the first step); zero the entropy terms, and clear LBLIOLD (PREVCZ_MOD), so that LIBURN does not reuse the previous kind card's end-of-step rates for the surface zone.

4. Model loop (DO 100 MODELN = 1, NMODLS(NK)), the evolution step, below.

5. End of a kind card (label 110): PUTSTORE the last model if asked; for a solar calibration, every third kind card CHKCAL compares $L$, $R$ (and the surface $Z/X$) with the targets and either ends (LSOL) or starts another cycle; the counter ICONV starts at 1 and the calibration gives up, silently, when it reaches 15 (the comment still says 10). The radius comes from $L = 4\pi R^2\sigma T_\mathrm{eff}^4$ (RLL, with BL = $\log(L/L_\odot)$): $$\log\frac{R}{R_\odot} = \tfrac12\Bigl[\log\frac{L}{L_\odot} + \log L_\odot - \log 4\pi - \log\sigma - 4\log T_\mathrm{eff}\Bigr] - \log R_\odot;$$ for LCALST, end when the star is matched (LSTAR, on an even kind card).

6. After the runs (label 250): for a converged solar calibration with LSNU, its summary to ISNU (the iterations and derivatives, the neutrino fluxes, the central values, the initial composition, the convection zone's base CZsurf_R, CZsurf_M and the surface composition, the energy-generation fractions), then STOP.

The evolution step (one pass of the model loop).

1. Output switches: rewind ISHORT if LRWSH; turn on pulsation output on the last model, at chosen ages (AGEOUT, JVS 02/11) or at the final age; PDIST if LPOUT.
2. After a failed model (label 15). If the previous attempt diverged or did not converge (LARGE), call STARIN again. Despite the comment "STARIN called here for timestep cutting," STARIN does not cut the step: with LARGE set it prints "RUN STOPPED" and stops (starin.f, lines 225–232). So a failed model ends the run.

3. Ageing? LAGE is true when MODEL ≥ 0, unless the kind card only relaxes a rescaled model (IRESCA = 2, except for pre-main-sequence models with $\log T_c$ < 6.6, which are aged as a failsafe unless L_ZERODT_PREMS is set). With L_ZERODT_PREMS (a CONTROL namelist switch, default false), a rescale card (KINDRN = 2) is a true zero-time-step rescale for a pre-main-sequence model too; the comment notes that a static solution then exists only if nuclear burning can carry the luminosity (as deuterium burning does at the deuterium birthline), and that the next evolved run should fix its step with SETDT, since HTIMER would otherwise take the stored step, which can be far too long.

4. Mass loss and mixing (if LAGE): MASSLOSS; save the composition in HCOMPP; MIX (convective mixing, nuclear burning through KEMCOM, settling through GRSETT); advance the age: DDAGE = $\Delta t$/CSECYR, DAGE += $10^{-9}$DDAGE (Gyr).

5. Rezoning: HPOINT (not in He-flash runs, which only save the old model for PTIME).
6. Light elements at the start (if LEXCOM): the convection zone's base JENV0, the overshoot distance $\alpha_e H_p$ with $H_p = P r^2/(\rho G m)$ (PSCAL0), and LIRATE88 (J = 2: the start-of-step rates, from the saved start-of-step model).

7. First guess of the structure. If LNEWS and the step is positive, extrapolate each point from the previous model's rates of change: $$\log T \mathrel{+}= \dot{(\log T)}_0\,\Delta t,\quad \log P \mathrel{+}= \dot{(\log P)}_0\,\Delta t,\quad L \mathrel{+}= L\,\dot{(\ln L)}_0\,\Delta t,\quad \log r \mathrel{+}= \dot{(\log r)}_0\,\Delta t;$$ otherwise zero these terms. A non-ageing step is made negative (a relaxation).

8. Structure: four levels of CRRECT (the Henyey iteration), each with its own iteration limit: level 1 (NITER1) with the previous envelope triangle; level 2 (NITER2) checking the surface boundary condition (LSBC); then, after storing the new rates of change $$\dot{(\log T)}_0 = \frac{\Delta\log T}{\Delta t},\quad \dot{(\ln L)}_0 = \frac{2(L - L_\mathrm{old})}{(L + L_\mathrm{old})\,\Delta t},\ \dots$$ for the next step's guess, level 3 (NITER3) and level 4 (NITER4), which repeats the iteration between mixing and structure (its comment: not for models where semiconvection matters; for high-precision work such as solar models). Any divergence (LARGE), or no convergence after level 4, goes back to step 2, which stops the run. Levels 1 and 2 keep the composition MIX gave at step 4; levels 3 and 4 reset it to the start-of-step composition HCOMPP and burn and mix again with the new structure (CRRECT, step 2), so the composition is never among the Henyey unknowns: it is advanced apart from the structure (an operator split; see the question below). With rotation, the block from level 4 through step 9 (level 4, the convergence test, MIXCZ, GETW, FPFT) can be iterated ITDIF1 times, restoring the start-of-step composition HCOMPP and $j$ each time.

9. After convergence: MIXCZ (fully mix the convection zones; G. Somers 6/14 also computes the turnover time there); SAGE = age; for rotating models GETW (the new $\Omega$: angular-momentum transport and wind loss) and FPFT (the distortion terms).

10. FINDSH again.
11. Light-element burning (if LEXCOM, MODEL ≥ 0 and the step is positive, and not rotating with rotational mixing; in that case GETW calls LIRATE88, LIBURN and LIBURN2 itself): CONVEC for the zone's base at the end, PSCAL1, LIRATE88 (J = 1: the end-of-step rates, from the new model), LIBURN.
12. After accretion (LNEWFIT): GETNEWENV restores the envelope fitting point, then FPFT.
13. Next step: HTIMER. When the final age is reached, the step written with the model is the previous step, floored at $10^{-3}$ of the age, DELTS = MAX(DELTS0, 1.0D-3*DAGE*1.0D9*CSECYR) (DAGE in Gyr), instead of the clipped last one, so that the sequence can be continued.
14. MODEL += 1 (unless relaxing); WRTOUT.
15. Stop conditions (MHP 10/24): end the kind card when the final age is reached, or the central D, $X$ or $Y$ falls below its threshold; for LCALST, CHKSCAL on even kind cards from the second model on (on success, LPASSR, it jumps to label 200, which skips PUTSTORE).

Which routines do what, in brief. Settings: SETVERSION, PARMIN, SETUPS. Starting model: STARIN (with RSCALE). Time step: HTIMER. Mass change: MASSLOSS, MDOT. Mixing, burning, settling: MIX, KEMCOM, GRSETT, MIXCZ, CONVEC. Rezoning: HPOINT. Structure: CRRECT (with COEFFT, ENVINT, HSOLVE and the physics: EOS, opacities, energy generation ENGEB). Rotation: GETW, FPFT, MOMI. Light elements: LIRATE88, LIBURN. Output: WRTHEAD, WRTOUT, PUTSTORE, WRTLST, PDIST. Calibration: SETCAL, CHKCAL, SETSCAL, CHKSCAL.

The program is one long routine with GOTOs across loop levels (back to label 15 after a failed model, which stops the run, forward to 110, 200 and 250). The state of the evolution lives in local arrays passed to every routine, and in modules. Some variables are passed before main.f gives them a value and are set by the callees: GL by WRTOUT, LMIX, LOCOND, LATMO, LDERIV by COEFFT through CRRECT. TEFFI is set and never used (line 1120↗; unused?). So is LIMIX (line 975↗; unused?), no longer passed to MIXCZ since KC's 2025-05-31 change. A neutrino-flux table block (LNUTAB, lines 461–550↗) is switched off by a hard-coded .FALSE. (unused?): it writes to unit 76, which nothing opens, passes IU and XH2, which main.f never sets, to ENGEB (IU replaced the loop index I in KC's 2025-05-31 change, lines 506–508↗), and stops at the first shell with $\log T < 6$ and at the first with $T_6 < 5$ (line 492↗, line 534↗), so it assumes that $T$ falls outward. The comment at the top still says to change the array size by a global change of "JSON=2000"; JSON = 5000 is now set once, in DIMENSIONS_MOD.