It seems that there are 2 potential issues with implementation of biharmonic viscosity when used on cubed-sphere grid:

1. in abstract_scalar_biharmonic_diffusivity_closure.jl, we compute "ζ★ᶠᶠᶜ" (starting at line 86) directly without the special cases for the cubed-sphere corners (as we do for the vorticity), so I think this will not be right at the cube corners. An easy fix could be to use the vorticity function instead.
2. the second issue is related to how we evaluate the horizontal laplacian of u and v, i.e., "∇²hfcᶜ" and "∇²hᶜᶠᶜ", that are used for biharmonic viscosity. Currently, the way these operators are defined in Operators/laplacian_operators.jl is more like a flux-form expression that does not work at the cubed-sphere corners. A way to fix this would be to use more like a vector-invariant analog where the laplacian is evaluated from the vorticity and divergence of the horizontal velocity component.