Hello Chenyu,
FYI I am not a professor, just a doctor. But you can just call me Chi-Tuan.
Anyway, I have suggested you to proceed step by step to isolate the issue. If the first 3 steps are OK, that means your implementation in Fortran files is suspicious. You should then investigate deeper in what you try.
First classical question: does it work in scalar mode/sequential mode (only 1 core)? If yes and not in parallel, this is a parallel issue. If no, you should start with less nodes at the boundary to test your own implementation. Maybe you can start with a coarser mesh or another mesh just to check your implementation.
Looking very quickly to what you have written, I have first quick remarks:
- please you the latest release, you seem to use a quick old one. There is no assistance for old releases,
- when I read:
IF (NCSIZE.GT.1) K = MESH%KNOLG%I(N)
K = N
I can say that the first line is useful for parallel computations but is useless as the 2nd line always changes the value K to N.
Be sure what you have implemented is exactly what you want.
Hope this helps,
Chi-Tuan