- 积分
- 3274
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2018-9-22
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
转换代码如下,大多数的数据都能正常运行,只有1988.4.17和1989.4.16这两天无法正常运行,对于1988.4.17这一天来说,转出来的02:00会显示成03:00,输出如下,HYSPLIT就会认为间隔是两小时,无法正常运行,但详细检查3d grib文件中的时间是正常的,想请教一下大家这是为什么呢?
Time index: 1 1988-04-17 01:00 Time index: 2 1988-04-17 03:00 Time index: 3 1988-04-17 03:00 Time index: 4 1988-04-17 04:00
# Convert ERA5 GRIB data to ARL data
import os
import datetime
#---- Set data folder
datadir = '/Volumes/MyPassport/ERA5/'
era5dir = '/Volumes/MyPassport/ERA5/'
for year in range(1988,1989):
begin = datetime.date(year,4,17)
end = datetime.date(year,4,17)
d = begin
delta = datetime.timedelta(days=1)
while d <= end:
cdate=d.strftime("%Y%m%d")
#print(cdate)
outfn = os.path.join(datadir, cdate+'.arl')
if os.path.exists(outfn):
d+=delta
continue
else:
print(cdate)
infn3d = addfile(era5dir+'ERA5_3d_'+cdate+'.grib')
infn2d = addfile(era5dir+'ERA5_2d_'+cdate+'.grib')
arlf = addfile(outfn, 'c', dtype='arl')
gvar2d = ['Surface_pressure_surface','2_metre_temperature_surface','10_metre_U_wind_component_surface',\
'10_metre_V_wind_component_surface','Boundary_layer_height_surface','Convective_available_potential_energy_surface',\
'Instantaneous_eastward_turbulent_surface_stress_surface','Instantaneous_northward_turbulent_surface_stress_surface']
gvar3d = ['Geopotential_isobaric','Temperature_isobaric','Vertical_velocity_isobaric',\
'U_component_of_wind_isobaric','V_component_of_wind_isobaric','Relative_humidity_isobaric']
#---- Corresponding variable names in ARL data file
avar2d = ['PRSS','T02M','U10M','V10M','PBLH','CAPE','UMOF','VMOF']
avar3d = ['HGTS','TEMP','WWND','UWND','VWND','RELH']
#--- Add DIFF fields - difference between the original data and the packed data
avar3d_diff = list(avar3d)
avar3d_diff.append('DIFW')
#---- Set parameters of ARL data file
gv = infn3d['Geopotential_isobaric']
nx = gv.dimlen(gv.ndim - 1)
ny = gv.dimlen(gv.ndim - 2)
levels = gv.dimvalue(gv.ndim - 3)[::-1]
print(levels)
nz = len(levels)
arlf.setlevels(levels)
arlf.set2dvar(avar2d)
for l in levels:
arlf.set3dvar(avar3d_diff)
arlf.setx(gv.dimvalue(gv.ndim - 1))
arlf.sety(gv.dimvalue(gv.ndim - 2))
tNum = infn3d.timenum()
fhour = 0
for t in range(0, tNum):
print 'Time index: ' + str(t)
atime = infn3d.gettime(t)
print atime.strftime('%Y-%m-%d %H:00')
dhead = arlf.getdatahead(infn3d.proj, 'RSMC', 2, fhour)
#Pre-write index record without checksum - will be over-write latter
arlf.writeindexrec(atime, dhead)
#Checksum list
ksumlist = []
# Write 2d variables
ksums = []
for avname,gvname in zip(avar2d, gvar2d):
gdata = infn2d[gvname][t,:,:]
if avname == 'PRSS':
gdata = gdata * 0.01
ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)
ksums.append(ksum)
ksumlist.append(ksums)
# Write 3d variables
for lidx in range(0, nz):
ksums = []
llidx = nz - lidx - 1
print(lidx,llidx)
for avname,gvname in zip(avar3d, gvar3d):
gdata = infn3d[gvname][t,llidx,:,:]
if avname == 'WWND':
gdata = gdata * 0.01
difw = arlf.diff_origin_pack(gdata)
elif avname == 'SPHU':
gdata = gdata * 1000.
elif avname == 'HGTS':
gdata = gdata / 9.80665
ksum = arlf.writedatarec(atime, lidx + 1, avname, fhour, 99, gdata)
ksums.append(ksum)
ksum = arlf.writedatarec(atime, lidx + 1, 'DIFW', fhour, 99, difw)
ksums.append(ksum)
ksumlist.append(ksums)
#Re-write index record with checksum
print(ksumlist)
arlf.writeindexrec(atime, dhead, ksumlist)
fhour += 1
arlf.close()
d+=delta
print 'Finished!'
|
|