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

[Ryan Abernathey]

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/20150605/c7bd8d92/attachment.htm>