@RolfSander @jimmielin: I've been looking at using the Backwards Euler integrator but I think I might have stumbled upon a bug in int/beuler.f90.

Backward Euler is implemented as an option in the SDIRK integrator. But there is a separate integrator file int/beuler.f90, which is basically a copy of int/sdirk.f90. In the int/sdirk.f90 we have the selection of method based on the value of ICNTRL(3):

!~~~> ICNTRL(3) - method selection       
      SELECT CASE (ICNTRL(3))
      CASE (0,1)
         CALL Sdirk2a
      CASE (2)
         CALL Sdirk2b
      CASE (3)
         CALL Sdirk3a
      CASE (4)
         CALL Sdirk4a
      CASE (5)
         CALL Sdirk4b
      CASE (6)
         CALL BEuler
      CASE DEFAULT
         CALL Sdirk2a
      END SELECT

But we also have this same CASE statement in int/beuler.f90. I would have thought we would want to only have the CALL BEuler routine there and skip calling the other Sdirk* routines.

If you agree, I can open a PR to fix this so that we only call the BEuler routine to set the proper coefficients for the integration.