Bug in surface roughness ZNT/Z0 inconsistency when using urban physics

[cenlinhe]

Describe the bug
This is a potential bug for surface roughness calculation/diagnostic inconsistency when activating urban physics as reported in the WRF forum: https://forum.mmm.ucar.edu/threads/why-does-the-drag-actually-decrease-and-wind-speed-increase-after-enabling-the-urban-canopy-model.25191/#post-59296

The wrfout ZNT/Z0 variables are the values from land model calculations assuming "nature" vegetation type when UCM is activated. The ZNT/Z0 values after passing to urban physics are updated to be Z0C (e.g., from urban parameter table) within urban physics modules for urban physics calculation. However, these updated urban ZNT/Z0 is not computed to be an updated grid-mean values to pass back to other wrf physics or wrfout modules, which appears to be a bug. At least, there should be a weighted average grid-level mean ZNT/Z0 as a diagnostic output variable to be used by other WRF physics. These does not affect the urban physics calculations though.

Proposed Solution
One solution would be to update the grid-mean ZNT/Z0 by weighted averaging urban and non-urban ZNT/Z0 after urban physics call so that other WRF physics and wrf output can use the correct surface roughness values.

[cenlinhe]

@tslin2 @dudhia @chenghaow Could you please take a look at this issue and confirm if this is a bug?

[tslin2]

The interpretation of BEP scheme is not directly related to Z0 because it estimates the momentum sink with a drag force that depends of the urban morphology.
For bulk (no UCM) and SLUCM, they use a roughness length that is not directly dependent on the urban morphology.

I am fine with the changes, this was made in Noah mosaic version:

WRF/phys/module_sf_noahdrv.F

Line 3953 in 1957ed7

ZNT(I,J)= EXP(FRC_URB2D(I,J)*ALOG(ZNT_URB)+(1-FRC_URB2D(I,J))* ALOG(ZNT(I,J))) ! ADD BY DAN

Also some previous related discussions:
DanLi-BU#5 (comment)
#2093