The reaction rate of heterogeneous reaction is 0

MFiX
,
2

A primitive, but highly effective, method of debugging is to insert WRITE statements into your code.

      if ((K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)) <= 0) THEN
         write(*,*) "CASE1"
         DES_RATES(REACTION1) = 0
      else if (x1 < 1d-3) THEN
         write(*,*) "CASE2"
         DES_RATES(REACTION1) = 0
      else
         DES_RATES(REACTION1) = x1**3*s*RO_g(IJK)/w*((a3*(a2+b2+b4)+(a2+b2)*b4)*(K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a3*(b2+b4)+b2*b4)*(K2*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a2*b4)*(K3*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)))
         write(*,*) "CASE3", DES_RATES(REACTION1)
         if (DES_RATES(REACTION1) < 0) THEN
            write(*,*) "NEGATIVE"
            DES_RATES(REACTION1) = 0
         ENDIF

I’m seeing a lot of negative values when I run this.

This

x1**3*s*RO_g(IJK)/w*((a3*(a2+b2+b4)+(a2+b2)*b4)*(K1*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a3*(b2+b4)+b2*b4)*(K2*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2))-(a2*b4)*(K3*x_g(IJK,CO)/mw_g(CO)-x_g(IJK,CO2)/mw_g(CO2)))

is a very complex expression, and I’m not going to try to analyze it. I would suggest break this down into some intermediate results, which you can print out and check.