Bug fixes for stability function and energy balance in single canopy …

…urban model (#2038)

TYPE: Bug fix

KEYWORDS: Stability function, 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).

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).

LIST OF MODIFIED FILES: 
M     phys/module_sf_urban.F

TESTS CONDUCTED: 
1. No tests are conducted. 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. 
2. The Jenkins tests are all passing.

RELEASE NOTE: Fixes for the calculations for PSIM in the mos subroutine and the temperature used in surface energy balance for the green roof. The later follows the OSU1DPBL document located at:  
[https://ftp.emc.ncep.noaa.gov/mmb/gcp/ldas/nceplsm/OSU1DPBL/OSU1DPBL-userguide.pdf]

phys/module_sf_urban.F

@@ -1184,7 +1184,7 @@ SUBROUTINE urban(LSOLAR,                                           & ! L
           
        CALL TDFCND (DF1,SMR(KZ), QUARTZ, SMCMAX )
        DF1 = DF1 * EXP(-2.0 * SHDFAC)  
-       RGR = EPSV*(RX-SIG*(TGRP**4.)/60.)
+       RGR = EPSV*(RX-SIG*(TA**4.)/60.)
        RGRR= (SGR+RGR) * 697.7 * 60.
        RCH = RHOO*CPP*CHGR
        RR1 = EPSV*(TA**4) * 6.48E-8 / (PS* CHGR) + 1.0
@@ -1727,7 +1727,7 @@ SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
         PSIM=ALOG((Z+Z0)/Z0) &
             -ALOG((X+1.)**2.*(X**2.+1.)) &
             +2.*ATAN(X) &
-            +ALOG((X+1.)**2.*(X0**2.+1.)) &
+            +ALOG((X0+1.)**2.*(X0**2.+1.)) &
             -2.*ATAN(X0)
         FAIH=1./SQRT(1.-16.*XXX)
         PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &