anchor.diff:      integer, parameter :: num_deriv_rpar = 0
       
   contains
  -
       subroutine do_integrate_crust(NScool_id,ierr)
           use iso_fortran_env, only : error_unit
           use constants_def
  
positive.code:         end if
          
          itol = 0
          ! use the core temperature as a reference for relative accuracy
          rtol = s% integration_tolerance
          atol = s% integration_tolerance * minval(z)
          
  
negative.code:         nullify(ipar)
          nullify(rpar)
  
      end subroutine do_integrate_crust
  
      subroutine evaluate_timestep(nr, xold, x, n, y, rwork_y, iwork_y, interp_y, lrpar, rpar, lipar, ipar, irtrn)
         use iso_fortran_env, only : output_unit, error_unit
         use NScool_terminal, only : do_write_terminal
         use NScool_history, only : do_write_history
         use NScool_profile, only : do_write_profile
         integer, intent(in) :: nr, n, lrpar, lipar
         real(dp), intent(in) :: xold, x
         real(dp), intent(inout) :: y(:)
         ! y can be modified if necessary to keep it in valid range of possible solutions.
         real(dp), intent(inout), target :: rwork_y(*)
         integer, intent(inout), target :: iwork_y(*)
         integer, intent(inout), pointer :: ipar(:) ! (lipar)
         real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
         interface
            include 'num_interp_y.dek'
         end interface
         integer, intent(out) :: irtrn ! < 0 causes solver to return to calling program.
         type(NScool_info), pointer :: s
         integer :: ierr
         logical :: print_terminal_header
         character(len=256) :: filename      
    
         irtrn = 0
         ierr = 0
  
         ! get the crust information information
         call get_NScool_info_ptr(ipar(i_id), s, ierr)
         if (ierr /= 0) then
            write (error_unit,*) 'unable to access NScool info in get_derivatives'
            irtrn = -1
            return
         end if
         if (n /= s% nz) then
            write (error_unit,*) 'wrong number of equations in solver'
            irtrn = -2
            return
         end if
  
         s% model = nr + s% starting_number_for_profile
         s% tsec = x + s% epoch_start_time
         s% dt = x-xold
    
         s% lnT(1:n) = y(1:n)
         s% T(1:n) = exp(s% lnT(1:n))
         call interpolate_temps(s)
    
         call get_coefficients(s,ierr)
         if (ierr /= 0) then
            write (error_unit,*) 'error while interpolating coefficients'
            irtrn = -3
            return
         end if
    
         call evaluate_luminosity(s, ierr)
         if (ierr /= 0) then
            write (error_unit,*) 'error while evaluating luminosity'
            irtrn = -3
            return
         end if
  
         ! update terminal information
         if (mod(s% model,s% write_interval_for_terminal) == 0) then
            if (mod(ipar(i_num_terminal_writes),s% write_interval_for_terminal_header) == 0) then
               print_terminal_header = .TRUE.
            else
               print_terminal_header = .FALSE.
            end if
            call do_write_terminal(ipar(i_id), ierr, print_terminal_header)
            if (ierr /= 0) then
               write(error_unit,*) 'failure writing to terminal'
               return
            end if
            ipar(i_num_terminal_writes) = ipar(i_num_terminal_writes) + 1
         end if
    
         ! update history log
         if (mod(s% model, s% write_interval_for_history) == 0) then
            call do_write_history(ipar(i_id),ierr)
            if (ierr /= 0) then
               write (error_unit,*) 'failure writing history log'
               return
            end if
            write (error_unit,'(a,i0,a)') 'saving model ',s% model,' to history log'
         end if
    
         ! update profile log
         if (mod(s% model, s% write_interval_for_profile) == 0) then
            write(filename,'(a,i4.4)') 'profile',s% model
            call do_write_profile(ipar(i_id),ierr)
            if (ierr /= 0) then
               write (error_unit,*) 'failure writing profile log'
               return
            end if
            write (error_unit,'(a,i0,a)') 'saving model ',s% model,' to profile log'
         end if
      end subroutine evaluate_timestep
  
      subroutine get_derivatives(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
          use const_def, only: dp
          integer, intent(in) :: n, lrpar, lipar
          real(dp), intent(in) :: x, h
  