Add an optional side drag at walls

[angus-g]

This adds a SIDE_DRAG parameter which applies a stress proportional to velocity at the wall, through the drag coefficient SIDE_DRAG_CD. I've been using this in idealised experiments with flat topography to confirm some WBC separation results with varying lateral boundary conditions. This scheme is disabled by default so there should be no answer changes.

[angus-g]

Whoops, this doesn't build because I got the rebase a little bit wrong.

[angus-g]

I'm not sure how much this overlaps with BOTTOMDRAGLAW either. It seemed to make more sense to do it this way for what I was doing.

[Hallberg-NOAA]

Hi Angus,

    Thank  you for this contribution. I can see it as being very
  valuable, especially for certain idealized configurations.

    This would seem to have significant overlap with the bottom drag
  that is imposed on layers that would intersect with the bottom
  within a grid cell via the existing CHANNEL_DRAG option. The real
  ocean has no sides, so CHANNEL_DRAG can operate even when you are
  not a a coastline, and because CHANNEL_DRAG is imposed implicitly,
  there are no stability limits as there could be for a strong
  lateral viscosity. Since the SIDE_DRAG keys off of the land-masks,
  it will only apply at coastlines.

    We should exercise caution if both CHANNEL_DRAG and SIDE_DRAG
  are being used (e.g., issue a warning note?), and add something to
  the documenting get_param call for SIDE_DRAG noting that as coded
  it operates at coastlines.

   The expressions (2.0 * G%mask2dT(i,j) -
    1.0) on lines 549 and 550 do
  not appear to me to have the correct staggering for a corner
  point. Please confirm that they will work as intended for
  northern, southern, eastern, and western coastlines, and at both
  convex and concave corners with all 4 possible orientations.

    I have one other coding suggestion. In MOM6 we exploit the
  case-invariance of F90 and use capital indices to indicate faces.
  (I like the mnemonic, J= j+1/2 which gets rounded down to j
  because it is integer math.) So at line 542 of your new code, we
  would prefer that you use "v_loc = v(i,J,k) + v(i+1,J,k)" instead
  of "v_loc = v(i,j,k) + v(i+1,j,k)". Once you get used to this
  convention, you can almost immediately spot when an expression is
  not appropriately staggered.  At one point Alistair even had
  written a set of CPP macros that could use this index-case
  convention to automatically shift between northeast and southwest
  grid indexing conventions and thereby automatically detect
  horizontal indexing errors.

    - Bob Hallberg

  On 05/04/2015 09:57 PM, Angus Gibson wrote:


  I'm not sure how much this overlaps with BOTTOMDRAGLAW
    either. It seemed to make more sense to do it this way for what
    I was doing.
  —
    Reply to this email directly or view
      it on GitHub.

[adcroft]

Important clarification: when Bob said "Since the SIDE_DRAG keys off of the
land-masks, it will only apply at coastlines." he meant "at MASKED
coastlines" (i.e. where the 2d MASK turns off calculations).

Dr Alistair Adcroft (Alistair.Adcroft@noaa.gov)
Princeton University Tel: (609) 987-5073
NOAA/GFDL, 201 Forrestal Road, Princeton, NJ 08540

On Tue, May 5, 2015 at 9:55 AM, Robert Hallberg notifications@github.com
wrote:

Hi Angus,

Thank you for this contribution. I can see it as being very
valuable, especially for certain idealized configurations.

This would seem to have significant overlap with the bottom drag
that is imposed on layers that would intersect with the bottom
within a grid cell via the existing CHANNEL_DRAG option. The real
ocean has no sides, so CHANNEL_DRAG can operate even when you are
not a a coastline, and because CHANNEL_DRAG is imposed implicitly,
there are no stability limits as there could be for a strong
lateral viscosity. Since the SIDE_DRAG keys off of the land-masks,
it will only apply at coastlines.

We should exercise caution if both CHANNEL_DRAG and SIDE_DRAG
are being used (e.g., issue a warning note?), and add something to
the documenting get_param call for SIDE_DRAG noting that as coded
it operates at coastlines.

