Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Forcing and output in SLUCM: actual temperature or potential temperature #2093

Open
xuelingbo opened this issue Aug 7, 2024 · 22 comments
Open

Comments

@xuelingbo
Copy link

In the forcing of subroutine urban, temperature TA is from T3D which is actual temperature:
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3393
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3499
After passing T3D to subroutine urban,
Firstly, TS is calculated according to TA and sensible heat flux :

TS = TA + FLXTH/CHS ! surface potential temp (flux temp)

Then, TH2 is calculated according to TA and TS:
TH2 = TS + (TA-TS) *(CHS/CHS2)

According to these codes, TA, TS, and TH2 should all be actual temperatures.
But in the diagnostic 2m air temperature T2, TH2 is treated as 2m potential temperature:
(TH2_URB2D(i,j)/((1.E5/PSFC(i,j))**RCP))*FRC_URB2D(I,J)

This seems somewhat inconsistent
Additionally, another inconsistent part is the comments for TA and TH2, both of which indicate they are potential temperatures. :
! TA [K] : Potential temperature at 1st wrf level (absolute temp)

! TH2 [K] : Diagnostic potential temperature at 2 m

In subroutine urban, TA is also used in the following codes, which I think actual temperature should be used:

TAV=TA*(1.+0.61*QA)

TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)

HR=RHO*CP*CHR*UA*(TRP-TA)*100.

TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP

XXX = 0.4*9.81*Z*TST/TA/UST/UST

@xuelingbo
Copy link
Author

We did some simulations using HRLDAS, and compare the diagnostic T2 of methods used by WRF and method that TH2_URB2D is considered as air temperature: xuelingbo/LSP-DS#5.

@cenlinhe
Copy link
Contributor

Looks like based on @xuelingbo 's tests, for low-altitude cities, the difference is very small, but for high-altitude urban regions, the difference is huge (up to ~4 degC). My understanding is that TH2 inside urban module code is the actual 2-m air temperature instead of potential temperature.
@joshi994 @chenghaow @tslin2

@tslin2
Copy link

tslin2 commented Aug 27, 2024

I guess some of the equations in SLUCM are required to use air temperature, so they probably made the assumption that air temperature is similar to the potential temperature in SLUCM. I will do a test case in a city for this modification.

@tslin2
Copy link

tslin2 commented Aug 27, 2024

another question I had is

! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K]

how about this line of the code, it was commented out, what is the difference between (CHS/CHS2) and (PSIT2/PSIT)

@joshi994
Copy link
Contributor

@cenlinhe @tslin2 @xuelingbo I agree with @tslin2 on the use of potential temperature. However, the diagnostic potential temperature (TH2) basically denoted the compression/expansion effects due to pressure change and within 2m height ground from the ground, pressure change is assumed to be having no effect. In addition, diagnostic potential temperature is a reference temperature.

@tslin2 There is a slight difference between PSIT and CHS variables if you go by their definition. PSIT I believe indicate the similarity function for the temperature whereas CHS denotes the exchange coefficient derived using those similarity functions. The exchange coefficient equation also include the urban roughness effects by including friction velocity in its definition hence, CHS appears to be a better option where both thermal and mechanical roughness impacts can be included while calculating the potential temperature. Let me know your views.

@xuelingbo
Copy link
Author

@tslin2 @joshi994 I understand they probably made the assumption that air temperature is similar to the potential temperatur as 2m height is near the ground. But TA is actual temperature at 1st atmospheric layer (~30m), I am concerned if it's correct to use this assumption.
I guess the assumption may have been used at the time to reduce the amount of calculations. So is it now possible to pass both the actual temperature(t_phy) and potential temperature(theta_urban) to the SLUCM, choose the appropriate one in the calculation and end up outputting a 2m actual temperature t2 that is consistent with t2mv, t2mb.

@tslin2
Copy link

tslin2 commented Aug 29, 2024

This is not directly relevant,

WRF/phys/module_sf_urban.F

Lines 428 to 429 in 0a11865

REAL, INTENT(OUT) :: U10 ! u at 10m [m/s]
REAL, INTENT(OUT) :: V10 ! u at 10m [m/s]

U10, V10 should be INTENT(INOUT) due to the use for Calculating Wind Direction and Assign Appropriae lf_urb

WDR = (180.0/PI)*ATAN2(U10,V10)

or use U1, V1 to calculate WDR

