Here’s the trouble:
! First-order method
IF (INTG_EULER) THEN
WHERE(PARTICLE_STATE(:MAX_PIP) == NORMAL_PARTICLE) &
DES_T_s(:MAX_PIP) = DES_T_s(:MAX_PIP) + &
DTSOLID*(Q_Source(:MAX_PIP)/(PMASS(:MAX_PIP)* &
DES_C_ps(:MAX_PIP)))
(gdb) p intg_euler
$8 = .TRUE.
(gdb) p max_pip
$9 = 4
(gdb) p particle_state
$10 = (1, 2, 0, 0)
(gdb) p des_t_s
$11 = (783, 783, 0, 0)
(gdb) p dtsolid
$12 = 2.4791425728385559e-05
(gdb) p q_source
$13 = (-7.1522768060022806e-05, 0, 0, 0)
(gdb) p pmass
$14 = (3.1101767270538953e-05, 3.1101767270538953e-05, 3.4028234663852886e+38, 3.4028234663852886e+38)
(gdb) p des_c_ps
$15 = (0, 0, 0, 0)
So, it’s division-by-zero because of des_c_ps being in the denominator.
Also note unusually large values for pmass (3.4028234663852886e+38)