@rakowsk
happens here:

I skip this writing bit am sure it is not the final solution :)

Hi Dima,
oh, that's a part I have to dive into...
But, as the code just warns, not crash: does it work in the end? As a very first test, just comment the print and proceed.

My very first guess is: The matrix structure of the ILU preconditioner is reused (which saves a lot of setup time), but computing new values would result in fill-in in a position that is not available, thus it has to be neglected. This is not critical - just the convergence would be slightly better with the non-zero value.
However, I may be wrong. It could also be a debug statement that the developer hoped would never be reached... again, the code should not break, but the convergence could be sub-optimal to not sufficient.

Fingers crossed that ignoring / commenting the print is sufficient...

Cheers, Natalja

About the new values, I am wrong. It is just a copy. How can it fail? What does the programmer want to tell us?
To be continued...

@rakowsk
by new values a new preconditioner was meant:
I just followed the instructions from here:

Yes, for the ice model, I expect that

• the nonzero pattern of the sparse matrix depends on the mesh structure and is static
• but the values might change a lot and the preconditioner must be recomputed
  And in FESOM, we did not recompute the preconditioner. Never tested feature -> could be a bug.

yes, it is about sea ice and we need a new preconditioner every ice time step. is it not possible with PARMs?

Ok, found the culprit: pARMS does not like matrix coefficients with value 0. It is not a problem when initially setting up the preconditioner (0 entries are ignored or removed - I did not check), but when re-using the ILU-structure, it fails with the cryptic "jw-1" message.
In the ice model, however, zeros are frequently found on the off-diagonal if the mesh node corresponding to the matrix row is not covered by ice, see subroutine VPbc.

I do not see a simple solution, but I will check if we can make pARMS digest zero off-diagonals, or if divide by zero have to be excluded.

A hack to keep on working: Choose a very small off-diagonal value in VPbc - and check if the result is still meaningful or results in ice at the equator:
where (ssh_stiff%colind_loc(is:ie)==row)
ice_stiff_values(is:ie)=1.0_8
elsewhere
ice_stiff_values(is:ie)= -1.e-6 !0.0_8
end where

Two minor bug fixes to follow, but they do not solve the problem, it's pure cosmetics.

Remark on the hack above: Please choose a "small" value from your ice modeling expertise, the 1.e-6 is just a guess! And I have not checked in this hack - all on your own risk! ;-)

Plus a fix to make it compile with gfortran. Problem: Fortran standard does not allow to have a self-referencing interface.

Less hacky intermediate solution: Skip the preconditioner. pARMS is now augmented by BiCGstab w/o preconditioner.
It is not optimal, but better than to re-use a no-longer-fitting ILU from the previous / first timestep: less and faster iterations.

Added a bug fix: There were no halo exchanges at all in ice_VP.F90
One for the solver right hand side is required for sure, I have added it. Please check if other arrays are missing a halo exchange, too!

Next issue: the values for u_ice, v_ice on no-ice-nodes, with rhs=0 and local identity matrix explicitly set in VPbc, may leave parms with u_ice, v_ice /=0.
In most cases, the absolute values are small. As pARMS is an iterative solver and has no knowledge that there is a trivial local solution, it is not incorrect as such.
I will add a correction after the call to pARMS.