Leith+E bugfix

[iangrooms]

Gustavo Marques noticed that with LEITH_CK = 0.0 we were still getting negative Laplacian viscosity coefficients southwest of land. I.e. we're getting backscatter when we shouldn't be, which is a bug.

This code is designed to limit the magnitude of the backscatter so that it doesn't overcome the biharmonic damping:

              if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) < abs(vort_xy_smooth(i,j))) then
                m_leithy(i,j) = CS%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2
              else
                m_leithy(i,j) = CS%m_leithy_max(i,j)
              endif

Note that CS%m_const_leithy has a factor of CS%c_K in the definition. When c_K = 0.0 and when smoothed vorticity at the q point is zero, this if statement ironically causes m_leithy to equal its largest possible value. (Note that the backscatter coefficient is m_leithy times the biharmonic coefficient.) Smoothed vorticity at the q point is always zero on land, so cells whose northwest corner are on land always have the max backscatter coefficient (which then gets smoothed into the interior). This commit zeros out the backscatter in those cells.

I have run one month of alpha05c with this to verify that it zeros out the offending backscatter coefficient. On Derecho see

/glade/work/igrooms/leithy_bugfix
/glade/derecho/scratch/igrooms/g.e30_a05c.G_JRA.TL319_t232_hycom1_N75.2025.999

In simulations where backscatter is turned on (c_K>0) the current code sets m_leithy to its max value in points southwest of land, which is not intended behavior. We do want backscatter, but not the max value. The bugfix zeros out m_leithy in points southwest of land. This is at least consistent with having small values of backscatter near land, which is what we want. (All the theory behind Leith+E is about eddies, not boundary currents.)

PS, one might be tempted to just change the if statement to

(CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) <= abs(vort_xy_smooth(i,j))

The problem there is that if both sides are equal to 0 then you use an expression that produces a NaN.

[gustavo-marques]

Thanks for the fix!
This PR passes the pr_mom tests.