@tslin2
Copy link

tslin2 commented Aug 29, 2024

@xuelingbo @joshi994 thanks for your views.

@xuelingbo agree on your view to have a consistent method for Noah-MP and SLUCM, this can also be related to the potential comparison with the BEP method for T2 using potential temperature.

@joshi994 the comment DanLi-BU#5 (comment) by @DanLi-BU somehow justify the use of CHS from the surface layer scheme that considers momentum and thermal roughness based on Noah-SLUCM.

In that sense, should Q2 be calculated from CHS2 and CHS?,

Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]

due to

WRF/phys/module_sf_urban.F

Lines 1628 to 1629 in 0a11865

TS = TA + FLXTH/CHS ! surface potential temp (flux temp)
QS = QA + FLXHUM/CHS ! surface humidity

@DanLi-BU
Copy link
Contributor

I did look into this issue a few months ago while I was writing a paper. Although I was not using Noah-MP (I was using Noah), I did spend sometime thinking about the calculation of T2.

Noah SLUCM

I'll start with Noah-SLUCM. In Noah-SLUCM, the TH2 computed by module_sf_urban.F is simply not used. One can verify this from the following code:

IF(SF_URBAN_PHYSICS.eq.1) THEN
DO j=j_start(ij),j_end(ij) !urban
DO i=i_start(ij),i_end(ij) !urban
IF(IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 )THEN
U10(I,J) = U10_URB2D(I,J) !urban
V10(I,J) = V10_URB2D(I,J) !urban
PSIM(I,J) = PSIM_URB2D(I,J) !urban
PSIH(I,J) = PSIH_URB2D(I,J) !urban
GZ1OZ0(I,J) = GZ1OZ0_URB2D(I,J) !urban
!m AKHS(I,J) = AKHS_URB2D(I,J) !urban
AKHS(I,J) = CHS(I,J) !urban
AKMS(I,J) = AKMS_URB2D(I,J) !urban
END IF !urban
ENDDO !urban
ENDDO !urban
ENDIF

In other words, the T2 of urban grid cells is still computed by module_sf_sfcdiags.F. In this case, whether TH2 is an air temperature or potential temperature doesn't make a difference.

Noah-MP SLUCM

Noah-MP does things differently. As @xuelingbo showed, T2 is computed as some sort of area-average:

IF(SF_URBAN_PHYSICS.eq.1) THEN
DO j=j_start(ij),j_end(ij) !urban
DO i=i_start(ij),i_end(ij) !urban
IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 )THEN
Q2(I,J) = (FVEGXY(I,J)*Q2MVXY(I,J) + (1.-FVEGXY(I,J))*Q2MBXY(I,J))*(1.-FRC_URB2D(I,J)) + &
Q2_URB2D(I,J)*FRC_URB2D(I,J)
T2(I,J) = (FVEGXY(I,J)*T2MVXY(I,J) + (1.-FVEGXY(I,J))*T2MBXY(I,J))*(1.-FRC_URB2D(I,J)) + &
(TH2_URB2D(i,j)/((1.E5/PSFC(i,j))**RCP))*FRC_URB2D(I,J)
TH2(I,J) = T2(i,j)*(1.E5/PSFC(i,j))**RCP
U10(I,J) = U10_URB2D(I,J) !urban
V10(I,J) = V10_URB2D(I,J) !urban
PSIM(I,J) = PSIM_URB2D(I,J) !urban
PSIH(I,J) = PSIH_URB2D(I,J) !urban
GZ1OZ0(I,J) = GZ1OZ0_URB2D(I,J) !urban
AKHS(I,J) = CHS(I,J) !urban
AKMS(I,J) = AKMS_URB2D(I,J) !urban
END IF !urban
ENDDO !urban
ENDDO !urban
ENDIF

In this case, whether TH2 is an air temperature or potential temperature does make a difference. How TH2 is computed makes perhaps an even more important difference.

Q1: is TH2 an air temperature or potential temperature

Long story short, it is an air temperature. So the (1.E5/PSFC(i,j))**RCP) correction in the above code (line 3392) is incorrect.

