[mesa-users] Show the evolution of the model at each step of the solver iteration

[Mathieu]

Hi everyone,

I was wondering if there is a way to look at the structure of a model at 
each iteration of the newton solver, before the solution for the next 
timestep is found.
I'd like to make some plots of the structure before the convergence to a 
solution, e.g. T(\rho), M(R), etc...
The purpose of this is to show how the numerical procedure "shakes" the 
star at every timestep to get it to the hydrostatic solution, in order 
to illustrate how the numerics of MESA work.
Of course, I would need it only for a few timesteps, just to have a few 
examples.
So, is there a way to save something like a profile*.data for each 
iteration of the newton solver? In "evolve.f" the profile is saved at 
the end of the step, so I guess I will not be able to get a true profile 
for each iteration. I think what I need is what is stored in the 
structure pointed by x in the routine newton in 
"$MESA_DIR/star/private/star_newton.f":

     subroutine newton( &
          s, nz, nvar, x, xold, sparse, &
          lrd, rpar_decsol, lid, ipar_decsol, which_decsol, &
          tol_correction_norm, &
          xscale, equ, work, lwork, iwork, liwork, AF, &
          lrpar, rpar, lipar, ipar, convergence_failure, ierr)

          type (star_info), pointer :: s
          ! the primary variables
          integer, intent(in) :: nz ! number of zones
          integer, intent(in) :: nvar ! number of variables per zone
          ! the total number of primary variables is neq
          real(dp), pointer, dimension(:) :: x ! =(nvar,nz)
          ! new vector of primaries
          real(dp), pointer, dimension(:) :: xold ! =(nvar,nz)
          ! old vector of primaries

So maybe I could use the routine newton_core_dump defined in star_newton.f .
Am I correct to try to get the values pointed by x if I want to get the 
variables for each zone at each iteration of the solver? Is there a 
function in the run_star_extras which is called at each iteration of the 
solver that I could use?

I appreciate any hint!
Thanks,

Mathieu