[MITgcm-devel] changes in pkg/layers/layers_fluxcalc.F

[Ryan Abernathey]

Ok, my latest commit just puts the interpolation back to the way it was
before. Long story short: I convinced myself there was a problem where
there was none.

FYI, for everyone watching this thread, the main point of my recent changes
to layers was not to screw up the interpolation but rather to introduce
lots of new features related to water mass transformation. Basically, in
addition to the isopycnal velocities, layers now can compute the diapycnal
velocity using the T and S tendency terms. This feature is still very
experimental and doesn't work perfectly, but it is pretty cool.
Documentation and examples will be created over the summer.

-Ryan



On Fri, Jun 5, 2015 at 11:59 AM, Ryan Abernathey <ryan.abernathey at gmail.com>
wrote:

> Hi,
>
> The original point of these changes was to fix a legitimate problem with
> interpolation that occurred when there was a completely filled cell below.
> The filled cells have tracer values of 0, and using these values in
> interpolation caused spurious values to be produced and the layers
> transports to be assigned incorrectly.
>
> In the process of correcting this, I realized that the interpolation
> coefficients were incorrect in the presence of hFac < 1. This is because
> they are precomputed in layers_init_fixed.F as 1D functions of z,
> neglecting the spatially variables changes in cell thickness.  The code I
> checked in tried to fix both these problems, but in retrospect it actually
> only fixed the first one. hFac = 0  is a singular limit, and you want
> different behavior for 0 < hFac < 1.
>
> I propose the following. For now we just do:
>            mfac = MapFact(kk)
>            IF (hFacS(i,j,k+1,bi,bj).EQ.0.0) THEN
>              mfac = 1.0
>            ENDIF
> This will use the same old interpolation as before (neglecting the
> thickness changes), except when the bottom cell is filled.
>
> If you approve, I will check this in, along with a reversion to the
> previous horizontal interpolation.
>
> For the future, I see two options to fix the vertical interpolation:
> - make MapFact a 3D function to account for spatially variable hFaC (could
> potentially require lots of memory)
> - compute MapFact on the fly inside layers_fluxcalc (additional
> computation cost)
>
> I would be interested in what approach you think is best.
>
> -Ryan
>
>
>
> On Thu, Jun 4, 2015 at 6:34 PM, Jean-Michel Campin <jmc at ocean.mit.edu>
> wrote:
>
>> Hi Ryan,
>>
>> When I changed the minus sign to a plus, I recover the previous
>> pkg/layers output where hFac is identically =1 over the column.
>>
>> And regarding the 2nd point (hFacC weighted average), I was not
>> too much concerned (in the previous implementation) if one (i-1 or i /
>> j-1 or j)
>> of the hFacC is zero since, in this case, the flow (uVel/vVel)
>> and hFacW/hFacS are zero.
>>
>> Cheers,
>> Jean-Michel
>>
>> On Thu, Jun 04, 2015 at 12:11:14PM -0400, Ryan Abernathey wrote:
>> > >
>> > > I run a simple 2-D (y-z) test with the updated pkg/layers and the
>> layers
>> > > output (VH & Hs) are different, even where the bottom is flat and all
>> > > hFac = 1 through the full column.
>> > > I look into your changes in layers_fluxcalc.F and the code line
>> 303-304
>> > > seems wrong to me: When hFacS=1, I expect to find the same weighting
>> as
>> > > before
>> > > (MapFact, 1-MapFact) but instead I am getting (2-MapFact, MapFact-1).
>> > > Can you confirm this (I might have missed something) ?
>> > >
>> >
>> > What you say appears to be true. Currently I have
>> >  mfac = 1.0 - hFacS(i,j,k+1,bi,bj)*(MapFact(kk) - 1.0)
>> > but it seems like instead it should be
>> >  mfac = 1.0 + hFacS(i,j,k+1,bi,bj)*(MapFact(kk) - 1.0)
>> > i.e. a simple sign error.
>> > This seems like such a big mistake that there is no way the code should
>> > have worked at all. But I have been using it successfully for months. I
>> > will have a closer look this afternoon and try to sort out what is
>> going on.
>> >
>> >
>> > > The second thing is that I am not sure that we want to consider an
>> hFacC
>> > > weighted averaged when computing the temperature/salt/rho at the
>> > > velocity location. In the case where the model uses the default
>> > > 2nd order centered  advection scheme, the previous layers averaging
>> (no
>> > > hFacC
>> > > weight) was close to the model advection ; but with the new layers
>> > > averaging
>> > > (with hFacC weight) this is no longer the case.
>> > > In other words, I am not convinced that this hFacC weighted averaged
>> is
>> > > more accurate for pkg/layers diagnostics.
>> > > I am curious to know why your made this choice (hFacC weight).
>> > >
>> >
>> > I see your point here too. My thinking was, when there is topography in
>> the
>> > bordering cell, the previous averaging would simply interpolate using
>> the
>> > zero in that cell, producing a highly spurious value at the velocity
>> point.
>> > So there is no doubt that the previous averaging was wrong in this case.
>> > The hFacC weighting was a way of dealing with this by de-weighting the
>> > neighbor cell, eventually to zero for full topography. But I think you
>> are
>> > saying it would be better to just use the masks here.
>> >
>> >
>> > >
>> > > Cheers,
>> > > Jean-Michel
>> > >
>> > > _______________________________________________
>> > > MITgcm-devel mailing list
>> > > MITgcm-devel at mitgcm.org
>> > > http://mitgcm.org/mailman/listinfo/mitgcm-devel
>> > >
>>
>> > _______________________________________________
>> > MITgcm-devel mailing list
>> > MITgcm-devel at mitgcm.org
>> > http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20150606/bcd357ae/attachment.htm>