Mismatch between calculated relatively humidity and CLM hist output

[yifanc17]

What version of the code are you using?
CESM 5.2.015


Have you made any changes to files in the source tree?
No


Describe your problem or question:
I forgot to include RH2M_U in CLM hist output in one of my land-only simulations, thus I want to use PBOT, TSA_U and Q2M_U to derive RH2M_U. I tested the manual calculation on another CLM hist output with both Q2M_U and RH2M_U available.

After checking the CTSM source code (QSatMod.F90, HumanIndexMod.F90, and UrbanFluxesMod.F90), I found two different implementations of the saturation specific humidity calculation:

• QSat(T, p, qs, es, qsdT, esdT) in QSatMod.F90, which is used in UrbanFluxesMod.F90 for calculating RH2M and RH2M_U.
• QSat_2(T_k, p_t, es_mb, de_mbdT, dlnes_mbdT, rs, rsdT, foftk, fdT) in HumanIndexMod.F90, which seems to be used only for diagnostic variables like wet bulb temperature.

I tried to translate both fortran subroutines into python functions (Qsat and Qsat_2 below) to reproduce the RH2M_U values. However, I found neither Qsat nor Qsat_2 can fully reproduce the model output for RH2M_U. The difference maps below (calculated minus model output) showed that Qsat_2 gives results closer to the model. My input variables' ranges are:

• PBOT_U: [50707.625, 106931.5] Pa
• TSA_U: [227.73, 331.94] K (approximately -45.4 to 58.8°C)
• Q2M_U: [4.0e-5, 0.0307] kg/kg
  which all seemed to be within valid physical limits and within the coefficient ranges (-75~100°C) used in the polynomial fits.

Thus I am not sure whether the discrepancy is due to issues in my python code. Additionally, do you think that QSat_2 is more appropriate for urban canopy cases? Any insights into how RH is computed in CLM for urban area or any tips for reproducing RH would be greatly appreciated! Thanks!

Below is my python code:

def Qsat(T_K,p_Pa):
    #water(T>=0°C)
    a=[6.11213476,0.444007856,0.143064234e-01,0.264461437e-03,
    0.305903558e-05,0.196237241e-07,0.892344772e-10,
    -0.373208410e-12,0.209339997e-15]
    
    #ice(-75°C<=T<0°C)
    c=[6.11123516,0.503109514,0.188369801e-0,0.420547422e-03,
    0.614396778e-05,0.602780717e-07,0.387940929e-09,
    0.149436277e-11,0.262655803e-14]
    
    T_C=T_K-273.15
    es_hPa=xr.where(
        T_C>=0.0,
        a[0]+T_C*(a[1]+T_C*(a[2]+T_C*(a[3]+T_C*(a[4]+T_C*(a[5]+
        T_C*(a[6]+T_C*(a[7]+T_C*a[8]))))))),
        c[0]+T_C*(c[1]+T_C*(c[2]+T_C*(c[3]+T_C*(c[4]+T_C*(c[5]+
        T_C*(c[6]+T_C*(c[7]+T_C*c[8])))))))
    )
    es_Pa=es_hPa*100
    qs=0.622*es_Pa/(p_Pa-0.378*es_Pa)
    return qs

def QSat_2(T_K, p_Pa):
    lambd_a=3.504
    alpha=17.67
    beta=243.5
    epsilon=0.6220
    es_C=6.112
    vkp=0.2854
    Cf=273.15
    refpres=1000.0
    
    p_tmb=p_Pa*0.01 # unit:mb
    tcfbdiff=T_K-Cf+beta
    es_mb=es_C*np.exp(alpha*(T_K-Cf)/tcfbdiff)
    ndimpress=(p_tmb/refpres)**vkp
    p0ndplam=refpres*ndimpress**lambd_a
    rs=epsilon*es_mb/(p0ndplam-es_mb) 
    
    return rs

temp_qs_U_Qsat=Qsat(model_TSA_U,model_PBOT)
temp_qs_U_Qsat2=QSat_2(model_TSA_U,model_PBOT)
calc_model_RH2M_U_Qsat=xr.where(model_Q2M_U/temp_qs_U_Qsat*100>100,100,model_Q2M_U/temp_qs_U_Qsat*100) #%
calc_model_RH2M_U_Qsat2=xr.where(model_Q2M_U/temp_qs_U_Qsat2*100>100,100,model_Q2M_U/temp_qs_U_Qsat2*100) #%

I have also attached the original .F90 code below (can't paste it here due to length limitation).

 

[oleson]

Hi Yifan, I see this in your ice coefficients:

0.188369801e-0

but the code has

0.188369801e-01

Since the differences you are seeing appear to be at high latitudes maybe this is part or all of the problem?

The QSat code is consistent with how the rest of the model, e.g., vegetated canopies, calculates saturated specific humidity.

The QSat_2 is advertised as being more accurate and is used in the Wet_Bulb subroutine. However, the Wet_Bulb subroutine is only invoked if the full set of human stress indices is selected (calc_human_stress_indices = 'ALL'). In that case it replaces the standard Stull version of Wet_Bulb (Wet_BulbS). The QSat_2 has not been ported to the rest of the model landunit types (vegetated, lake, etc.). It is also a bit slower than QSat_2.

 

[yifanc17]

Hi Yifan, I see this in your ice coefficients:

0.188369801e-0

but the code has

0.188369801e-01

Since the differences you are seeing appear to be at high latitudes maybe this is part or all of the problem?

The QSat code is consistent with how the rest of the model, e.g., vegetated canopies, calculates saturated specific humidity.

The QSat_2 is advertised as being more accurate and is used in the Wet_Bulb subroutine. However, the Wet_Bulb subroutine is only invoked if the full set of human stress indices is selected (calc_human_stress_indices = 'ALL'). In that case it replaces the standard Stull version of Wet_Bulb (Wet_BulbS). The QSat_2 has not been ported to the rest of the model landunit types (vegetated, lake, etc.). It is also a bit slower than QSat_2.

Hi Keith,

Yes thank you! Sorry I should've checked my functions more carefully. It was that coffecient causing the discrepancy in Qsat, now the difference (%) is minimal. And thanks for the clarification on Qsat_2 versus Qsat, that's really helpful!