More broadly, the use of air temeprature throughout module_sf_urban.F is concerning. For example, sensible heat flux is supposedly defined based on potential temperature, not air temeprature. Air and potential temperatures can be quite different at the 1st atmospheric model level, so @xuelingbo 's concern is well justified. In my opinion, we should use potential temperature for all the sensible heat flux calculations in module_sf_urban.F. Air temperature is still needed for diagnostic purposes and for computing saturated water vapor pressure. This will require some major changes in the module_sf_urban.F code. I'm not sure whether Noah/Noah-MP define sensible heat flux based on potential temperature or not.

Q2: how should T2 be diagnosed for a grid cell conceptualized as part of grass (handled by Noah or Noah-MP) and part of impervious land (handled by SLUCM)

This is a bigger issue and I thought about it a lot. The challenge is associated with conceptualizing the atmospheric surface layer between the land surface and the 1st atmospheric level. One can think of two scenarios/approaches: 1) turbulence in the atmospheric surface layer is so mixed that all patches see one turbulent transfer coefficient and one flux, and 2) turbulence is isolated that different patches see different turbulent transfer coefficients and different fluxes. Both are consistent with the assumption that the 1st atmospheric model level receives the total flux, but they have vastly different implications for the temperature profiles in the atmospheric surface layer and thus TSK/T2.

As far as I know, Noah-SLUCM uses the first approach. This is why TS and QS of the impervious patch are calculated using CHS which is a turbulent transfer coefficient computed outside of module_sf_urban.Fand used by the grass patch. More details can be found in DanLi-BU#5 (comment). A side note: in this case, the CHS should be computed with roughness lengths representing the entire landscape (some combination of grass and impervious patch roughness lengths). However, in the current code it's computed using the roughness lengths of grass. This is why some people reported that urban roughness lengths are smaller than those of cropland (https://forum.mmm.ucar.edu/threads/why-roughness-znt-of-urban-is-lower-than-croplands-in-ucm.11006/). And this is technically a bug.

This also explains why T2/Q2 are computed by module_sf_sfcdiags.F using the average flux, and why TH2/Q2 computed in module_sf_urban.F are not used (since they are computed using impervious patch only fluxes, Q2 is even computed using impervious patch only turbulent transfer coefficient, the PSIT and PSIT2).

As far as I can see, Noah-MP-SLUCM diagnoses TSK in a similar way as Noah-SLUCM. In this case, Noah-MP-SLUCM has implicitly chosen the 1st approach. That being said, the area average of T2 is inconsistent with the spirit of the 1st approach.

If we were to implement the 2nd approach say in Noah-MP, we would have to diagnose TS and T2 using the respective fluxes and transfer coefficients. To do so, we need to replace CHS and CHS2 in the calculation of TS and T2 in module_sf_urban.F by the impervious land only transfer coefficients. In short, we diagnose TSK/T2 for each patch separately (using respective fluxes and transfer coefficients) and then area average them.

@tslin2
Copy link

tslin2 commented Aug 30, 2024

Thanks @DanLi-BU for comments, agree on the sensible heat of using potential temperature for SLCUM

@xuelingbo
Copy link
Author

xuelingbo commented Sep 2, 2024

@DanLi-BU Thank you for your comments, I fully understand your discussion about how should we diagnostic T2 in Noah-MP-SLUCM. Just for reminder, in HRLDAS, CHS and CHS2 is calculated inside SLUCM because it's offline, so I think it's proper to diagnostic T2 in HRLDAS use three-tiles-average. @tslin2
https://github.com/NCAR/hrldas/blob/e6853efc84d884836122e42e02a9ea2e748ec75b/urban/wrf/module_sf_urban.F#L1537
https://github.com/NCAR/hrldas/blob/e6853efc84d884836122e42e02a9ea2e748ec75b/urban/wrf/module_sf_urban.F#L1561

@chenghaow
Copy link

another question I had is

! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K]

how about this line of the code, it was commented out, what is the difference between (CHS/CHS2) and (PSIT2/PSIT)

A related issue is Q2. Regardless of whether CHS or PSIT should be used for Q2, the following statement:

Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]

should be written as

Q2 = QS + (QA-QS)*(PSIT/PSIT2)

so that it is consistent with the 2-m air temperature calculation:

TH2 = TS + (TA-TS) *(CHS/CHS2)

@chenghaow
Copy link

As @xuelingbo and @cenlinhe noted, the difference in H can become much bigger with increasing elevation. But to my knowledge, many UCMs still use air temperature to calculate sensible heat flux.

I just did a quick search in Noah-MP. Interestingly, for sensible heat flux from the surface, it also uses surface temperature and air temperature (rather than potential temperature). See:

