anchor.diff:             ! Recompute variables
              call compute_variables(s)
   
  +           iteration = iteration + 1
              ! Iterate while T hasn't converged
  -           j = 0
              T_converged = .false.
              do while (.not. T_converged)
                 ! Integrate T
  
positive.code:   ! FIXME
    S_crit = 1.e99_x_precision
    ! FIXME
    delta_S_max = 1e-4
    ! FIXME I love it
    delta_T_max = 1e-4
  
    ! Initial time = 0
    t = 0._x_precision
  
negative.code:               
                ! Recompute variables
                call compute_variables(s)
                if (mod(j, 1000) == 0) then
                   call snapshot(s, iteration, t, 13)
                   print*,'snapshot', iteration, t
                end if
  
                ! print *, j, log10(s%T(50)*state_0%T_0), log10(s%S(50)*state_0%S_0), t
                iteration = iteration + 1
             end do
             ! Output things here
          end if
  