Stability function correction for momentum in mos subroutine of SLUCM.

[joshi994]

The first line should be a single-line "purpose" for this change

TYPE: Bug fix

KEYWORDS: Stability function, stability parameter, mos subroutine, Monin-Obukhov similarity theory, Urban canopy
model, surface energy balance, Green roof temperature.

SOURCE: Parag Joshi (Brookhaven National Laboratory), Katia Lamer (Brookhaven National Laboratory)

DESCRIPTION OF CHANGES:
Problem: The calculations for PSIM in the mos subroutine of module-sf_urban.F has been corrected. The earlier equations has a typo where X (= (1-16\zeta)^1/4) is mentioned instead of X0. In addition, following the formulation of surface energy balance of OSU1DPBL scheme by Ek and Mahrt 1991, it appears that when the following definition is used $\sigma * \theta_s^4 = \sigma T_0^4 (1 + 4((\theta_s - T_0)/T_0)))$ is used, the temperature in the calculations of net flux at the surface should use T_0 (temperature at first atmospheric level). However, in the module TGRP (temperature of the green roof surface at previous time step is used)is used instead of T_0 (please refer to the attached screenshots of the formulation).

Also I am not sure why the term (\theta_0 - T_0) is neglected in the definition of YY (Line/command 1145 WRFv-4.5.2) whereas the equation for the $\theta_s$ in attached screenshot shows the term. Due to insignficant impact?

Solution:
Following the formulation, the temperature variable has been updated according. Moreover, the mos subroutine is also corrected by replacing X by X0 where X0 is =X*z0/(z+z0) (as per the WRF/phys/module_sf_urban.F).

ISSUE: For use when this PR closes an issue.
Fixes #123

LIST OF MODIFIED FILES: module_sf_urban.F

TESTS CONDUCTED:

1. A test case could be run to compare the impact of changes in the TGR and stability parameter value as well as CD, US, and ALPHA.

For more details, refer to the pages # 47-50 of the document link in the release note.

RELEASE NOTE: The surface energy balance for the green roof is modified following the OSU1DPBL scheme. The document can be located at:
[https://ftp.emc.ncep.noaa.gov/mmb/gcp/ldas/nceplsm/OSU1DPBL/OSU1DPBL-userguide.pdf]

]

[weiwangncar]

The regression test results:

Test Type              | Expected  | Received |  Failed
= = = = = = = = = = = = = = = = = = = = = = = =  = = = =
Number of Tests        : 23           24
Number of Builds       : 60           57
Number of Simulations  : 158           150        0
Number of Comparisons  : 95           86        0

Failed Simulations are: 
None
Which comparisons are not bit-for-bit: 
None

[cenlinhe]

I think here the original code is correct. We need to use "TGRP" for the green roof surface temperature to compute the outgoing longwave radiation for the greenroof surface not the air. But this RGR does lead to the inconsistence in the "YY" calculation using "TA" a few lines below. I am not quite sure about this yet. this is the reference for this greenroof development: https://link.springer.com/article/10.1007/s10546-014-9991-6
maybe you could contact the original developer: Prof. Jiachuan Yang (HKUST): cejcyang@ust.hk

[joshi994]

@cenlinhe I have tried explaining my part in the attachment here. Please let me know if you need further clarification. Also, please correct me if I am terribly wrong.
TA_TGRP_Clarification.pdf

[joshi994]

@cenlinhe I believe that for consistency in estimating YY, TA needs to be replaced by TGRP while keeping TGRP in the calculations of RGR (as done for the roof case where TRP is used instead), since an iterative approach is implemented. This may provide better results.

[cenlinhe]

@cenlinhe I believe that for consistency in estimating YY, TA needs to be replaced by TGRP while keeping TGRP in the calculations of RGR (as done for the roof case where TRP is used instead), since an iterative approach is implemented. This may provide better results.

I agree that TA should be replaced by TGRP in YY calculation

[cenlinhe]

@cenlinhe I have tried explaining my part in the attachment here. Please let me know if you need further clarification. Also, please correct me if I am terribly wrong. TA_TGRP_Clarification.pdf

I agree that current code has some inconsistency in RGR, RR1, and YY calculations. From the greenroof surface energy balance perspective, I think TGRP should be used instead of TA in all these three places. Also, the term (theta0 - T0) is missing but I do not know the reason, likely because the assumption that they are the same.

[joshi994]

I think the iterative approach is required to achieve the greenroof surface energy flux balance convergence, which means this equation needs to be inside the iteration loop. Could you please do a test to compare the simulations with TGRP and replacing TGRP by TA (based on your proposed solution) to see how much difference it makes?

@cenlinhe Unfortunately, I do not have computing facilities to run a test simulation. I tried on my local computer but it didn't work well. I will keep you posted if anything changes.

[jyang104]

Thanks for pointing out the issue. After carefully reviewing the code and the PDF, to be in consistent with the equations, TA shall be used instead of TGRP in the calculation of RGR, RR1, and YY.

For the term (theta0 - T0), their values shall be close at the first atmospheric level, so I ignore that term during the coding process.

[joshi994]

Thanks for pointing out the issue. After carefully reviewing the code and the PDF, to be in consistent with the equations, TA shall be used instead of TGRP in the calculation of RGR, RR1, and YY.

For the term (theta0 - T0), their values shall be close at the first atmospheric level, so I ignore that term during the coding process.

@jyang104 Thank you so much for your response and for confirming the inconsistency. Now this makes me wonder that when a similar energy balance is applied to a common roof, wall, and ground shouldn't the $\sigma*\theta_s^4$, where $\theta_s$ can represent the surface temperature of roof, wall, and ground, must use a similar approximation $\sigma*\theta_s^4 = \sigma*\theta_A^4 (1 + 4(\theta_s - T_A)/T_A$? I believe that the nonlinear nature of the resulting balance equation in $\theta_s$ is the sole reason for such an approximation. Hence directly using TRP, TWP, or TGP (even if it is from the previous time step) in the calculation of the longwave radiation part of the energy balance causes further inconsistencies. Please correct me if I am wrong. I am extremely new to WRF modeling and trying to understand the physics.

[jyang104]

@joshi994 You are correct. Part of the reason for the inconsistency is that the calculation of outgoing longwave in the original sf_urban (i.e, for wall and ground and roof) does not involve evapotranspiration that the equation of \sigma*\theta_s^4 is used. When I implemented the green roof, I followed the sf_noahlsm in general that the surface energy balance used the approximation formula. In this case, it will be easier to just correct the green roof part of the code, rather than fixing the ones for TRP, TWP, or TGP.

[joshi994]

@jyang104 @cenlinhe It make sense corecting the green roof part only. However, it would be nice to see the temperature of the green roof calculated using YY and ZZ1, the linear approximation $\sigma*\theta_s^4 = \sigma*\theta_A^4 (1 + 4(\theta_s - T_A)/T_A$ and using the Newton-Raphson approach right below it. I am leaving it over to the experts to decide and close the request. Thanks for the input @jyang104 and @cenlinhe for providing me the opportunity to speak with @jyang104.

[cenlinhe]

This typo is indeed a bug.

[cenlinhe]

I am also bringing @barlage here just in case he has any insights on this problem.