The expressions (2.0 * G%mask2dT(i,j) -
1.0) on lines 549 and 550 do
not appear to me to have the correct staggering for a corner
point. Please confirm that they will work as intended for
northern, southern, eastern, and western coastlines, and at both
convex and concave corners with all 4 possible orientations.

I have one other coding suggestion. In MOM6 we exploit the
case-invariance of F90 and use capital indices to indicate faces.
(I like the mnemonic, J= j+1/2 which gets rounded down to j
because it is integer math.) So at line 542 of your new code, we
would prefer that you use "v_loc = v(i,J,k) + v(i+1,J,k)" instead
of "v_loc = v(i,j,k) + v(i+1,j,k)". Once you get used to this
convention, you can almost immediately spot when an expression is
not appropriately staggered. At one point Alistair even had
written a set of CPP macros that could use this index-case
convention to automatically shift between northeast and southwest
grid indexing conventions and thereby automatically detect
horizontal indexing errors.

• Bob Hallberg

On 05/04/2015 09:57 PM, Angus Gibson wrote:

I'm not sure how much this overlaps with BOTTOMDRAGLAW
either. It seemed to make more sense to do it this way for what
I was doing.
—
Reply to this email directly or view
it on GitHub.

—
Reply to this email directly or view it on GitHub
#169 (comment)
.

[angus-g]

Thanks for the feedback. I didn't realise CHANNEL_DRAG would do something similar to what I wanted. Indeed, in my simple test of a channel with uniform flow it replicates the expected behaviour, although the effect is much stronger for a given drag coefficient. With that in mind, I don't think it's essential that this is incorporated into MOM6, and it would have to be mutually exclusive with CHANNEL_DRAG if it was. This has been a good introduction to the C-grid staggering and things like masking for me. I'd be happy to contribute to documentation so the possible uses of the bottom drag types are more apparent.

With regard to the masking, my intention was to calculate the stress only at q-points lying on the coastline, with the need to flip the sign of the stress depending on whether this q-point is oriented to the left or right of the v-point to which it applies (for the str_yx component, the same argument applies for the str_xy component). I agree that it would probably be more appropriate to use mask2dCu and mask2dCv for this. This would still require masking on mask2dBu to ensure it was applied only at coastlines. Perhaps just dropping the abs from Alistair's suggestion:

drag_stress_u = vel_gam * u_loc * (G%mask2dCu(I,j) - G%mask2dCu(I,j+1))
drag_stress_v = vel_gam * v_loc * (G%mask2dCv(i,J) - G%mask2dCv(i+1,J))

This seems correct to me, but I'll follow up this comment with more evidence. At first glance, no stress will be applied in a concave corner (but there are no unique velocity points to which it would be applied in this case anyway). At a convex corner, uvmag will have contributions from both velocity components and stress will be applied to both components as well.

Finally, thanks for the pointer on the indexing convention. I had noticed it and tried to work it out, but it seems so obvious now! Maybe that could be added to the wiki somewhere?

[StephenGriffies]

Angus,

We are aiming to document a great deal of MOM6 via a software known as
Doxygen. This approach embeds LaTeX into the fortran code itself. The
parameterizations/vertical/MOM_KPP.F90 module has a reasonably mature set
of this documentation (look in particular at the bottom of the module).

I recommend that for your documentation work, that you also use Doxygen.
Furthermore, Alistair noted that there is a way to import .eps figures to
the Doxygen documentation. Reading through the present discussion makes it
clear that a suite of figures will go a long way towards helping to get the
boundary conditions correct and fully documented. This sort of
documentation is more robust than the Wiki, or a separate pdf file ala
MOM5, since the documentation lives directly in the Fortran file.

Best,
Steve

On Tue, May 5, 2015 at 7:34 PM, Angus Gibson notifications@github.com
wrote:

Thanks for the feedback. I didn't realise CHANNEL_DRAG would do something
similar to what I wanted. Indeed, in my simple test of a channel with
uniform flow it replicates the expected behaviour, although the effect is
much stronger for a given drag coefficient. With that in mind, I don't
think it's essential that this is incorporated into MOM6, and it would have
to be mutually exclusive with CHANNEL_DRAG if it was. This has been a
good introduction to the C-grid staggering and things like masking for me.
I'd be happy to contribute to documentation so the possible uses of the
bottom drag types are more apparent.

