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

modify create_time_series to save RESTOM time series file #240

Merged
merged 10 commits into from
Jan 9, 2024
30 changes: 30 additions & 0 deletions lib/adf_diag.py
Original file line number Diff line number Diff line change
Expand Up @@ -433,6 +433,11 @@ def call_ncrcat(cmd):
#Note: could use `open_mfdataset`, but that can become very slow;
# This approach effectively assumes that all files contain the same variables.

#Check if the files contain information for RESTOM:
lwnet_name = "FLNT"
swnet_name = "FSNT"
has_restom_data = (lwnet_name in hist_file_var_list) and (swnet_name in hist_file_var_list) #RESTOM = FSNT-FLNT

#Check what kind of vertical coordinate (if any) is being used for this model run:
#------------------------
if 'lev' in hist_file_ds:
Expand Down Expand Up @@ -585,6 +590,31 @@ def call_ncrcat(cmd):
_ = mpool.map(call_ncrcat, list_of_commands)
#End with

#Make a time series file for RESTOM
if has_restom_data:
import netCDF4 # netCDF4.default_fillvals is a dict of default fill values
longwavenet = ts_dir[case_idx] + os.sep + \
".".join([case_name, hist_str, lwnet_name, time_string, "nc" ])
shortwavenet = ts_dir[case_idx] + os.sep + \
".".join([case_name, hist_str, swnet_name, time_string, "nc" ])
restom_fil = ts_dir[case_idx] + os.sep + \
".".join([case_name, hist_str, "RESTOM", time_string, "nc" ])
# combine into a single file
tmp1 = xr.open_dataset(longwavenet)
tmp2 = xr.open_dataset(shortwavenet)
tmp3 = tmp2[swnet_name]-tmp1[lwnet_name]
tmp3.name = "RESTOM"
tmp3.attrs['long_name'] = "residual radiative flux at top of model"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would just add the units for RESTOM here too

tmp3.attrs['units'] = "W/m2"

Copy link
Collaborator Author

@brianpm brianpm Jan 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These lines are gone now because @nusbaume asked for me to use the implementation following PRECT.

enc = {}
if not isinstance(tmp3, xr.Dataset):
tmp3 = tmp3.to_dataset()
for dv in tmp3.data_vars:
if tmp3[dv].dtype == 'float32':
enc[dv] = {"_FillValue":netCDF4.default_fillvals['f8'],
"zlib": False}
for cv in tmp3.coords:
enc[cv] = {'zlib': False, '_FillValue': None}
tmp3.to_netcdf(restom_fil, encoding=enc)
#End cases loop

#Notify user that script has ended:
Expand Down