Dear Users,

1. An almost stupid bug was found to forget to change the seed of random numbers in tracking for each thread. This causes degenerated random numbers in TrackParticles and DynamicPaertureSurvey when NPARA > 1 && RAD. It is temporarily tweaked by adding 2, 4, 6, .. for the seed of each thread, but a more rigorous solution is expected.

> Dear Users,
>
> 1. An almost stupid bug was found to forget to change the seed of random numbers in tracking for each thread. This causes degenerated random numbers in TrackParticles and DynamicPaertureSurvey when NPARA > 1 && RAD. It is temporarily tweaked by adding 2, 4, 6, .. for the seed of each thread, but a more rigorous solution is expected.
>
src/trackd.f, src/tftrack.f内で疑似乱数系列の調整に使っている tfaddseed()関数が、特定の疑似乱数生成器プラグインの実装に依存しています

例えば、SFMT PRNGプラグインだと SeedRandom[]の返すリストの第2要素は、乱数系列の固有識別子になっています
また、数列にシリアル化される PRNGプラグインの内部状態は、実装依存で Checksumや識別用の Magic Numberが
埋め込まれていることも有ります。

個人的には、NPARAが変わると TrackParticles[]の出力の同一性が損なわれるという点は、
計算の再現性という意味で evilだと思います。

例えば、疑似乱数発生を Master Process側で処理して、Shared Memory経由で Worker Processが自分の担当範囲に
対応する疑似乱数を拾うとか…
(Worker Process内で、粒子の並べ替えを行う場合、乱数を拾う場所も並べ替える必要が有りますが)

Dear Users,

1. Another strange behavior was found in the logic of EPS in BEND with RAD. The division number did not necessarily increase/decrease for smaller/larger values opt EPS.

> Dear Users,
>
> 1. Another strange behavior was found in the logic of EPS in BEND with RAD. The division number did not necessarily increase/decrease for smaller/larger values opt EPS.

以下、重箱の隅

この部分のコードで、分割数の上限(旧4/新1024)が設定されていますが、根拠は有るのでしょうか？
tbend()自身に再帰しているので、最終的には ndiv == 1を満足するところまで再帰分割されるので、
call stack depthを考えると一度に分割した方が効率が良いように思えます。

もし、上限を設定するとすれば、Integer型の表現域の上限が根拠になると思います
Fortran95的には変数から、 上限 RADIX(ndiv)*(RADIX(ndiv)**(DIGITS(ndiv)-1) - 1)+1を導出できます
# RADIX(ndiv)**DIGITS(ndiv) - 1が途中の計算でオーバーフローするのを避けるために面倒な表現に…
# コンパイル時に確定する定数なので、まともなコンパイラなら実行時の負荷は増えません

また、上限の制約を int()で変換された後に行いますが、Real型の値
(1 + (abs(phi0) * max(h0/100.d0,1.d0/rphidiv) / eps1))がInteger型の
表現域を越えている場合、どのような値に変換されるか保証はされていなかったと思います
したがって、max()はReal型でとる必要が有ります int(max(dble(ndivmax), ...))
さらに、例外状況を考えると、eps1で除算した項がNaNである可能性が排除出来ないので、
NaNを含む入力に max()が何を返すか保証がなければ、int()にNaNが渡る可能性を考慮する必要があります。

現状の実装で起こりうるケースとしては、EPSが 十分小さな正の値の時、
(... / eps1)が Integer型の表現域を越える大きな値 or +Infinityになり、
int()での変換結果が INT_MAXに張り付かない実装では、意図通りに
再帰分割が働かないことになります。

例えば FreeBSD/amd64環境の gfortran 4.5.4 20110721では、
* int(1.d0 / 1.d-300)は -2147483648になります
* int(NaN)は、-2147483648になります

同様に、BEND B = (EPS = 0/0);で EPSに NaNを設定可能です。
意図的に設定しなくても、何らかの計算結果から EPSパラメータをセットする場合に NaNになる可能性を排除出来ません