With regard to the masking, my intention was to calculate the stress only
at q-points lying on the coastline, with the need to flip the sign of the
stress depending on whether this q-point is oriented to the left or right
of the v-point to which it applies (for the str_yx component, the same
argument applies for the str_xy component). I agree that it would
probably be more appropriate to use mask2dCu and mask2dCv for this. This
would still require masking on mask2dBu to ensure it was applied only at
coastlines. Perhaps just dropping the abs from Alistair's suggestion:

drag_stress_u = vel_gam * u_loc * (G%mask2dCu(I,j) - G%mask2dCu(I,j+1))
drag_stress_v = vel_gam * v_loc * (G%mask2dCv(i,J) - G%mask2dCv(i+1,J))

This seems correct to me, but I'll follow up this comment with more
evidence. At first glance, no stress will be applied in a concave corner
(but there are no unique velocity points to which it would be applied in
this case anyway). At a convex corner, uvmag will have contributions from
both velocity components and stress will be applied to both components as
well.

Finally, thanks for the pointer on the indexing convention. I had noticed
it and tried to work it out, but it seems so obvious now! Maybe that could
be added to the wiki somewhere?

—
Reply to this email directly or view it on GitHub
#169 (comment)
.

Dr. Stephen M. Griffies
NOAA Geophysical Fluid Dynamics Lab
201 Forrestal Road
Princeton, NJ 08542
USA

[Hallberg-NOAA]

Angus, 

    The form for drag_stress_[uv] you suggest below makes sense to
  me and should satisfy all of our symmetry requirements.

    The reason why channel drag gives a stronger effect than your
  drag_stress than is that it makes a different set of assumptions
  about the scales over which the velocity profiles drop off near
  the edges (a subgridscale solution based on a Laplacian
  Smagorinsky viscosity) in place of the grid spacing.  I have a
  half-written document about this that I never found the time to
  complete, but if you are interested in pursuing this question
  further I could sent it to you.  Write to me at
  Robert.Hallberg@noaa.gov outside of the gitHub comment loop if you
  are interested.

    - Bob Hallberg


  On 05/05/2015 07:34 PM, Angus Gibson wrote:


  Thanks for the feedback. I didn't realise CHANNEL_DRAG
    would do something similar to what I wanted. Indeed, in my
    simple test of a channel with uniform flow it replicates the
    expected behaviour, although the effect is much stronger for a
    given drag coefficient. With that in mind, I don't think it's
    essential that this is incorporated into MOM6, and it would have
    to be mutually exclusive with CHANNEL_DRAG if it
    was. This has been a good introduction to the C-grid staggering
    and things like masking for me. I'd be happy to contribute to
    documentation so the possible uses of the bottom drag types are
    more apparent.
  With regard to the masking, my intention was to calculate the
    stress only at q-points lying on the coastline, with the need to
    flip the sign of the stress depending on whether this q-point is
    oriented to the left or right of the v-point to which it applies
    (for the str_yx component, the same argument
    applies for the str_xy component). I agree that it
    would probably be more appropriate to use mask2dCu
    and mask2dCv for this. This would still require
    masking on mask2dBu to ensure it was applied only
    at coastlines. Perhaps just dropping the abs from
    Alistair's suggestion:

    drag_stress_u = vel_gam * u_loc * (G%mask2dCu(I,j) - G%mask2dCu(I,j+1))

drag_stress_v = vel_gam * v_loc * (G%mask2dCv(i,J) - G%mask2dCv(i+1,J))

  This seems correct to me, but I'll follow up this comment with
    more evidence. At first glance, no stress will be applied in a
    concave corner (but there are no unique velocity points to which
    it would be applied in this case anyway). At a convex corner, uvmag
    will have contributions from both velocity components and stress
    will be applied to both components as well.
  Finally, thanks for the pointer on the indexing convention. I
    had noticed it and tried to work it out, but it seems so obvious
    now! Maybe that could be added to the wiki somewhere?
  —
    Reply to this email directly or view
      it on GitHub.

[adcroft]

Angus, Based on the email exchange, I'll close this PR and wait to see what you work out with Bob. -A.