TYPE: bug fix

KEYWORDS: rho, sfs_opt, m_opt, cal_titau

SOURCE: Robert Arthur and Jeff Mirocha (LLNL)

DESCRIPTION OF CHANGES:
Problem:
In the cal_titau subroutines in module_diffusion_em, rho was not always interpolated to the correct location on the 
staggered grid. This issue affected the subgrid stress terms when sfs_opt .gt. 1. Additionally, when sfs_opt .eq. 0 and 
m_opt .eq. 1, the actual subgrid stresses applied in the code were correct, but the output was wrong. This is a 
follow-on from issue #1214 by @matzegoebel, which addressed the surface stress output from subroutine 
vertical_diffusion_2 when sfs_opt .gt. 0 or m_opt .eq. 1, but did not address the cal_titau subroutines.

Solution:
A rhoavg variable was added to cal_titau_12_21, cal_titau_13_31, and cal_titau_23_32 and used instead of the 
cell-centered rho variable. Additionally, because cal_titau_12_21 now uses "corner" ghost points for rho, a new 
halo/bc communication was added to solve_em.F to ensure that rho is specified correctly at those points. Correct 
specification of the corner points also requires the updates to subroutine set_physical_bc3d in #1314, included in 
the v4.2.2 release. Without the latter two updates for corner points, bit-for-bit errors occurred in regression tests for 
idealized, doubly-periodic domains.

LIST OF MODIFIED FILES: 
dyn_em/module_diffusion_em.F
dyn_em/solve_em.F
Registry/Registry.EM_COMMON

TESTS CONDUCTED:
1. Idealized, doubly-periodic tests with various processor breakdowns show bit-for-bit identical results.
2. Jenkins is all pass.

RELEASE NOTE: Bug fix to use correct rho value in cal_titau subroutines. In the cal_titau subroutines in module_diffusion_em, rho was not always interpolated to the correct location on the staggered grid. This issue affected the subgrid stress terms when sfs_opt .gt. 1. Additionally, when sfs_opt .eq. 0 and m_opt .eq. 1, the actual subgrid stresses applied in the code were correct, but the output was wrong. A rhoavg variable was added to cal_titau_12_21, cal_titau_13_31, and cal_titau_23_32 and used instead of the cell-centered rho variable.

…del#1314)

TYPE: bug fix

KEYWORDS: boundary conditions, symmetric, corners

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES:
Problem:
In idealized simulations, the routine `set_physical_bc3d` sets the boundary zone points for periodic, symmetric, and open boundary conditions.
At first, the routine does not set the corner points (e.g., (ids-1,jds-1)): For the y-boundaries, it iterates i only from `MAX(ids,its-1)` to  `MIN(ite+1,ide+istag)` and analogously for the x-boundaries (note that for open BC it actually also tries to set the corner points (see line 896), but not effectively, because the RHS values in lines 903-905 are still 0 there, as the y-BC has not been applied yet).

For doubly-periodic BC, the corner points are set correctly afterward, starting at line 1100. In any other case, however, it seems that the corner points are not assigned and appear as zeros. 

These corner points are later used, e.g. in the SGS UV-flux (titau_21_12) in the diffusion module:
https://github.com/wrf-model/WRF/blob/f311cd5e136631ebf3ebaa02b4b7be3816ed171f/dyn_em/module_diffusion_em.F#L5511-L5517
Thus, because `rho(i_start-1,k,j_start-1)=0` and `xkx(i_start-1,k,j_start-1)=0`, we have `xkxavg(i_start, k, j_start)` much smaller than at the other locations.

This is usually not a problem, since the deformation defor12 and thus also the flux titau_21_12 should be 0 at these corner points when using symmetric boundary conditions on at least one of these boundaries.
I write "should" because that's a theoretical argument. In practice, only the deformation at the ids or jds boundaries is zero. Due to an inconsistency in the symmetric BC for the deformation, the deformation at the ide or jde boundaries is not 0. This problem is described in the following.

The deformation (variable = 'd') is treated like an unstaggered variable: https://github.com/wrf-model/WRF/blob/f311cd5e136631ebf3ebaa02b4b7be3816ed171f/share/module_bc.F#L1003-L1011
Thus, the deformation which was initially zero at the boundary is set to a non-zero value by `dat(i,k,jde) = dat(i,k,jde-1)`.
All staggered variables except the normal velocity (so also the deformation) should be treated equally in the symmetric boundaries since the location of the symmetric boundary is on the staggered grid. Therefore I adjusted the relevant if statement in the second commit.

If open boundary conditions are used on any of the boundaries, only the fluxes on the interior are calculated, so the corner points also don't matter. However, to avoid any current or future side effects due to the missing corner points, corner points are set anyway as it is done in `set_physical_bc2d`.

Solution:

1. commit:
The solution to set the corner points is the same as is already applied in the 2d routine `set_physical_bc2d`.

2. commit:
To treat all variables staggered in a certain direction equally, the variables `jstag` and `istag` are used in the if statement instead of checking the value of `variable`. This was done only for the symmetric BC. For the open BC, it might be necessary, too.

3. commit:
The corner point assignment for doubly-periodic BC is now redundant and was removed.

ISSUE:
Fixes wrf-model#1311

LIST OF MODIFIED FILES: share/module_bc.F

TESTS CONDUCTED:
With the first commit, I tested, whether the corner points for doubly-periodic BC are set in the same way as in the original code. This test allowed removing this redundant part of the code in commit 3.

To compare the different model versions, I ran an idealized simulation with km_opt=2, diff_opt=2, tke_heat_flux=0.2, U = 5 m/s, periodic in x, symmetric in y.
The following figure shows the titau_21_12 flux at k=kds after 20 minutes for three model versions:
left: original, middle: fixed corners (commit 1), right: fixed corners and fixed symmetric BC for the deformation (commit 1 and 2)
![fuy_sgs](https://user-images.githubusercontent.com/17001470/98154601-505be880-1ed5-11eb-9f0a-6eef1175e730.png)

All versions have titau_21_12=0 at j=jds. In the left and middle panel, the fluxes at jde and jde-1 are almost identical. However, the original version has fluxes with lower magnitudes at the upper left and right corners (as the corners of xkmh and rho are not set as explained above). In the right panel, the flux at j=jde is now also zero as it should be.

The following figure shows a time series of U at the upper left corner:
![u_timeseries](https://user-images.githubusercontent.com/17001470/98154606-52be4280-1ed5-11eb-93d3-563062274e5c.png)

We can see that fixing the symmetric BC to have zero fluxes at j=jde does have a significant effect.
Fixing the corners (commit 1) only has a minor effect, which is maximized at the corner points (j=jde, i=ids, and ide) as can be seen here:
![u_diff](https://user-images.githubusercontent.com/17001470/98154624-55b93300-1ed5-11eb-906d-b8aaef7f1596.png)

RELEASE NOTE: 
Missing corner points are set when setting boundary conditions of 3D variables and all staggered variables are treated equally in symmetric BC. Thereby the periodic BC become cleaner and more consistent with set_physical_bc2d and the symmetric BC are fixed as horizontal deformation is treated now as a staggered variable which leads to correct BC for SGS momentum fluxes at the northern boundaries.

…1396)

TYPE: bug fix

KEYWORDS: rho, sfs_opt, m_opt, cal_titau

SOURCE: Robert Arthur and Jeff Mirocha (LLNL)

DESCRIPTION OF CHANGES:
Problem:
In the cal_titau subroutines in module_diffusion_em, rho was not always interpolated to the correct location on the 
staggered grid. This issue affected the subgrid stress terms when sfs_opt .gt. 1. Additionally, when sfs_opt .eq. 0 and 
m_opt .eq. 1, the actual subgrid stresses applied in the code were correct, but the output was wrong. This is a 
follow-on from issue wrf-model#1214 by @matzegoebel, which addressed the surface stress output from subroutine 
vertical_diffusion_2 when sfs_opt .gt. 0 or m_opt .eq. 1, but did not address the cal_titau subroutines.

Solution:
A rhoavg variable was added to cal_titau_12_21, cal_titau_13_31, and cal_titau_23_32 and used instead of the 
cell-centered rho variable. Additionally, because cal_titau_12_21 now uses "corner" ghost points for rho, a new 
halo/bc communication was added to solve_em.F to ensure that rho is specified correctly at those points. Correct 
specification of the corner points also requires the updates to subroutine set_physical_bc3d in wrf-model#1314, included in 
the v4.2.2 release. Without the latter two updates for corner points, bit-for-bit errors occurred in regression tests for 
idealized, doubly-periodic domains.

LIST OF MODIFIED FILES: 
dyn_em/module_diffusion_em.F
dyn_em/solve_em.F
Registry/Registry.EM_COMMON

TESTS CONDUCTED:
1. Idealized, doubly-periodic tests with various processor breakdowns show bit-for-bit identical results.
2. Jenkins is all pass.

RELEASE NOTE: Bug fix to use correct rho value in cal_titau subroutines. In the cal_titau subroutines in module_diffusion_em, rho was not always interpolated to the correct location on the staggered grid. This issue affected the subgrid stress terms when sfs_opt .gt. 1. Additionally, when sfs_opt .eq. 0 and m_opt .eq. 1, the actual subgrid stresses applied in the code were correct, but the output was wrong. A rhoavg variable was added to cal_titau_12_21, cal_titau_13_31, and cal_titau_23_32 and used instead of the cell-centered rho variable.