HeatSensibleBareGrd = ShCoeff * (TemperatureGrdBare - TemperatureAirRefHeight)
Source: https://github.com/NCAR/noahmp/blob/fa2d2098bac1e437a21ccc15494b152996d5bc84/src/SurfaceEnergyFluxBareGroundMod.F90#L153

TemperatureAir2mBare = TemperatureGrdBare - HeatSensibleBareGrd / & (DensityAirRefHeight*ConstHeatCapacAir) * 1.0 / ExchCoeffSh2mBare
Source: https://github.com/NCAR/noahmp/blob/fa2d2098bac1e437a21ccc15494b152996d5bc84/src/SurfaceEnergyFluxBareGroundMod.F90#L212-L213

TVIR = (1.0 + 0.61*SpecHumidityRefHeight) * TemperatureAirRefHeight
Source: https://github.com/NCAR/noahmp/blob/fa2d2098bac1e437a21ccc15494b152996d5bc84/src/ResistanceBareGroundMostMod.F90#L94

@cenlinhe Thoughts?

@DanLi-BU
Copy link
Contributor

DanLi-BU commented Sep 4, 2024

@xuelingbo I took a quick look at the modifications in HRLDAS. You are right that this seems consistent with the second approach I described (diagnosing T2/TS for each patch separately and then aggregate them). However, the counter-argument is if so, then why don't we diagnose T2/TS for roof and canyon separately? Remember, in the current SLUCM, the fluxes from roof and canyon are calculated separately. In other words, the roof and canyon are treated as different patches.

Another way of looking at this issue is to ask the question why is the calculation of CHS/CHS2 only dependent on Z0HC, or why is Z0HR playing no role in the calculation of CHS/CHS2? We are trying to diagnose a T2/TS for the entire impervious land which composes of a roof and a canyon, I can't think of a reason why Z0HC should be used to represent the heat transfer capability of the entire impervious land. Note here we are talking about heat transfer, not momentum transfer. Maybe the role of roof can be neglected for momentum transfer, but I'm not sure about heat transfer.

I don't have good answers to these questions though :)

@DanLi-BU
Copy link
Contributor

DanLi-BU commented Sep 4, 2024

another question I had is

! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K]

how about this line of the code, it was commented out, what is the difference between (CHS/CHS2) and (PSIT2/PSIT)

A related issue is Q2. Regardless of whether CHS or PSIT should be used for Q2, the following statement:

Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]

should be written as

Q2 = QS + (QA-QS)*(PSIT/PSIT2)

so that it is consistent with the 2-m air temperature calculation:

TH2 = TS + (TA-TS) *(CHS/CHS2)

@chenghaow I guess CHS/CHS2 are conductances, PSIT/PSIT2 are resistances.

@joshi994
Copy link
Contributor

joshi994 commented Sep 4, 2024

In the forcing of subroutine urban, temperature TA is from T3D which is actual temperature: https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3393 https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3499 After passing T3D to subroutine urban, Firstly, TS is calculated according to TA and sensible heat flux :

TS = TA + FLXTH/CHS ! surface potential temp (flux temp)

Then, TH2 is calculated according to TA and TS:

TH2 = TS + (TA-TS) *(CHS/CHS2)

According to these codes, TA, TS, and TH2 should all be actual temperatures.
But in the diagnostic 2m air temperature T2, TH2 is treated as 2m potential temperature:

(TH2_URB2D(i,j)/((1.E5/PSFC(i,j))**RCP))*FRC_URB2D(I,J)

This seems somewhat inconsistent
Additionally, another inconsistent part is the comments for TA and TH2, both of which indicate they are potential temperatures. :

! TA [K] : Potential temperature at 1st wrf level (absolute temp)

! TH2 [K] : Diagnostic potential temperature at 2 m

In subroutine urban, TA is also used in the following codes, which I think actual temperature should be used:

TAV=TA*(1.+0.61*QA)

TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)

HR=RHO*CP*CHR*UA*(TRP-TA)*100.

TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP

XXX = 0.4*9.81*Z*TST/TA/UST/UST

It appears that the only variable T3D is the actual temperature. TA is potential temperature. Rest of the variables such as TH2, TS, etc. calculated using TA are potential temperatures. T2 in surface_driver.F using TH2_URB2D is actual temperature. In the subsequent line/command (3393) TH2 is again converted to the potential temperature.

