[MITgcm-devel] exf: u/v-wind interpolation and cubed sphere

[Ed Hill]

On Fri, 2006-06-23 at 11:10 -0400, chris hill wrote:
> Martin,
> 
>   Hello.
> 
>   When is Dimitris with you? He and Benny did the interpolation stuff in 
> exf_ so you could talk to him. It does look weird at 90N/S.
> 
>   Certainly it does not make sense to have velocity at 90N/S on a 
> lat-lon mesh, since directions are undefined at that point. It may make 
> sense to calculate a stream function along with the u/v inputs and 
> interpolate that to the target grid zeta points (in a 3d, unstructured 
> cartesian space if need be) and then get u and v from gradients on the 
> target grid in the region of the input grid pole(s). This will work for 
> non-divergent inputs, not sure how to deal with the integration constant 
> for divergent inputs.
> 
>   BTW - Ed is creating a next generation regrid package that we hope 
> will do this and much more, but maybe not this year (Ed commenst?).

Hi folks,

I agree with Martin and don't see how its possible that the U/V
interpolation code in exf_set_uv.F is correct.  It looks like its trying
to do some sort of extrapolation by using values away from the pole but
that seems fishy and, at best, is likely to introduce a discontinuity.

So, here are some comments:

 - As Chris said, lat-lon grids are degenerate (non-isomorphic 
   to the sphere surface) and one result is that velocities at 
   the poles are undefined.

 - In my boundless spare time I'm working on a new regrid scheme 
   but it won't be ready for a few months.  Sorry...

 - The SCRIP code (which is only meant for scalar quantities) 
   does a coordinate transform (essentially, it rotates the domain) 
   when working within some neighborhood of the poles.  This is 
   very similar to what Chris describes above but he mentions 3D 
   and SCRIP does it in 2D.  Coordinate transforms are OK if you 
   also correctly transform (rotate) the velocity components and 
   if you have a reasonably smooth way of matching values at the 
   edge of where the transform begins/ends (which is not easy).

 - If I understand correctly, Martin wants to employ an approach 
   where he uses the current scheme but ignores values within some 
   small neighborhood of the pole and then fills them in afterward 
   using interpolation.  This is a really quick and simple approach 
   and it should work OK provided that the neighborhood is small 
   enough (and perhaps located inside a tile) so that it doesn't
   completely overlap and fill an exchange region.  If it does 
   overlap (and "fill") an exchange width then it will be more 
   difficult to write an interpolation algorithm in MITgcm that 
   is consistent across face edges.

Oh, and add the usual disclaimer that all "simple" interpolation schemes
inevitably act as non-physical sources/sinks for quantities such mass,
energy, momentum, etc...

Ed

-- 
Edward H. Hill III, PhD
office:  MIT Dept. of EAPS 54-1424;  77 Massachusetts Ave.
             Cambridge, MA 02139-4307
emails:  eh3 at mit.edu                ed at eh3.com
URLs:    http://web.mit.edu/eh3/    http://eh3.com/
phone:   617-253-0098
fax:     617-253-4464