标签:
需要温度和风场U/V分量格点数据,计算中主要用到cdiff函数,结果正确与否未经验证(希望有人用GrADS验证一下)。
脚本程序:
print ‘Open data files...‘ f_air = addfile(‘D:/Temp/nc/air.2011.nc‘) f_uwnd = addfile(‘D:/Temp/nc/uwnd.2011.nc‘) f_vwnd = addfile(‘D:/Temp/nc/vwnd.2011.nc‘) print ‘Read data array...‘ tidx = 173 # Jun 23, 2011 t = f_air.gettime(tidx) lidx = 1 # 1000 hPa air = f_air[‘air‘][tidx,lidx,::-1,:] uwnd = f_uwnd[‘uwnd‘][tidx,lidx,::-1,:] vwnd = f_vwnd[‘vwnd‘][tidx,lidx,::-1,:] lon = f_air[‘lon‘][:] lat = f_air[‘lat‘][:] lon, lat = meshgrid(lon, lat) # Calculate print ‘Calculate...‘ dtx = cdiff(air, ‘x‘) dty = cdiff(air, ‘y‘) dx = cdiff(lon, ‘x‘) * pi / 180 dy = cdiff(lat, ‘y‘) * pi / 180 tadv = -1*((uwnd*dtx)/(cos(lat*pi/180)*dx)+vwnd*dty/dy)/6.37e6 #Plot print ‘Plot...‘ axesm() mlayer = shaperead(‘D:/Temp/map/country1.shp‘) geoshow(mlayer, edgecolor=‘black‘) layer = contourfm(tadv, cmap=‘grads_rainbow‘) title(‘Temporature advection (‘ + t.strftime(‘%Y-%m-%d‘) + ‘)‘) colorbar(layer) xlim(0, 360) ylim(-90, 90)
标签:
原文地址:http://www.cnblogs.com/yaqiang/p/4641226.html