Let's go by the definition of different temperature measures, in order to calculate the virtual potential temperature (TAV), the potential temperature is used. Hence TA has to be the potential temperature not the actual temperature. For simulations over Riyadh, if you do not consider TH2_URB2D as potential temperature, I believe you are then looking at potential temperature instead of actual temperature (which is T2). Line/command 587 in module_sf_noaddrv.F mentioned TH2_URB2D as potential temperature. @xuelingbo how did you modified the TH2_URB2D to actual temperature?

Also, to calculate the Q2, I believe that @DanLi-BU explanation on using CHS from Noah for calculating the QS in SLUCM make sense as the pervious and impervious surfaces have their own contribution. Your views @cenlinhe, @tslin2

@tslin2
Copy link

tslin2 commented Sep 4, 2024

@xuelingbo I took a quick look at the modifications in HRLDAS. You are right that this seems consistent with the second approach I described (diagnosing T2/TS for each patch separately and then aggregate them). However, the counter-argument is if so, then why don't we diagnose T2/TS for roof and canyon separately? Remember, in the current SLUCM, the fluxes from roof and canyon are calculated separately. In other words, the roof and canyon are treated as different patches.

Another way of looking at this issue is to ask the question why is the calculation of CHS/CHS2 only dependent on Z0HC, or why is Z0HR playing no role in the calculation of CHS/CHS2? We are trying to diagnose a T2/TS for the entire impervious land which composes of a roof and a canyon, I can't think of a reason why Z0HC should be used to represent the heat transfer capability of the entire impervious land. Note here we are talking about heat transfer, not momentum transfer. Maybe the role of roof can be neglected for momentum transfer, but I'm not sure about heat transfer.

I don't have good answers to these questions though :)

related to this, NCAR/hrldas#79

I did test this by uncomment

! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]

It did change the T2 results if using radiative temperature for TS; resulted in a higher temperature for both TS and T2. This method is more similar to the Noah-MP and BEP/BEM methods for calculating the TSK using the Stefan-Boltzmann law, but not exactly the same.

NoahMP-SLUCM:
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3519
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3527

NoahMP-BEP/BEM:
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3740-L3744

@tslin2
Copy link

tslin2 commented Sep 4, 2024

@joshi994 I agree with you that SLUCM needs to use potential temperature in module_sf_urban.F, but current TA in SLUCM is read from WRF 1st layer actual atmospheric temperature (T3D/T_PHY), not potential temperature. So the codes need to read the TH_PHY, which is the potential temperature in WRF. Just need to use actual temperature for some calculations in SLUCM, which can be converted based on potential temperature. The potential temperature was used in BEP/BEM for sensible heat flux calculation as well and converted to actual temperature for other purposes, such as calculating saturated water vapor pressure.

This involves another issue that @DanLi-BU and @chenghaow raised about the potential temperature in the Noah and Noah-MP. Previously, I did test it in a case using NoahMP-BEP by simply using the potential temperature to diagnose T2 for all grids, like BEP/BEM for T2 for urban grids in module_surface_driver.F, but this is an indirect way but not the best as Noah-MP used actual temperature, not potential temperature. I think using potential temperature is more suitable for urban models.

T2(I,J) = TH_PHY(i,1,j)/((1.E5/PSFC(I,J))**RCP) !urban

@xuelingbo
Copy link
Author

xuelingbo commented Sep 5, 2024

@xuelingbo I took a quick look at the modifications in HRLDAS. You are right that this seems consistent with the second approach I described (diagnosing T2/TS for each patch separately and then aggregate them). However, the counter-argument is if so, then why don't we diagnose T2/TS for roof and canyon separately? Remember, in the current SLUCM, the fluxes from roof and canyon are calculated separately. In other words, the roof and canyon are treated as different patches.
Another way of looking at this issue is to ask the question why is the calculation of CHS/CHS2 only dependent on Z0HC, or why is Z0HR playing no role in the calculation of CHS/CHS2? We are trying to diagnose a T2/TS for the entire impervious land which composes of a roof and a canyon, I can't think of a reason why Z0HC should be used to represent the heat transfer capability of the entire impervious land. Note here we are talking about heat transfer, not momentum transfer. Maybe the role of roof can be neglected for momentum transfer, but I'm not sure about heat transfer.
I don't have good answers to these questions though :)

