[MITgcm-devel] some thoughts on Pressure-Coordinate development

[Jean-Michel Campin]

Hi Martin,

Thanks for the comments. I will try to clarify the points you raised.

Regarding:
> > 1) Implement Mass-coordinate
> currently, mass2rUnit = gravity for p-coords. The conversion factor is rhoConst*gravity in some cases
> (i.e. when ever it is not possible to use the isomorphism, for example in ggl90 or with qsw in apply_forcing_t). 
> Would you have a mass2rUnit(k) = rhoConst*gravity(k) if usingMassCoord?

We should use mass2rUnit wherever we can (instead of explicitly "gravity"), so that when we switch
to Mass-coord it will still be valid, with mass2rUnit = 1 in this case (yes, this is simple).
But when we need to get the depth at some location (i,j,k), e.g., in ggl90, we need something 
else and this was my second bullet point.

> > 2) Estimating the depth in Ocean in P-coord or Mass-coord:
> Not sure if I understand how it would work with phiHyd (as this is a changing variable). 
> I do think, that the full non-linear free surface should be the ???default??? mode. 
> As a linear free surface is already an approximation, is it really important to have the 
> best approximation of depth in this case?

I made several mistakes here and the result is very confusing:
a) there are 2 types of conversion, a local one to get dz (from dp or dm) and an absolute
 one to get the full height position z. This was not clear before.
b) I meant "totPhiHyd" and not phiHyd, and it is used in pkg/atm_phys but only works
  when gravity is constant. In this case, from the full geopotential we could get the 
  full height by dividing by gravity (see: pkg/atm_phys/atm_phys_dyn2phys.F, lines 113-115).
  Obviously, there is more work to do to get something valid for variable gravity(k).
c) you are right that "totPhiHyd" is a changing variable so that if we use it to 
 get the depth it's more complicated. But the other issue is that it does integrate
 some approximations so that the position of the sea-surface (SSH) that we get from
 totPhiHyd (or more precisely "phiHydLow") might be a bit off, and this could be annoying.
d) the "compromise" solution (i.e, a conversion factor that only depends on k) is 
  almost already in place to convert dp (or dm) to dz at level interface rF(k) if
  we use "wUnit2rVel" (= gravity*rhoLoc(k)); but it's not defined for center level position (easy
  to add); and this does not help when we need the full height (one could think of using phiRef/g
  but, currently, it's not more precise that just rC(k)/(rhoConst*g) and not great for Mass-coord 
  neither) so we would need to compute this and store it (as a 1-D vector).

Cheers,
Jean-Michel