Hi all

I'm trying to implement a simulation of the Strait of Gibraltar (South of Spain), with a domain spanning from the westernmost limit of the Iberian peninsula (San Vicente Cape) to approximately half of the  southern Spain Coast length (Cartagena). It is few hundreds of km, hence a quite big domain.
The grid is curvilinear, with spherical coordinates, and resolution is of order of 1-2 km in the large areas and few hundreds of meters in the very Strait.
I interpolated bathymetry from a variety of sources and then I smoothed it by iteratively computing the slope factor, so to guarantee it is never higher than abs(0.2).
The vertical coordinates are sigma, with 30 uneven layers.
My intention is to improve horizontal and vertical resolution, and run it as a decomposition domain, but first I want to verify that it works adequately with only one domain.
The model is forced with astronomic water level conditions at both open boundaries (the western and the eastern one). Yes, I know that the same boundary condition is not recommended, but the domain should be big enough to not be an issue. For both boundaries I defined two conditions, the southern and the northern extremes, with tidal constants of 8 constituents, interpolated by the TPXO9-atlas model. M2 amplitude is of order 1 m and 2 cm at western and eastern boundary, respectively, while the other constituents (only semidiurnal and diurnal) are somehow smaller.
The initial condition is interpolated from the MedSea Reanalysis Model from Copernicus.
In a first instance I only configured sea level and current, nor salinity neither temperature, just to make the most simple approximation.

The run goes very well during the first 6 days, and then suddenly crashes. The more struggling fact is that the problem originates only at the northernmost cell of the western boundary and after 6 days of smooth and expected behavior. At that cell, despite the water level matches very well with the one imposed by the astronomic forcing, the current velocity behaves quite weird. It is not harmonic and, after 6 days it suddenly (with no previous changes of the sea level) starts increasing up to completely unreal values (of order of tens of m/s). Obviously, after few iteration, the CFL condition starts failing, and the model ends crashing.

I tried different configurations of viscosity, roughness and 3D turbulence scheme, but nothing works.
Any idea of what is wrong?

I attach the .bca, and .mdf files and few screenshots of sea level and current of the western boundary. The blue line is the northern cell and the red one is the southern one.

Thank you in advance.
Best regards
Simone