Isolating the Coriolis term in momentum budget

[suyashbire1]

I am trying to carry out the thickness weighted averaged momentum budget as derived in Young(2012). It requires separate advection and Coriolis terms. So I was wondering if it is a good idea to isolate the Coriolis term from CAu, gKEu and rvxv.

fv = CAu - gKEu - rvxv
fvbrute = f*(v(i,j) + v(i+1,j) + v(i,j+1) + v(i+1,j+1))/4

These two methods give slightly different results. Is this supposed to happen?

These figures show the zonal variation of aforementioned terms at a fixed y (meridional midpoint of the domain) and a fixed zl away from the surface. We're using the Sadourny energy conserving scheme.

[adcroft]

These two methods give slightly different results. Is this supposed to happen?

Yes. The expression f*(v(i,j) + v(i+1,j) + v(i,j+1) + v(i+1,j+1))/4 is not how the linear Coriolis term is discretized. There multiple discretizations (you are probably using SADOURNY_ENERGY?) but in all of them the linear contribution is hidden inside the PV advection term such as ( q(i,j) * ( vh(i+1,j) + vh(i,j) ) + q(i,j-1) * ( vh(i+1,j-1) + vh(i,j-1) ) )/2 where q(i,j) is the absolute vorticity over an area-weighted thickness (see L335 of MOM_CoriolisAdv.F90).

I think the linear term can be found from fv = CAu - rvxv. The term gKEu is the d/dx KE term.

Edit: Corrected mistake in formula above where I had written v instead of vh. It should be vh as now shown.

[suyashbire1]

you are probably using SADOURNY_ENERGY?

Yes.

I think the linear term can be found from fv = CAu - rvxv. The term gKEu is the d/dx KE term.

I think we would still have to use fv = CAu - rvxv - gKEu because rvxv = v(v_x - u_y) and gKEu = - (vv_x + uu_x).

We have a set up where a number of layers outcrop. For example in the following figure, meridional variation of various quantities is shown at a fixed x and for a fixed layer. This particular layer outcrops at around 37N (h goes to zero). CAu, gKEu and rvxv all have sharp kinks at this point.

I think the problem lies with dividing by h.

When we calculate fv as (f/h)*vh like in L513, for outcropping layers, we're potentially multiplying and dividing by a small quantity h. How is vh handled in such cases?

Also, I think the error may arise from averaging of layer thicknesses in L331 from neighbouring T points to Bu points. What happens when, say, two of the T points have zero thickness and other two do not?

[adcroft]

Yes, you are right: CAu = fv + rvxv + gKEu. I had thought the KE term was with the pressure but I'm thinking of another model I once worked. :)

When we calculate fv as (f/h)*vh like in L513, for outcropping layers, we're potentially multiplying and dividing by a small quantity h. How is vh handled in such cases?

vh returned from the continuity solver (layer thickness equation) and is an input to this subroutine but can be assumed to be well-behaved; i.e. zero when upstream layers vanish. Most importantly, if all four h points around the q-point are vanished then the two vh that are guaranteed to vanish.

To make the 1/h in f/h behave so that f/h * vh is finite, we actually calculate f/(h+h_neglect) where h_neglect is so small that if added to an Angstrom in real*8 it doesn't change the value of Angstrom. This is a trick used throughout the code to avoid division by zero. In this case, the only time that f/h would "feel" the value of h_neglect is when it will be multiplied by vanishing (zero) values for vh so the product approaches zero when the layers actually vanish. See line L333 for the h_neglect in the q term.

Also, I think the error may arise from averaging of layer thicknesses in L331 from neighbouring T points to Bu points. What happens when, say, two of the T points have zero thickness and other two do not?

I think L333 answers that? This reason we like SADOURNY75_ENERGY is because it handles these situations well.

If you are calculating fv yourself off-line, then it is possible a discrepancy arises if you are using time-averages. Correlations between 1/h and vh mean that <v> != <1/h><vh>. To make your term balance perfect, you could add the fv diagnostic into the model. Hopefully then CAu = fv + rvxv + gKEu exactly. The other gotcha is that the Coriolis term might be bounded or using Bob's "energy dissipation" method which replaces and bounds the vh's. Check that CORIOLIS_EN_DIS is False.

[suyashbire1]

Thank you for the clarifications. :)

Correlations between 1/h and vh mean that <v> != <1/h><vh>.

I'll need a more frequent output so I'll check this next time I run the model again.

CORIOLIS_EN_DIS is false.

[Hallberg-NOAA]

Hi Suyash,

  Your issue is more than just temporal variability. At each
  timestep, vh is determined by the continuity solver based on an
  positive-definite, upwind biased integral over a timestep of a
  subgrid-scale distribution (usually a parabolic fit) of h. A
  simple 2nd-order averaging of the two adjacent values of h onto
  the v point does not lead to a positive definite continuity
  solver, and hence can not be used. In fact there is no scheme for
  interpolating h onto v that is independent of v and is positive
  definite for h. Within the range of thickness within the stencil
  of the calculation of vh, there is a value of h (call it ~h), such
  that vh = ~h * v, but where you will find ~h depends on v via the
  continuity solver.
  - Bob Hallberg

On 05/04/2016 02:52 PM, Suyash Bire wrote:

  Thank you for the clarifications. :)

    Correlations between 1/h and vh
      mean that <v> != <1/h><vh>.

  I'll need a more frequent output so I'll check this next time I
    run the model again.
  CORIOLIS_EN_DIS is false.
  —
    You are receiving this because you are subscribed to this
    thread.
    Reply to this email directly or view
      it on GitHub