refactoring ERA5 data mode

[uturuncoglu]

Description of changes

This PR aims to restructure ERA5 data mode to support:
(1) Using ERA5 as a replacement of JRA55 under CESM framework. This work does not include any necessary change needs to be done in CIME interface.
(2) Adding new data component configuration option for data atmosphere (DATM) as skip_field_check. The default value is .false. and does not change the exiting behavior. By setting it .true. the field check part is skipped and user is able to provide just set of fields supported by the data mode.
(3) The changes related to support UFS WM RTs with help of (2) and also some other changes.

Specific notes

Contributors other than yourself, if any:

CDEPS Issues Fixed (include github issue #):

Are there dependencies on other component PRs (if so list):
No

Are changes expected to change answers (bfb, different to roundoff, more substantial):
No

Any User Interface Changes (namelist or namelist defaults changes):
New optional skip_field_check option for datm

Testing performed (e.g. aux_cdeps, CESM prealpha, etc):
UFS WM RTs:

datm_cdeps_lnd_era5
datm_cdeps_lnd_era5_rst
hafs_regional_datm_cdeps

UFS Coastal:

DATM-DOCN-CICE6 and DATM-SCHSIM-CICE6 regional sea-ice configuration

Hashes used for testing:

[uturuncoglu]

@billsacks I am still working on it. Once tests are passed I'll let you know.
@NickSzapiro-NOAA This will also address the failed land RTs. skip_field_check needs to be used for those. I'll test them and let you know when all UFS WM tests are passed. There are some changes also in the UFS WM level CDEPS CMake interface since some files are added and removed from CDEPS. Once UFS WM sync CDEPS, we could address those issues.

[NickSzapiro-NOAA]

Hi @uturuncoglu @billsacks . May I ask if this is still planned?

I'm interested so I ran some tests in UFS.
datm_cdeps_lnd_era5_intel RT has answer changes but runs.
This would be nice to get EMC/CDEPS sync to ESCOMP unstuck.

There is also interest in using ERA5 for datm-ocn-ice work, like for RTOFS development.

fyi, AI came back with some review notes (NickSzapiro-NOAA#11 (review)).
Please let me know if you'd like contributions

[sanAkel]

@uturuncoglu I got a link to this (draft) PR from @NickSzapiro-NOAA in the context of UFS DATM with the "data" atmosphere coming from ERA-5. 2 questions for either/both of you:

1. What would be a good ETA for the PR to go from draft ➡️ ready? And eventually merged.
2. Reg the splitting of (net) downward shortwave into the four bands, would it be possible to please not hard code those fractions? Suggestion: read via yaml, default to these values. or to those in ACCESS-NRI.

[uturuncoglu]

@sanAkel @NickSzapiro-NOAA The only thing that I need to do is to add new formulas to valuate saturation mixing ratio. I am plaining to add new option for it. We could also introduce namelist option to get those fractions from datm_in. I am also taking off starting from next week and back to US end of the month. So, I am not sure if I could do everything before love since I am also trying to finalize things for other projects but if I could do, I'll update you through the PR and and mark tis as ready to review.

[sanAkel]

@sanAkel @NickSzapiro-NOAA The only thing that I need to do is to add new formulas to valuate saturation mixing ratio. I am plaining to add new option for it. We could also introduce namelist option to get those fractions from datm_in. I am also taking off starting from next week and back to US end of the month. So, I am not sure if I could do everything before love since I am also trying to finalize things for other projects but if I could do, I'll update you through the PR and and mark tis as ready to review.

Thanks @uturuncoglu for ⬆️ update. If you have it ready before your leave, great! If not, we will wait or make do whatever we got!

[sanAkel]

I am also taking off starting from next week and back to US end of the month.

@uturuncoglu are you back? Hope you had a nice break!

Checking if you have time to finish this work? 🙏

[uturuncoglu]

@sanAkel Yes, I am back. I am trying to catch other things but I could work on this next week and we could finalize it.

[uturuncoglu]

@sanAkel BTW, i am assuming this is tested without any issue in your end without the remaining additional things.

[sanAkel]

@sanAkel BTW, i am assuming this is tested without any issue in your end without the remaining additional things.

@uturuncoglu
Glad to hear all went well with your trip, welcome back!

@raphaeldussin is going to set up a DATM run with the ERA5 forcing.

If you happen to have a set of sample forcing files (from the ERA5 reanalysis), that you can point us to, we can test it at our end.

Thanks again!

[uturuncoglu]

@sanAkel i am trying the finalize the implementation but I think new algorithm to calculate qsat is wrong and give negative and very large values. I tested following code under Python like following,

qsat = compute_q_sat(300, 102315)

and it is giving -2.6563530902204677. The temperature is in K and pressure is in Pa in this case. Unless, if there is a unit issue in here, there must be some issue with the codes. I also tried to use it with CDEPS and I am getting weird results.

def compute_e_sw(t):
    """Saturation vapor pressure over water (hPa)"""
    return 10 ** ((0.7859 + 0.03477 * t) / (1 + 0.00412 * t))

def compute_f_w(t, p):
    """Pressure correction factor f_w (unitless)"""
    return 1 + 1e-6 * p * (4.5 + 0.0006 * t**2)

def compute_q_sat(t, p):
    """Saturation specific humidity (unitless)"""
    correction_factor_salt_water = 0.98 # not used in Gill but used in Tsujino et al (2018)
    e_sw = compute_e_sw(t)
    f_w = compute_f_w(t, p)
    e_sw_prime = f_w * e_sw *correction_factor_salt_water
    omega = 0.62197
    return (omega * e_sw_prime) / (p - (1 - omega) * e_sw_prime)

[NickSzapiro-NOAA]

AI suggests it is a units issue ... ie

def compute_e_sw(t_c):
    """Saturation vapor pressure over water (hPa). Expects Celsius."""
    return 10 ** ((0.7859 + 0.03477 * t_c) / (1 + 0.00412 * t_c))

def compute_f_w(t_c, p_hpa):
    """Pressure correction factor f_w (unitless). Expects Celsius and hPa."""
    return 1 + 1e-6 * p_hpa * (4.5 + 0.0006 * t_c**2)

def compute_q_sat(t_k, p_pa):
    """Saturation specific humidity (kg/kg). Inputs: Kelvin and Pascals."""
    
    # 1. Convert inputs to empirical requirements (Celsius and hPa/mb)
    t_c = t_k - 273.15
    p_hpa = p_pa / 100.0
    
    # 2. Compute empirical terms
    correction_factor_salt_water = 0.98 
    e_sw = compute_e_sw(t_c)
    f_w = compute_f_w(t_c, p_hpa)
    
    e_sw_prime_hpa = f_w * e_sw * correction_factor_salt_water
    omega = 0.62197
    
    # 3. Compute specific humidity
    # p_hpa and e_sw_prime_hpa are now both in hPa
    q_sat = (omega * e_sw_prime_hpa) / (p_hpa - (1 - omega) * e_sw_prime_hpa)
    
    return q_sat

[raphaeldussin]

I used the Bolton approximation (which is fairly accurate at the observed range of temperatures), you can compare with my code here https://github.com/raphaeldussin/newERAtools/blob/e29d89707058cca52abb5d8015a85894756f9c37/compute_q2_ERA5.f90#L154

[sanAkel]

@sanAkel i am trying the finalize the implementation but I think new algorithm to calculate qsat is wrong and give negative and very large values. I tested following code under Python like following,

qsat = compute_q_sat(300, 102315)

and it is giving -2.6563530902204677. The temperature is in K and pressure is in Pa in this case. Unless, if there is a unit issue in here, there must be some issue with the codes. I also tried to use it with CDEPS and I am getting weird results.

def compute_e_sw(t):
    """Saturation vapor pressure over water (hPa)"""
    return 10 ** ((0.7859 + 0.03477 * t) / (1 + 0.00412 * t))

def compute_f_w(t, p):
    """Pressure correction factor f_w (unitless)"""
    return 1 + 1e-6 * p * (4.5 + 0.0006 * t**2)

def compute_q_sat(t, p):
    """Saturation specific humidity (unitless)"""
    correction_factor_salt_water = 0.98 # not used in Gill but used in Tsujino et al (2018)
    e_sw = compute_e_sw(t)
    f_w = compute_f_w(t, p)
    e_sw_prime = f_w * e_sw *correction_factor_salt_water
    omega = 0.62197
    return (omega * e_sw_prime) / (p - (1 - omega) * e_sw_prime)

@uturuncoglu Please see ⬆️ response from @raphaeldussin.

If your error persists, please check with @raphaeldussin.

[uturuncoglu]

@sanAkel @raphaeldussin @NickSzapiro-NOAA I fixed the unit issue and pushed it to the fork. Now you could select type of algorithm in datm_in using qsat_algorithm option and it gets 0 and 1 values and they are creating spatially similar results but option 1 values are slightly higher as you can see from the plot. The left plot is created from option 0 (old version and default) and the right one is option 1 (new implementation).

Anyway, I tested this with UFS WM datm_cdeps_lnd_era5 RT without any issue (the test is passing). @NickSzapiro-NOAA if you don't mind could you also check in your end (you need to modify datm_in.IN in param folder and add skip_field_check = .true.). Other please also check in your end and let me know if you need any further change. You could create PR to my fork and branch for any change. If you are okay, I'll make this PR ready for review and sync with head.

[uturuncoglu]

I also sync the branch with main.

[sanAkel]

That's awesome @uturuncoglu 👏 !

@NickSzapiro-NOAA Can we run the same RT test @uturuncoglu ran on wcoss-2? To rephrase my question, do we have the same datasets for that test available?

[NickSzapiro-NOAA]

Thanks @uturuncoglu . I will start some tests tonight

The same RT input data should be available on all platforms @sanAkel
https://ufs-weather-model.readthedocs.io/en/develop/BuildingAndRunning.html#get-data
but note that datm_cdeps_lnd_era5 is a "land" RT

[sanAkel]

Thanks @uturuncoglu . I will start some tests tonight

Thanks @NickSzapiro-NOAA

The same RT input data should be available on all platforms @sanAkel https://ufs-weather-model.readthedocs.io/en/develop/BuildingAndRunning.html#get-data but note that datm_cdeps_lnd_era5 is a "land" RT

What do we need to do make into an "ocean" RT? Please list steps.. 🙏

[NickSzapiro-NOAA]

What do we need to do make into an "ocean" RT

Model configuration and input data with forcing, ICs, and fix files.
This PR lets CDEPS datm in UFS serve ERA5 to ocean-ice components again (if all goes well)

I think there is quite some overlap with @uturuncoglu in datm-ocean-ice efforts, but I'm not sure this ESCOMP/CDEPS thread is right place to continue this topic. Happy to continue in ufs-community/ufs-weather-model#3294 , NOAA-EMC/RTOFS_GLO#147 , or wherever else as well

[sdurski]

To re-reiterate the ERA5 precipitation issue I mentioned at the UFS-coastal technical tag up today:

This may be rather particular to our current DATA-CICE and CICE-SCHISM setups but,

If ERA5 provides precipitation in units of meters (cumulative meters of rain per hour), rather than kg/m^2/s, it may be misinterpreted as a rate leading to 3.6X more precip than is actually occurring.

(Code/line numbers are from datm/datm_datamode_era5_mod.F90 as of head e301493)
precmax is set once, on first_timestep , from the rain field's global domain-max L519-525:
then
681 | ! convert m to kg/m^2/s
682 | if (precmax < 0.01_r8) then
683 | if (mainproc) write(logunit,) trim(subname),' precipitation related variables are already in kg/m^2/s unit!'
684 | else
685 | if (associated(Faxa_rain)) Faxa_rain(:) = Faxa_rain(:)/3600.0_r8rhofw
686 | if (associated(Faxa_rainc)) Faxa_rainc(:) = Faxa_rainc(:)/3600.0_r8rhofw
687 | if (associated(Faxa_rainl)) Faxa_rainl(:) = Faxa_rainl(:)/3600.0_r8rhofw
688 | if (associated(Faxa_snowc)) Faxa_snowc(:) = Faxa_snowc(:)/3600.0_r8rhofw
689 | if (associated(Faxa_snowl)) Faxa_snowl(:) = Faxa_snowl(:)/3600.0_r8rhofw
690 | end if

But magnitude alone can't separate the two conventions the L682 test is choosing between. For true precip R mm/hr the file value is R/1000 as metres-of-1-hour-accumulation vs R/3600 as kg/m²/s — they differ by only 3.6×, and both sit in the same 0…~10⁻² range for ordinary precip. So the 0.01 cutoff effectively asks only "is there a ≥10 mm/hr cell somewhere in the domain at the first timestep?"

Model reports ""already in kg/m^2/s unit!" to output log on L683, even though CDEPS doesn't actually check the units in the ERA5 input file.

In our SCHISM CICE test case (prior to the implementation of Faxa_snow) this led to far too much snow falling on the sea ice, resulting in too thick ice that persisted too long into the spring.

NOTE ALSO currently Faxa_snow is left out of the above conversion (L685-689)

I think ideally CDEPS should read in units from the ERA5 file explicitly and only if units are not found should it make a guess. Barring that the guess should be reported in the log as an assumption that has some ambiguity associated with it.

[NickSzapiro-NOAA]

How do you decide which export state variables are initialized (to zero)?

Say, what about taux, tauy?

[uturuncoglu]

@NickSzapiro-NOAA As i see the fields are created in dshr_fldlist_realize() call in dshr/dshr_fldlist_mod.F90 and it seems that there is no any initial value assigned to them. Maybe ESMF is initializing them internally but I am not %100 sure about it. We could use ESMF_FieldFill to set zero initially in that level to be in the safe side. @billsacks let me know what do you think?

[billsacks]

I did a very quick look through other CDEPS components and it looks like this kind of zeroing generally isn't done. I haven't looked at this ERA5 datamode enough to know if it might be needed for this one specially. My initial thought is that these variables should be set in the call to datm_datamode_era5_advance, so there's probably no need to initially zero them, but you shouldn't take my word for it because I'm not very familiar with this.

[NickSzapiro-NOAA]

Quick note that tests pass on ursa and gaea.c6 .
datm RTs are too slow on wcoss2 in general.

But will need some response to #421 before this could be merged to EMC/CDEPS.

[uturuncoglu]

@NickSzapiro-NOAA It is good to know that tests are passing. Thanks for confirmation.

[uturuncoglu]

@sdurski regarding to the issue related with the unit conversion. I am not sure what needs to be done. I think the best approach is to assume that the stream will be on ERA5 original units. The we could remove the range check from all the fields and just apply unit conversion in all the cases. So, it would be the user responsibility to provide information in expected unit. This is same for all the data modes.

[uturuncoglu]

@sdurski The units for variable from ERA5 side are:

• m for precipitation will be converted to kg/m^2/s in all the base
• J/m2 for heat flux variables and will be converted to W/m2
• N/m^2 s (Pa ?) for momentum fluxes and will be converted to N/m2

So, we need to confirm those and I could remove the check for ranges from all the variables.

[sdurski]

@uturuncoglu The era5 file that I have been using for SCHISM-CICE runs and which was used for Joey's data-atm/data-ocn/CICE runs has the following units for variables

Variable Units
glbrad W m⁻² downward shortwave (global) radiation
dlwsfc W m⁻² downward longwave at surface
wndewd m s⁻¹ eastward 10 m wind
wndnwd m s⁻¹ northward 10 m wind
airtmp K 2 m air temperature
spchmd (none in file) — specific humidity, kg kg⁻¹ 2 m specific humidity
ttlpcp m total precipitation (signed: negative = snow, positive = rain)
Pair Pa surface / mean-sea-level pressure

I recently changed the variable ttlpcp in the input file to be a rate in kg m⁻² s⁻¹ to prevent the mis-identification in CDEPS. But both the precipitation rate (mtpr) and the total precipitation (tp) are available from Copernicus when extracting era5 forcing data. So original units can be ambiguous. If we're hardwiring an assumption in ttlpcp, meters matches the variable name better.

For the short wave and longwave there's a similar ambiguity. You can grab
msdwswrf → glbrad Mean surface downward short-wave radiation flux W m⁻²
msdwlwrf → dlwsfc Mean surface downward long-wave radiation flux W m⁻²

or
ssrd (surface solar radiation downwards) — J m⁻², 1-hour accumulation
strd (surface thermal radiation downwards) — J m⁻², 1-hour accumulation

Is CDEPS currently assuming J/m^2 and always converting? It seems like my W/m^2 values are getting through ok.

I think the answer is different if we're trying to set something up for any anonymous user who grabs their own ERA5 forcing, or if we are just trying to make sure CDEPS works for the current tests that involve ERA5. If the latter, just let me know what units you settle on and I can adjust my forcing file accordingly.

[uturuncoglu]

@sdurski I think best approach is to fix the units in CDEPS side and force the user to use that convention. If they want to bring another variable with another units, they need to do it in the CDEPS side with some customization (maybe they could create their own ERA5 mode). Current, design does not allow to be generic enough and cover all the cases. @sanAkel Let me your unit preference too. We need to find a common ground between UFS, UFS Coastal and CESM and fix the units to move on in this PR.

[sanAkel]

@sanAkel Let me your unit preference too. We need to find a common ground between UFS, UFS Coastal and CESM and fix the units to move on in this PR.

@uturuncoglu I will let @raphaeldussin respond to this question. 🙏