I have been having issues when running my nested model configuration while using the GLS vertical mixing scheme. The vertical diffusivity occasionally reaches huge values, and eventually blows up the model. Depending on the exact model settings I choose (such as changing the timestep or using ln_zad_Aimp), it can take a few days or a couple of months to blow up, but it always does eventually. The same configuration without AGRIF runs fine, and the same configuration with TKE runs fine.
This always occurs in the first few levels and the same few gridpoints of my child domain, close to the coast and the boundary with the parent domain.
I am using ln_zdfevd = .true., which should set the vertical diffusivity to 10m2/s when a layer is unstable, but the layers where this occurs do not have a negative N^2.
This issue was occurring in NEMO 4.2 and is basically unchanged when I switched to NEMO 5.
I’ve attached a picture of this, you can see the points with huge vertical diffusivity in the southwesternmost corner of the domain. I also tried closing these few points but the problem just moved to adjacent points.
We’ve also found this problem in our nested configuration in NEMO 4.2, the symptom was mxln was becoming too large in a few seemingly random cells in the connection zone on the child grid. Limiting mxln (to be <=1000) helped but there was still too much spurious vertical mixing. The issue occurred with the GLS scheme, and it went away with TKE.
Anyhow, where NEMO transfers u,v,t,s and ssh from parent to child, there is also code that transfers avm. We added code to transfer avt, and this eliminated the problem.
I don’t have a clear understanding about this yet (was avm transferred to fix a similar issue? is this transferring of avm/avt masking a bug that should be fixed elsewhere?) so I think it amounts to a workaround at the moment
Thanks for raising this. I experienced myself GLS generating weird values at the boundaries of a nested domain. At these locations, the turbulence scheme relies on large scale values (shear and stratification) provided by the parent domain, but is still fully prognostic in time. A better solution would be to set turbulence values (kinetic energy and dissipation) from the parent model, but there’s a one step time lag in NEMO for turbulence that prevents to do this correctly.
What is important to consider is that these boundary values should not interfere with the solution inside the domain. The only turbulent value that impacts the nested domain is avm, because it is used to compute the shear production term which is itself given as a spatial average of staggered in space velocity shears. Otherwise, everything is 1d, and located at tracer (T) points.
As a consequence, I don’t understand why setting avt at the boundary would impact the solution in the nested domain. It certainly overwrite “weird” values generated by the GLS scheme at these locations but should not have an impact on the solution elsewhere. Could you check this please ?
This said, if the model crashes because of these local, abnormal values, the trick makes sense.
