While testing the carboncycle simulation in GEOS-Chem (which uses the Forward-Euler integrator, int/feuler.f90), I was getting negative concentrations that were hard to pin down. Long story short, I traced this to this line of code:

but the generated Fun routine always takes the vector of fixed species concentrations (FIX) as the 2nd argument:

SUBROUTINE Fun ( V, F, RCT, Vdot, Aout, Vdotout )

! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)
! F - Concentrations of fixed species (local)
  REAL(kind=dp) :: F(NFIX)

This led to the F dummy-argument variable being set to a random negative number and propagating negatives through the code.

The fix is just to change

CALL Fun( Y, VAR, RCONST, Ydot ) 

to

CALL Fun( Y, FIX, RCONST, Ydot ) 

in int/feuler.f90. I will submit a PR soon for this for KPP 3.0.0.