fixed boundary issues in rhs_ph

[matzegoebel]

TYPE: bug fix

KEYWORDS: geopotential, index, horizontal advection, open, specified

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES:

When using advection order < 5 and open or specified boundary conditions, the upper limits in the horizontal advection loops of geopotential are one too small.

Let's look at the y-advection with advection order = 2, for instance and assume jte=jde. For open and specified BC the
lower limit is increased by 1 from jds to jds+1, but the upper limit is decreased by 2 from jde-1 to jtf=jde-1-2=jde-3. It should be jtf=jde-2, however. Starting at line 2092 the jde-1 point is set, but the jde-2 point is not handled anywhere.

The issue is analogous for x-advection and advection order = 4 (jtf=jde-4 instead of jtf=jde-3).

For advection order >= 5, the limits are set correctly, because a different notation is used: jtf = min(jtf,jde-4)

LIST OF MODIFIED FILES:
dyn_em/module_big_step_utilities_em.F

TESTING CONDUCTED:
I did an LES test case simulation, with periodic y-boundaries and open x-boundaries, mean x-wind of 10 m/s, h_sca_adv_order=3, and north-south oriented cosine ridge in the center of the domain.
Here is the perturbation geopotential after 1 minute integration at k=2:

and here the difference plot:

You can clearly see the higher geopotential in the new formulation at i=ide-3. This is the point where the horizontal advection is not carried out by the current algorithm.

After 10 minutes it looks like this:

RELEASE NOTES: When using advection order < 5, and either open or specified boundary conditions, the upper limit of the indexing in the horizontal advection loops of geopotential was previously incorrect. For horizontal advection order >= 5, the index limits were not incorrect.

[davegill]

@matzegoebel
Matthias,
This is a very clean PR. We very much appreciate when users contribute mods that assist a large segment of the community.

1. Would you please provide a method to allow us to verify that this is solving the stated problem. For example, a before vs after plot, or some additional diagnostic prints.

2. When supporting verification materials are provided, also edit the PR message a bit. Instead of

I think it should be jtf=jde-2, however.

the new statement will be

It should be jtf=jde-2, however.

[matzegoebel]

Ok, I did an LES test case simulation, with periodic y-boundaries and open x-boundaries, mean x-wind of 10 m/s, h_sca_adv_order=3, and north-south oriented cosine ridge in the center of the domain.
Here is the perturbation geopotential after 1 minute integration at k=2:

and here the difference plot:

You can clearly see the higher geopotential in the new formulation at i=ide-3. This is the point where the horizontal advection is not carried out by the current algorithm.

After 10 minutes it looks like this:

[davegill]

@matzegoebel
Matthias,
Would you please fix this line. It is a bit confusing:

but the upper limit is decreased by 2: jtf=jtf-2=jde-1-2=jde-3. 

Are you saying that jtf is originally jde-1, but then we reduce it further? So effectively jtf=jde-3?

[matzegoebel]

Are you saying that jtf is originally jde-1, but then we reduce it further? So effectively jtf=jde-3?

Yes correct. I tried to make this a bit more clear in the commit message now.