-
Notifications
You must be signed in to change notification settings - Fork 0
/
load_1d_model_output.py
400 lines (382 loc) · 24 KB
/
load_1d_model_output.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
import numpy as np
import xarray as xr
import pandas as pd
import glob, os
"""
-----------------------------------------------------------------------------------
Israel Silber
Last update: 7/16/2020
----------------------------------------------------------------------------------
Current methods:
load_1d_model_output - load multiple fields from multiple output types onto an
Xarray dataset object.
init_load_1d_model_output - initialize model output metadata (should match the
1D model code - updated as of 7/16/2020.
"""
#---------------------------------------------------------------------------------------------------------------
def load_1d_model_output(run_path, time_range=[np.nan, np.nan],
types_to_load=["prof", "flx", "tke", "mc", "sfc"], fields_to_load={}, **kwargs):
"""
----------------------
Israel Silber
Last update: 7/16/2020
----------------------
This method loads (Golaz) 1D model output. The method has default output lists and
descriptions, but can accept different output parameter sorting or lists (shoudl
match the model output subroutines, all located in /model/initialize.f
Parameters
----------
run_path:
path where output files are located.
time_range:
2 element list providing the time range in fractional hours of simulation time
(set both to np.nan to load everything).
type_to_load (load when true):
list containing output types to load (currently working with prof, flx, tke,
mc, and sfc.
fields_to_load:
a dictionary containing field names to load for every requested output type
(load all fields if empty, if output type not in keys, or if output type values
are empty).
Returns
-------
Output_1D:
xarray dataset gathering all model output types (if requested)
"""
"""
os.chdir("/home/meteo/ixs34/Python/1d_model")
run_path = '/home/meteo/ixs34/golaz_1d/runs/jyhRuns/results/awr_drz/'
"""
if not kwargs:
out_info = init_load_1d_model_output()
else:
out_info = init_load_1d_model_output(kwargs)
# -------------------------------------load data------------------------------------
output_1d = {}
for out_type in types_to_load:
for filename in sorted(glob.glob(run_path + out_type + "*.dat")):
if (float(filename[-9:-6]) + float(filename[-6:-4])/60. >= time_range[0] and \
float(filename[-9:-6]) + float(filename[-6:-4])/60. < time_range[1]) or \
np.isnan(time_range[0]):
output_tmp = pd.read_csv(filename, delim_whitespace=True, header=None,
names=out_info[out_type]["names"])
# grab the model domain height array
if "zt" in output_tmp.keys():
hgt = output_tmp["zt"].values.tolist()
output_tmp = output_tmp.drop(columns=['zt'])
if "zm" in output_tmp.keys():
hgt = output_tmp["zm"].values.tolist()
output_tmp = output_tmp.drop(columns=['zm'])
if len(fields_to_load) > 0:
if out_type in fields_to_load.keys():
if len(fields_to_load[out_type]) > 0:
fields_to_filter = []
for ff in fields_to_load[out_type]:
if ff in output_tmp.keys() and ff is not "zt" and ff is not "zm":
fields_to_filter.append(ff)
else:
print(ff + " not in " + out_type + " output type")
if fields_to_filter: # filter empty lists
output_tmp = output_tmp.filter(items=fields_to_filter)
# convert from pandas to xarray add append to time dim
output_tmp = output_tmp.to_xarray()
output_tmp = output_tmp.rename({'index': 'height'})
output_tmp = output_tmp.expand_dims({"time": \
[float(filename[-9:-6]) + float(filename[-6:-4])/60.]}, axis=1)
output_tmp = output_tmp.assign_coords({"height": hgt})
# add attributes.
for nn in out_info[out_type].keys():
if nn is not "names":
for ff in out_info[out_type][nn].keys():
if ff in output_tmp.keys():
output_tmp[ff] = \
output_tmp[ff].assign_attrs({nn: out_info[out_type][nn][ff]})
# concatenate and arrange
if out_type not in output_1d.keys():
output_1d[out_type] = output_tmp
else:
output_1d[out_type] = xr.concat([output_1d[out_type], output_tmp], "time")
print("done loading " + filename)
del output_tmp
else:
print(filename + " outside the requested time range (%.2f-%.2f h)" \
% (tuple(time_range)))
return output_1d
#---------------------------------------------------------------------------------------------------------------
def init_load_1d_model_output(sim_version=True, **kwargs):
"""
This method checks whether output field["names"] other then the default were requested
and sets them accordingly (note that if names are different, "long_name" and "units"
should be different as well).
Parameters
----------
sim_version:
If True then using the simulator version mc output namelist order. Non-simulator
version should still be updated.
"""
out_info = {}
# -------------------------------------init prof------------------------------------
out_info["prof"] = {}
if 'out_prof_names' in kwargs.keys():
out_info["prof"]["names"] = kwargs['out_prof_names']
else:
out_info["prof"]["names"] = ["zm", "thetail", "u", "v", "Km", "Kh", "thv_flx", "TKE",
"Khh", "lx", "rt", "rl", "ri", "CF", "T", "p", "rho",
"thv", "htrt", "htrt_lw", "htrt_sw", "w", "hiv", "h"]
if 'out_prof_units' in kwargs.keys():
out_info["prof"]["units"] = kwargs['out_prof_units']
else:
out_info["prof"]["units"] = {"zm": "m", "thetail": "K", "u": "m/s", "v": "m/s", "Km": "",
"Kh": "", "thv_flx": "J m^{-3}", "TKE": "m^2 s^{-2}", "Khh": "",
"lx": "m", "rt": "kg/kg", "rl": "kg/kg", "ri": "kg/kg", "CF": "",
"T": "K", "p": "Pa", "rho": "kg m^{-3}",
"thv": "K", "htrt": "K/h", "htrt_lw": "K/h",
"htrt_sw": "K/h", "w": "m/s", "hiv": "J", "h": "J"}
if 'out_prof_long_name' in kwargs.keys():
out_info["prof"]["long_name"] = kwargs['out_prof_long_name']
else:
out_info["prof"]["long_name"] = {"zm": "Mean grid height",
"thetail": "Mean ice liquid potential temperature",
"u": "u wind", "v": "v wind", "Km": "Momentum eddy diffusivity",
"Kh": "Heat eddy diffusivity", "thv_flx": "Cp*rho*thvw",
"TKE": "Turbulent kinetik energy",
"Khh": "Coefficient for second order terms",
"lx": "Master length scale variable for Galperin's MY2.5 scheme",
"rt": "Total mixing ratio", "rl": "Liquid mixing ratio (cld+rain)",
"ri": "Ice mixing raito", "CF": "Fractional cloudiness",
"T": "Temperature", "p": "Pressure", "rho": "Density",
"thv": "Virtual potential temperature",
"htrt": "Radiative heating rate",
"htrt_lw": "LW radiative heating rate",
"htrt_sw": "SW radiative heating rate", "w": "vertical velocity",
"hiv": "ice vapor static energy", "h": "moist static energy"}
if 'out_prof_scale' in kwargs.keys():
out_info["prof"]["scale"] = kwargs['out_prof_scale']
else:
out_info["prof"]["scale"] = {"zm": "linear", "thetail": "linear", "u": "linear",
"v": "linear", "Km": "linear",
"Kh": "linear", "thv_flx": "linear", "TKE": "linear",
"Khh": "linear", "lx": "linear",
"rt": "log", "rl": "log", "ri": "log", "CF": "linear",
"T": "linear", "p": "log", "rho": "linear",
"thv": "linear", "htrt": "linear", "htrt_lw": "linear",
"htrt_sw": "linear", "w": "linear", "hiv": "linear",
"h": "linear"}
# -------------------------------------init flux------------------------------------
out_info["flx"] = {}
if 'out_flx_names' in kwargs.keys():
out_info["flx"]["names"] = kwargs['out_flx_names']
else:
out_info["flx"]["names"] = ["zt", "Km", "Kh", "Khh", "u_turb_flx", "v_turb_flx", "thil_heat_flx",
"rt_heat_flx", "thv_flx", "swfu", "swfd", "lwfu", "lwfd", "SST"]
if 'out_flx_units' in kwargs.keys():
out_info["flx"]["units"] = kwargs['out_flx_units']
else:
out_info["flx"]["units"] = {"zt": "m", "Km": "", "Kh": "", "Khh": "", "u_turb_flx": "s^{-1}",
"v_turb_flx": "s^{-1}", "thil_heat_flx": "J m^{-4}",
"rt_heat_flx": "J K^{-1} m^{-4}", "thv_flx": "J m^{-3}",
"swfu": "W m^{-2}", "swfd": "W m^{-2}", "lwfu": "W m^{-2}",
"lwfd": "W m^{-2}", "SST": "K"}
if 'out_flx_long_name' in kwargs.keys():
out_info["flx"]["long_name"] = kwargs['out_flx_long_name']
else:
out_info["flx"]["long_name"] = {"zt": "Turbulent grid heights", "Km": "Momentum eddy diffusivity",
"Kh": "Heat eddy diffusivity",
"Khh": "Coefficient for second order terms",
"u_turb_flx": "u momentum turbulent flux",
"v_turb_flx": "v momentum turbulent flux",
"thil_heat_flx": "thetail heat flux",
"rt_heat_flx": "rt heat flux",
"thv_flx": "Cp*rho*thvw", "swfu": "Upwelling SW flux",
"swfd": "Downwelling SW flux",
"lwfu": "Upwelling LW flux", "lwfd": "Downwelling LW flux",
"SST": "Surface temperature"}
if 'out_flx_scale' in kwargs.keys():
out_info["flx"]["scale"] = kwargs['out_flx_scale']
else:
out_info["flx"]["scale"] = {"zt": "linear", "Km": "linear", "Kh": "linear",
"Khh": "linear", "u_turb_flx": "linear",
"v_turb_flx": "linear", "thil_heat_flx": "linear",
"rt_heat_flx": "symlog", "thv_flx": "linear",
"swfu": "linear", "swfd": "linear", "lwfu": "linear",
"lwfd": "linear", "SST": "linear"}
# -------------------------------------init TKE------------------------------------
out_info["tke"] = {}
if 'out_tke_names' in kwargs.keys():
out_info["tke"]["names"] = kwargs['out_tke_names']
else:
out_info["tke"]["names"] = ["zt", "TKE", "tke_shear", "tke_buoy", "tke_ttrspt", "tke_vtrspt",
"tke_dissip"]
if 'out_tke_units' in kwargs.keys():
out_info["tke"]["units"] = kwargs['out_tke_units']
else:
out_info["tke"]["units"] = {"zt": "m", "TKE": "m^2 s^{-2}", "tke_shear": "m^2 s^{-2}",
"tke_buoy": "m^2 s^{-2}",
"tke_ttrspt": "m^2 s^{-2}", "tke_vtrspt": "m^2 s^{-2}", "tke_dissip": "m^2 s^{-2}"}
if 'out_tke_long_name' in kwargs.keys():
out_info["tke"]["long_name"] = kwargs['out_tke_long_name']
else:
out_info["tke"]["long_name"] = {"zt": "Turbulent grid heights", "TKE": "Turbulent kinetik energy",
"tke_shear": "TKE shear production",
"tke_buoy": "TKE buoyancy production",
"tke_ttrspt": "TKE turbulent transport",
"tke_vtrspt": "TKE vertical transport",
"tke_dissip": "TKE dissipation"}
if 'out_tke_scale' in kwargs.keys():
out_info["tke"]["scale"] = kwargs['out_tke_scale']
else:
out_info["tke"]["scale"] = {"zt": "linear", "TKE": "linear", "tke_shear": "linear",
"tke_buoy": "linear",
"tke_ttrspt": "linear", "tke_vtrspt": "linear",
"tke_dissip": "linear"}
# -------------------------------------init microphysics---------------------------------
out_info["mc"] = {}
if 'out_mc_names' in kwargs.keys():
out_info["mc"]["names"] = kwargs['out_mc_names']
elif sim_version is True:
out_info["mc"]["names"] = ["zm", "rt", "rvap", "rcld", "rrain", "rpice", "rsnow", "raggr",
"rgrau", "rhail", "cpice", "amean", "cmean", "rhoiavg", "aspect",
"vthabn", "vthabm", "dmdtvap", "ssi", "zdr", "zh",
"drpp_dep", "drpp_sedim", "drrp_sedim", "drrp_collect", "drt_horizAdv",
"drc_ncphys", "drv_mcphys", "drc_turb", "drv_turb", "dthil_subs", "drt_subs"]
else:
out_info["mc"]["names"] = ["zm", "rt", "rvap", "rcld", "rrain", "rpice", "rsnow", "raggr",
"rgrau", "rhail", "cpice", "amean", "cmean", "rhoiavg", "aspect",
"vthabn", "vthabm", "dmdtvap" ,"ssi", "drpp_dep", "drpp_sedim",
"drrp_sedim", "drrp_collect", "drt_horizAdv", "drc_ncphys",
"drv_mcphys", "drc_turb", "drv_turb", "dthil_subs", "drt_subs"]
if 'out_mc_units' in kwargs.keys():
out_info["mc"]["units"] = kwargs['out_mc_units']
elif sim_version is True:
out_info["mc"]["units"] = {"zm": "m", "rt": "kg/kg", "rvap": "kg/kg", "rcld": "kg/kg",
"rrain": "kg/kg", "rpice": "kg/kg", "rsnow": "kg/kg", "raggr": "kg/kg",
"rgrau": "kg/kg", "rhail": "kg/kg", "cpice": "m^{-3}", "amean": "\mu m",
"cmean": "\mu m", "rhoiavg": "kg m^{-3}", "aspect": "", "vthabn": "m/s",
"vthabm": "m/s", "dmdtvap": "-" ,"ssi": "-", "zdr": "dB", "zh": "dBZ",
"drpp_dep": "kg kg^{-1} s^{-1}", "drpp_sedim": "kg kg^{-1} s^{-1}",
"drrp_sedim": "kg kg^{-1} s^{-1}", "drrp_collect": "kg kg^{-1} s^{-1}",
"drt_horizAdv": "kg kg^{-1} s^{-1}", "drc_ncphys": "kg kg^{-1} s^{-1}",
"drv_mcphys": "kg kg^{-1} s^{-1}", "drc_turb": "kg kg^{-1} s^{-1}",
"drv_turb": "kg kg^{-1} s^{-1}", "dthil_subs": "K s^{-1}",
"drt_subs": "kg kg^{-1} s^{-1}"}
else:
out_info["mc"]["units"] = {"zm": "m", "rt": "kg/kg", "rvap": "kg/kg", "rcld": "kg/kg",
"rrain": "kg/kg", "rpice": "kg/kg", "rsnow": "kg/kg", "raggr": "kg/kg",
"rgrau": "kg/kg", "rhail": "kg/kg", "cpice": "m^{-3}", "amean": "\mu m",
"cmean": "\mu m", "rhoiavg": "kg m^{-3}", "aspect": "", "vthabn": "m/s",
"vthabm": "m/s", "dmdtvap": "-" ,"ssi": "-",
"drpp_dep": "kg kg^{-1} s^{-1}", "drpp_sedim": "kg kg^{-1} s^{-1}",
"drrp_sedim": "kg kg^{-1} s^{-1}", "drrp_collect": "kg kg^{-1} s^{-1}",
"drt_horizAdv": "kg kg^{-1} s^{-1}", "drc_ncphys": "kg kg^{-1} s^{-1}",
"drv_mcphys": "kg kg^{-1} s^{-1}", "drc_turb": "kg kg^{-1} s^{-1}",
"drv_turb": "kg kg^{-1} s^{-1}", "dthil_subs": "K s^{-1}",
"drt_subs": "kg kg^{-1} s^{-1}"}
if 'out_mc_long_name' in kwargs.keys():
out_info["mc"]["long_name"] = kwargs['out_mc_long_name']
elif sim_version is True:
out_info["mc"]["long_name"] = {"zm": "Mean grid height", "rt": "Total mixing ratio",
"rvap": "Vapor mixing ratio", "rcld": "Cloud water mixing ratio",
"rrain": "Rain water mixing ratio",
"rpice": "Pristine ice mixing ratio", "rsnow": "Snow mixing ratio",
"raggr": "Aggregate mixing raito",
"rgrau": "graupel mixing ratio", "rhail": "hail mixing ratio",
"cpice": "Ice number concentration",
"amean": "mean a-axis length", "cmean": "mean a-axis length",
"rhoiavg": "Mean ice density", "aspect": "Aspect ratio",
"vthabn": "Number weighted terminal fall speed",
"vthabm": "Mass weighted terminal velocity",
"dmdtvap": "Vapor growth rate", "ssi": "Ice supersaturation",
"zdr": "Differential reflectivity factor (Radar)",
"zh": "Horizontally-polarized reflectivity factor (Radar)",
"drpp_dep": "Tendencies for depositional growth of pristine ice",
"drpp_sedim": "Tendencies for sedimentation of pristine ice",
"drrp_sedim": "Tendencies for sedimentation of rain",
"drrp_collect": "Tendencies for collection of cloud into rain",
"drt_horizAdv": "Tendencies in rt due to horizontal advection",
"drc_ncphys": "Tendencies in rc due to microphysics",
"drv_mcphys": "Tendencies in rv due to microphysics",
"drc_turb": "Tendencies in rc due to turbulence",
"drv_turb": "Tendencies in rv due to turbulence",
"dthil_subs": "Tendencies in thetail due to subsidence",
"drt_subs": "Tendencies in rt due to subsidence"}
else:
out_info["mc"]["long_name"] = {"zm": "Mean grid height", "rt": "Total mixing ratio",
"rvap": "Vapor mixing ratio", "rcld": "Cloud water mixing ratio",
"rrain": "Rain water mixing ratio",
"rpice": "Pristine ice mixing ratio", "rsnow": "Snow mixing ratio",
"raggr": "Aggregate mixing raito",
"rgrau": "graupel mixing ratio", "rhail": "hail mixing ratio",
"cpice": "Ice number concentration",
"amean": "mean a-axis length", "cmean": "mean a-axis length",
"rhoiavg": "Mean ice density", "aspect": "Aspect ratio",
"vthabn": "Number weighted terminal fall speed",
"vthabm": "Mass weighted terminal velocity",
"dmdtvap": "Vapor growth rate", "ssi": "Ice supersaturation",
"drpp_dep": "Tendencies for depositional growth of pristine ice",
"drpp_sedim": "Tendencies for sedimentation of pristine ice",
"drrp_sedim": "Tendencies for sedimentation of rain",
"drrp_collect": "Tendencies for collection of cloud into rain",
"drt_horizAdv": "Tendencies in rt due to horizontal advection",
"drc_ncphys": "Tendencies in rc due to microphysics",
"drv_mcphys": "Tendencies in rv due to microphysics",
"drc_turb": "Tendencies in rc due to turbulence",
"drv_turb": "Tendencies in rv due to turbulence",
"dthil_subs": "Tendencies in thetail due to subsidence",
"drt_subs": "Tendencies in rt due to subsidence"}
if 'out_mc_scale' in kwargs.keys():
out_info["mc"]["scale"] = kwargs['out_mc_scale']
elif sim_version is True:
out_info["mc"]["scale"] = {"zm": "linear", "rt": "linear", "rvap": "linear",
"rcld": "log",
"rrain": "log", "rpice": "log", "rsnow": "log", "raggr": "log",
"rgrau": "log", "rhail": "log", "cpice": "log", "amean": "log",
"cmean": "log", "rhoiavg": "linear", "aspect": "linear",
"vthabn": "linear",
"vthabm": "linear", "dmdtvap": "linear" ,"ssi": "linear",
"zdr": "linear", "zh", "linear",
"drpp_dep": "symlog", "drpp_sedim": "log",
"drrp_sedim": "log", "drrp_collect": "log",
"drt_horizAdv": "symlog", "drc_ncphys": "symlog",
"drv_mcphys": "symlog", "drc_turb": "symlog",
"drv_turb": "symlog", "dthil_subs": "linear",
"drt_subs": "symlog"}
else:
out_info["mc"]["scale"] = {"zm": "linear", "rt": "linear", "rvap": "linear",
"rcld": "log",
"rrain": "log", "rpice": "log", "rsnow": "log", "raggr": "log",
"rgrau": "log", "rhail": "log", "cpice": "log", "amean": "log",
"cmean": "log", "rhoiavg": "linear", "aspect": "linear",
"vthabn": "linear",
"vthabm": "linear", "dmdtvap": "linear" ,"ssi": "linear",
"drpp_dep": "symlog", "drpp_sedim": "log",
"drrp_sedim": "log", "drrp_collect": "log",
"drt_horizAdv": "symlog", "drc_ncphys": "symlog",
"drv_mcphys": "symlog", "drc_turb": "symlog",
"drv_turb": "symlog", "dthil_subs": "linear",
"drt_subs": "symlog"}
# -------------------------------------init sfc------------------------------------
out_info["sfc"] = {}
if 'out_sfc_names' in kwargs.keys():
out_info["sfc"]["names"] = kwargs['out_sfc_names']
else:
out_info["sfc"]["names"] = ["time", "SHF", "LHF", "swd", "swu", "lwd", "lwu", "Qi", "SST"]
if 'out_sfc_units' in kwargs.keys():
out_info["sfc"]["units"] = kwargs['out_sfc_units']
else:
out_info["sfc"]["units"] = {"time": "s", "SHF": "W m^{-2}", "LHF": "W m^{-2}", "swd": "W m^{-2}",
"swu": "W m^{-2}", "lwd": "W m^{-2}", "lwu": "W m^{-2}", "Qi": "W m^{-2}", "SST": "K"}
if 'out_sfc_long_name' in kwargs.keys():
out_info["sfc"]["long_name"] = kwargs['out_sfc_long_name']
else:
out_info["sfc"]["long_name"] = {"time": "time", "SHF": "Sensible heat flux",
"LHF": "Latent heat flux", "swd": "Downwelling SW flux", "swu": "Upwelling SW flux",
"lwd": "Downwelling LW flux", "lwu": "Upwelling LW flux",
"Qi": "Heat flux from ocean into ice", "SST": "Surface temperature"}
if 'out_sfc_scale' in kwargs.keys():
out_info["sfc"]["scale"] = kwargs['out_sfc_scale']
else:
out_info["sfc"]["scale"] = {"time": "linear", "SHF": "linear", "LHF": "linear",
"swd": "linear",
"swu": "linear", "lwd": "linear", "lwu": "linear",
"Qi": "linear", "SST": "linear"}
return out_info