-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcreate_LBC_HCLIM_grib.py
executable file
·305 lines (192 loc) · 8.49 KB
/
create_LBC_HCLIM_grib.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
import pygrib
import yaml
import numpy as np
import os
import datetime
from netCDF4 import Dataset
from tqdm import tqdm
from yaml.loader import SafeLoader
# ------------------------------------------------------------------------------
def is_year_leap(year):
"""
Check if a specific year is leap
Parameters:
year (int) : year of interest
Returns:
Logical : True (False) for leap (non-leap year)
"""
if (year % 400 == 0) and (year % 100 == 0):
return True
elif (year % 4 ==0) and (year % 100 != 0):
return True
else:
return False
def read_namelist(filename):
with open(filename) as f:
namelist = yaml.load(f, Loader=SafeLoader)
return namelist
def num_days_in_month(year, month):
daysinmonth = [29,31,28,31,30,31,30,31,31,30,31,30,31]
if is_year_leap(year) and month==2:
return daysinmonth[0]
else:
return daysinmonth[month]
def read_ERA5_ab(filename):
#Load constant data for vertical levels: ERA5 a & b on half-levels, PGW pressure levels
levm = np.loadtxt(filename)
alevm = levm[:,0]
blevm = levm[:,1]
return alevm, blevm
def read_delta_ab(namelist):
# Read sp deltas
nc_file = Dataset(namelist['dir_deltas']+'/PGW2/ensemble/delta_sp_PGW2_3h.nc','r')
am_delta = nc_file.variables['hyam'][:].data
bm_delta = nc_file.variables['hybm'][:].data
nc_file.close()
return am_delta, bm_delta
def read_delta_var(namelist,pgw,tstep,varname):
# Read sp deltas
nc_file = Dataset(namelist['dir_deltas']+'/'+pgw+'/ensemble/delta_'+varname+'_'+pgw+'_3h.nc','r')
if varname == 'sp':
delta_var = nc_file.variables['delta_'+varname][tstep,0,:,:].data
elif varname == 'skt' or varname == 'sst':
delta_var = nc_file.variables['delta_'+varname][tstep,:,:].data
else:
delta_var = nc_file.variables['delta_'+varname][tstep,:,:,:].data
nc_file.close()
return delta_var
def read_ablevels(namelist,nx,ny):
am_era5, bm_era5 = read_ERA5_ab(namelist['level_file'])
am_delta, bm_delta = read_delta_ab(namelist)
am_delta = np.expand_dims(am_delta,axis=(1,2))
am_delta = np.tile(am_delta, (1,nx,ny))
bm_delta = np.expand_dims(bm_delta,axis=(1,2))
bm_delta = np.tile(bm_delta, (1,nx,ny))
am_era5 = np.expand_dims(am_era5,axis=(1,2))
am_era5 = np.tile(am_era5, (1,nx,ny))
bm_era5 = np.expand_dims(bm_era5,axis=(1,2))
bm_era5 = np.tile(bm_era5, (1,nx,ny))
read_ab = False
return am_era5, bm_era5, am_delta, bm_delta, read_ab
def read_delta_vinterp(namelist, pgw, tstep, p_era5, p_delta, case):
delta_var = read_delta_var(namelist, pgw, tstep, case)
# Interpolate delta value
delta_var_interp = np.zeros(p_era5.shape) * np.nan
for jj in range(p_era5.shape[1]):
for ii in range(p_era5.shape[2]):
delta_var_interp[:,jj,ii] = np.interp(np.squeeze(p_era5[:,jj,ii]), p_delta[:,jj,ii], np.squeeze(delta_var[:,jj,ii]))
return delta_var_interp
def transform_q2rh(hus, ta, pa):
"""
Compute relative humidity from specific humidity.
"""
hur = 0.263 * pa * hus *(np.exp(17.67*(ta - 273.15)/(ta-29.65)))**(-1)
return(hur)
def transform_rh2q(hur, ta, pa):
"""
Compute specific humidity from relative humidity.
"""
hus = (hur * np.exp(17.67 * (ta - 273.15)/(ta - 29.65))) / (0.263 * pa)
return(hus)
def is_year_leap(year):
"""
Check if a specific year is leap
Parameters:
year (int) : year of interest
Returns:
Logical : True (False) for leap (non-leap year)
"""
if (year % 400 == 0) and (year % 100 == 0):
return True
elif (year % 4 ==0) and (year % 100 != 0):
return True
else:
return False
def calculate_DOY(year, month, day):
# Remove leap day, DOY same as 28th February
if is_year_leap(year) and month ==2 and day==29:
day_of_year = (datetime.datetime(1999, 2, 28) - datetime.datetime(1999, 1, 1)).days
else:
day_of_year = (datetime.datetime(1999, month, day) - datetime.datetime(1999, 1, 1)).days
return day_of_year
def calculate_tstep(year, month, day, hh):
doy = calculate_DOY(year,month,day)
return int(8*doy + hh/3)
# ------------------------------------------------------------------------------
# 1 - Read namelist
namelist = read_namelist('../nam/namelist_create_LBC_HCLIM_grib.yaml')
read_ab = True
# 2 - Loop over dates
for yy in namelist['years']:
for mm in [3]:#namelist['months']:
for dd in [1]:#range(1,num_days_in_month(yy,mm)+1):
for hh in [0]:#range(0,22,3):
tstep = calculate_tstep(yy,mm,dd,hh)
pgw1fname = namelist['input_dir']+'/'+str(yy)+'/ERA5_NEUR_{0:04d}{1:02d}{2:02d}{3:02d}00+000H00M'.format(yy,mm,dd,hh)
tempera5fname = namelist['template_dir']+namelist['template_file']
if os.path.isfile(pgw1fname):
# Open grib file
grbs_era = pygrib.open(pgw1fname)
grbs_tmp = pygrib.open(tempera5fname)
grb_var = grbs_era.message(703)
lnsp = grb_var.values[:,:]
sp = np.exp(lnsp)
# Read a/b levels
if read_ab: am_era5, bm_era5, am_delta, bm_delta, read_ab = read_ablevels(namelist, sp.shape[0],sp.shape[1])
for pgw in ['PGW2','PGW3']:
delta_sp = read_delta_var(namelist,pgw,tstep,'sp')
nsp = sp + delta_sp
# ERA5 pressure levels
p_era5 = am_era5 + bm_era5*nsp
# Delta pressure levels
p_delta = am_delta + bm_delta*nsp
nera5fname = namelist['output_dir']+'/'+str(yy)+'/ERA5_NEUR_'+pgw+'_{0:04d}{1:02d}{2:02d}{3:02d}00+000H00M'.format(yy,mm,dd,hh)
print(tstep, nera5fname)
grbs_nera = open(nera5fname,'wb')
#for igrb in tqdm(range(1,grbs_era.messages+1), desc='Variables'):
for igrb in range(1,grbs_era.messages+1):
grb_var = grbs_era.message(igrb)
grb_tmp = grbs_tmp.message(igrb)
if igrb == 3:
delta_sst = read_delta_var(namelist,pgw,tstep,'sst')
grb_tmp.values = grb_var.values + delta_sst
elif igrb in [11, 13, 15, 16, 17]:
delta_ts = read_delta_var(namelist,pgw,tstep,'skt')
grb_tmp.values = grb_var.values + delta_ts
elif igrb>=18 and igrb<=154:
if igrb == 18: delta_var_interp = read_delta_vinterp(namelist, pgw, tstep, p_era5, p_delta, 't')
if igrb == 18:
tem = np.zeros((p_era5.shape))
ntem = np.zeros((p_era5.shape))
ilevel = grb_var.level
grb_tmp.values = grb_var.values + delta_var_interp[ilevel-1,:,:]
tem[ilevel-1,:,:] = grb_var.values
ntem[ilevel-1,:,:] = grb_tmp.values
elif igrb>=155 and igrb<=291:
if igrb == 155: delta_var_interp = read_delta_vinterp(namelist, pgw, tstep, p_era5, p_delta, 'r')
ilevel = grb_var.level
rh = transform_q2rh(grb_var.values, tem[ilevel-1,:,:], am_era5[ilevel-1,:,:]+bm_era5[ilevel-1,:,:]*sp)
nrh = rh + delta_var_interp[ilevel-1,:,:]
nq = transform_rh2q(nrh, ntem[ilevel-1,:,:],p_era5[ilevel-1,:,:])
grb_tmp.values = nq
if igrb == 291:
del tem
del ntem
elif igrb>=292 and igrb<=428:
if igrb == 292: delta_var_interp = read_delta_vinterp(namelist, pgw, tstep, p_era5, p_delta, 'u')
ilevel = grb_var.level
grb_tmp.values = grb_var.values + delta_var_interp[ilevel-1,:,:]
elif igrb>=429 and igrb<=565:
if igrb == 429: delta_var_interp = read_delta_vinterp(namelist, pgw, tstep, p_era5, p_delta, 'v')
ilevel = grb_var.level
grb_tmp.values = grb_var.values + delta_var_interp[ilevel-1,:,:]
elif igrb == 703:
grb_tmp.values = np.log(nsp)
else:
grb_tmp.values = grb_var.values
grb_tmp.dataDate = yy*10000+mm*100+dd
msg = grb_tmp.tostring()
grbs_nera.write(msg)
grbs_nera.close()
grbs_era.close()
grbs_tmp.close()