NDSOLVE-Error ("For the method !("IDA"), only machine real code is available. "

I am trying to solve mass transfer equations in an adsorption column to obtain the breakthrough curves. Previously, I asked a question on this and received great help:NDSolve for system of PDEs -Error: fewer dependent variables!

Now, I have changed the equations a little bit to account for the properties of another reactor and material. This is the updated code:

(*From Vitaliy Kaurov for nice display of operators*)
pdConv[f_] := 
 TraditionalForm[
  f /. Derivative[inds__][g_][vars__] :> 
    Apply[Defer[D[g[vars], ##]] &, 
     Transpose[{{vars}, {inds}}] /. {{var_, 0} :> 
        Sequence[], {var_, 1} :> {var}}]]
(*PDEModels/tutorial/HeatTransfer/HeatTransferVerificationTests#
463435833*)
ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, [Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k}, 
  V = If[Velocity === "NoFlow", 
    0, [Rho]*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[VectorQ[a], a = DiagonalMatrix[a]];
  (*Note the-sign in the operator*)
  a = PiecewiseExpand[Piecewise[{{-a, True}}]];
  Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, [Rho]_, Cp_, 
  Velocity_, Source_] := [Rho]*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModel[T, X, k, [Rho], Cp, Velocity, Source]
(*Create Parametric PDE Operators for Fluid and Solid*)
(*Include Small Diffusive Term to Fluid Otherwise MMA Might Complain 
about Being Convection Dominated*)
parmfop = 
  TimeHeatTransferModel[c[t, x], t, {x}, {d}, 1, 
   1, {u}, -([Rho] (1 - [Epsilon])/[Epsilon]) Q];
parmsop = 
  TimeHeatTransferModel[q[t, x], t, {x}, {0}, 1, 1, "NoFlow", Q];
parmfop // pdConv
parmsop // pdConv

eps = 0.43;
qflow = 5* (10^-7); (*m^3/s*)
dint = 0.022;(*m*)
Across = [Pi] * (dint^2)/4;(*m^2*)
usup = qflow/Across;
ueps = usup/eps;
Tgas = 301;(*K*)
Rgas = 8.3145;(*J/mol/K*)
Pgas = 1.02;(*bar*)
yf = 0.2;
c0 = yf*Pgas*101325/(Rgas*Tgas);
KL = 0.226*60;(*1/min*)
rho = 1228.5;
qm = 5.09;(*mol/kg*)
nm = 0.429;
K0 = 4.31* (10^-4);(*1/bar*)
dH = -29.38;(*kJ/mol*)
Keq = (K0*Exp[-dH*1000/(Rgas*Tgas)])/101325;(*1/Pascal*)
tend = 10;(*min*)
xend = 0.171;(*m*)
Qsource = -KL * (q[t, 
     x] - (qm*
      Keq*(c[t, x]*Rgas*
         Tgas)/((1 + (Keq*(c[t, x]*Rgas*Tgas))^nm)^(1/nm))));
pdef = (parmfop == 0) /. {[Epsilon] -> eps, u -> ueps, [Rho] -> rho,
     d -> 1, Q -> Qsource};
pdes = (parmsop == 0) /. {Q -> Qsource};
dc = DirichletCondition[c[t, x] == c0 (1 - Exp[-40 t]), x == 0];
icf = c[0, x] == 0;
ics = q[0, x] == 0;
{cifn, qifn} = 
  NDSolveValue[{pdef, pdes, dc, icf, ics}, {c, q}, {t, 0, tend}, {x, 
    0, xend}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.5}}}];

However, after running the code I receive this error: NDSolveValue::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

I checked the equations and I think the main issue is most probably in the Qsource equation, as the code works fine when I replace it with a simple equation. On the other hand, that equation is obtained by fitting isotherm models to experimental data and should be used as it is.

I really appreciate if you could please help fix this issue.

Thanks