anchor.diff:      integer, parameter :: tov_volume = 5
   	integer, parameter :: num_tov_variables = 5
       
  -    integer, parameter :: tov_model_handle = 1
  +    integer, parameter :: tov_model_npts = 1
       integer, parameter :: num_tov_ipar = 1
   	
       integer, parameter :: tov_output_step_crust = 1
  
positive.code:     	integer :: ierr, i
      	real(dp) ::  lnP, lgP, xwant, lgx, r, a, m, phi, p, vol
      	real(dp) :: rho, eps, lgRho, dlgRho, lgEps, dlgEps
  	
          ierr = 0
          irtrn = 0
  		xwant = rpar(tov_last_recorded_step) - rpar(tov_output_step_crust)
  		do while (xwant > x)
  			
              lnP = xwant
  			r = interp_y(tov_radius, xwant, rwork_y, iwork_y, ierr) + rpar(tov_core_radius)
  			a = interp_y(tov_baryon, xwant, rwork_y, iwork_y, ierr) + rpar(tov_core_mass)
  
negative.code:         rpar(tov_last_recorded_step) = lnP + rpar(tov_output_step_crust)
          rpar(tov_core_mass) = Mcore
          rpar(tov_core_radius) = Rcore*1.0e5/length_g
          ipar(tov_model_npts) = 1
          npts = ceiling(lnPend - lnP)/rpar(tov_output_step_crust) + 1
          
          call alloc_tov_model(npts, s, ierr)
          if (ierr /= 0) return
          
  		call dop853(n,tov_derivs_crust,lnP,y,lnPend,h,tov_default_max_step_size,tov_default_max_steps, &
  			& rtol,atol,itol, tov_solout_crust, iout, work, lwork, iwork, liwork,  &
  