related to this, NCAR/hrldas#79

I did test this by uncomment

! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]

It did change the T2 results if using radiative temperature for TS; resulted in a higher temperature for both TS and T2. This method is more similar to the Noah-MP and BEP/BEM methods for calculating the TSK using the Stefan-Boltzmann law, but not exactly the same.
NoahMP-SLUCM: https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3519 https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3527

NoahMP-BEP/BEM: https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L3740-L3744

Actually, in my research (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2023JD040531), for the TSK issue, I uncomment this line

! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]
, so the ouput TSK is radiative temeprature, but the calculation of TH2 is before this line,
TH2 = TS + (TA-TS) *(CHS/CHS2)
, so the value of TH2_URB2D will not be affected.

@joshi994
Copy link
Contributor

joshi994 commented Sep 6, 2024

@joshi994 I agree with you that SLUCM needs to use potential temperature in module_sf_urban.F, but current TA in SLUCM is read from WRF 1st layer actual atmospheric temperature (T3D/T_PHY), not potential temperature. So the codes need to read the TH_PHY, which is the potential temperature in WRF. Just need to use actual temperature for some calculations in SLUCM, which can be converted based on potential temperature. The potential temperature was used in BEP/BEM for sensible heat flux calculation as well and converted to actual temperature for other purposes, such as calculating saturated water vapor pressure.

This involves another issue that @DanLi-BU and @chenghaow raised about the potential temperature in the Noah and Noah-MP. Previously, I did test it in a case using NoahMP-BEP by simply using the potential temperature to diagnose T2 for all grids, like BEP/BEM for T2 for urban grids in module_surface_driver.F, but this is an indirect way but not the best as Noah-MP used actual temperature, not potential temperature. I think using potential temperature is more suitable for urban models.

T2(I,J) = TH_PHY(i,1,j)/((1.E5/PSFC(I,J))**RCP) !urban

I digged deep into the governing flux-form equation that WRF solver solves. It appears that the output temperature form the WRF model itself is potential temperature and not the actual temperature. In addition, the WRF output list indicate that T2 is actual temperature (converted using the Clausius-clapeyron equation) and TH2 is infact potential temperature. Please refer to pages 8-9 (specifically eqs. 2.7-2.14 and the text in between.
WRF user guide.pdf
11
22
https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/output.html#wrf-output-fields
https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/output_variables.html

As per the govenering equations T3D = T (perturbation potential temperature) + Reference potential temperature, whereas t_phy is converted to physical/actual temperature by thermodynamics relations. For instance, t_phy is supplied in module SLAB but it is converted to potential temperature within the module following the use of potential temperature. However, when T3D is supplied directly, it is actually the potential temperature. Hence, there's a difference between T3D and T_PHY and WEF solves for potential temperature and actual/physical temperature is obtained later.
@cenlinhe @xuelingbo @DanLi-BU @chenghaow

@tslin2
Copy link

tslin2 commented Sep 6, 2024

@joshi994 ,
T from the wrfoutput file is perturbation potential temperature theta-t0

state real th_phy_m_t0 ikj dyn_em 1 - irhd "t" "perturbation potential temperature theta-t0" "K"

but T3D is different from T,
T3D in module_sf_noahmpdrv.F
https://github.com/NCAR/noahmp/blob/848f54ad3d28c4303151fe5ad83724e232694422/drivers/wrf/module_sf_noahmpdrv.F#L23

is T_PHY
in module_surface_driver.F

T_PHY, QV_CURR, U_PHY, V_PHY, SWDOWN, SWDDIR, &

@joshi994
Copy link
Contributor

I also looked into the "module_pbl_driver.F" (https://github.com/wrf-model/WRF/blob/master/phys/module_pbl_driver.F) where it says:
CALL ysu( &
U3D=u_phytmp,V3D=v_phytmp,T3D=t_phy

however, when I looked into the YSU PBL scheme (module_bl_ysu.F), I found the following:
call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
,tx=t3d(ims,kms,j)
and later
thx(i,k) = tx(i,k)/pi2d(i,k) where pi2d is the Exner function which is basically the ration of (P2/P1)^(R/C_p). It suggests that the ysu pbl scheme is treating t3d as a potential temperature rather than the actual temperature. Following the comments from @tslin2 if t3d/t_phy are the same and denote actual temperature then the temperature used in PBL schemes needs to be looked at again.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

6 participants