anchor.diff:    else
        call FARKSetIin('IMEX', 1, ier)
     end if
  -  call FARKSetIin('ADAPT_METHOD', adapt_method, ier)
  -  call FARKSetRin('ADAPT_CFL', cflfac, ier)
  -  call FARKSetRin('ADAPT_SAFETY', safety, ier)
  -  call FARKSetRin('ADAPT_BIAS', bias, ier)
  -  call FARKSetRin('ADAPT_GROWTH', growth, ier)
  -  call FARKSetRin('ADAPT_LB', hfixed_lb, ier)
  -  call FARKSetRin('ADAPT_UB', hfixed_ub, ier)
  -  call FARKSetRin('ADAPT_K1', k1, ier)
  -  call FARKSetRin('ADAPT_K2', k2, ier)
  -  call FARKSetRin('ADAPT_K3', k3, ier)
  -  call FARKSetRin('ADAPT_ETAMX1', etamx1, ier)
  -  call FARKSetRin('ADAPT_ETAMXF', etamxf, ier)
  -  call FARKSetRin('ADAPT_ETACF', etacf, ier)
  +  call FARKSetIin('ADAPT_METHOD',    adapt_method, ier)
  +  call FARKSetRin('ADAPT_CFL',       cflfac, ier)
  +  call FARKSetRin('ADAPT_SAFETY',    safety, ier)
  +  call FARKSetRin('ADAPT_BIAS',      bias, ier)
  +  call FARKSetRin('ADAPT_GROWTH',    growth, ier)
  +  call FARKSetRin('ADAPT_LB',        hfixed_lb, ier)
  +  call FARKSetRin('ADAPT_UB',        hfixed_ub, ier)
  +  call FARKSetRin('ADAPT_K1',        k1, ier)
  +  call FARKSetRin('ADAPT_K2',        k2, ier)
  +  call FARKSetRin('ADAPT_K3',        k3, ier)
  +  call FARKSetRin('ADAPT_ETAMX1',    etamx1, ier)
  +  call FARKSetRin('ADAPT_ETAMXF',    etamxf, ier)
  +  call FARKSetRin('ADAPT_ETACF',     etacf, ier)
     call FARKSetIin('ADAPT_SMALL_NEF', small_nef, ier)
  -  call FARKSetRin('NEWTON_CRDOWN', crdown, ier)
  -  call FARKSetRin('NEWTON_RDIV', rdiv, ier)
  -  call FARKSetRin('LSETUP_DGMAX', dgmax, ier)
  -  call FARKSetIin('LSETUP_MSBP', msbp, ier)
  -  call FARKSetIin('PREDICT_METHOD', predictor, ier)
  -  call FARKSetIin('MAX_NITERS', maxcor, ier)
  -  call FARKSetRin('NLCONV_COEF',nlscoef, ier)
  -  call FARKSetIin('MAX_NSTEPS', 10000, ier)
  +  call FARKSetRin('NEWTON_CRDOWN',   crdown, ier)
  +  call FARKSetRin('NEWTON_RDIV',     rdiv, ier)
  +  call FARKSetRin('LSETUP_DGMAX',    dgmax, ier)
  +  call FARKSetIin('LSETUP_MSBP',     msbp, ier)
  +  call FARKSetIin('PREDICT_METHOD',  predictor, ier)
  +  call FARKSetIin('MAX_NITERS',      maxcor, ier)
  +  call FARKSetRin('NLCONV_COEF',     nlscoef, ier)
  +  call FARKSetIin('MAX_NSTEPS',      10000, ier)
     
     ! output solver parameters to screen
     call FARKWriteParameters(ier)
   
  -  ! set dense linear solver
  +  ! specify use of dense linear solver
     call FARKDense(NEQ, ier)
     call FARKDenseSetJac(1, ier)
   
  
positive.code:   ipar(1) = NEQ
    allocate(y(neqn), ytrue(neqn), rtols(neqn), atols(neqn))
  
    ! set initial conditions
    call init(neqn, T0, y, ytrue, consis)
  
    ! set tolerances
    rtols = 0.d0
    atols = 0.d0
    call settolerances(neqn, rtols, atols, tolvec)
    if (.not. tolvec) then
       atols = atols(1)
    end if
    rtol = rtols(1)
    
    ! initialize vector module
    call FNVInitS(4, NEQ, ier)
  
    ! initialize ARKode solver
    call FARKMalloc(T0, y, imex, 2, rtol, atols, iout, rout, ipar, rpar, ier)
  
    ! set optional inputs
    if (order /= 0) then
  
negative.code: !-----------------------------------------------------------------
  
  
  subroutine farkdjac(neq,t,y,fy,DJac,h,ipar,rpar,wk1,wk2,wk3,ier)
  !-----------------------------------------------------------------
  ! Jacobian computation routine
  !-----------------------------------------------------------------
  