MeteoInfoLab脚本示例：计算湿位涡

良辰

MeteoInfo 发表于 2016-8-18 21:51
for j in range(nz):
    prs[:,j,:,:] = lev[j]

王老师想问一下，这个句话mpv[t,z-1,:,:] = pv.array想实现什么意思，如果我想单输出pv1，还需要写成这种格式吗？我修改了一段位涡（干位涡）的计算代码，就是假相当位温替换为位温，只输出pv1一项。但是不知道最后怎么给数组赋值！代码如下：不知道最后一句的写法！
# Calculate dry potential vorticity
print 'Calculate dry potential vorticity...'
pv = zeros([nt,nz,ny,nx], dtype='double')
pv = DimArray(pv, rh.dims)
pv.setdimvalue(1, lev[1:nz-1])
for t in range(nt):
    tt = meteof.gettime(t)
    print tt.strftime('%Y-%m-%d %H:00')
    for z in range(1, nz-1):
        f = zeros([ny,nx])
        f1 = 2*7.292*sin(lat*3.14159/180.0)*0.00001
        for i in range(nx):
            f[:,i] = f1
        g = 9.8
        dp = 100*(lev[z-1]-lev[z+1])
        dtheta = theta[t,z-1,:,:]-theta[t,z+1,:,:]
        du = uwnd[t,z-1,:,:]-uwnd[t,z+1,:,:]
        dv = vwnd[t,z-1,:,:]-vwnd[t,z+1,:,:]
        dx1 = 2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
        dx = zeros([ny,nx])
        for i in range(nx):
            dx[:,i] = dx1
        dy = 2.0*6370949.0*3.14159/180.0
        dtx = cdiff(eqt[t,z,:,:], 1)
        dty = cdiff(eqt[t,z,:,:], 0)
        pv = -g*(vort[t,z,:,:]+f)*dtheta/dp  
#        pv[t,z-1,:,:] = pv.array

已经解决！！！变量名重复了，做如下修改即可！
        pvd = -g*(vort[t,z,:,:]+f)*dtheta/dp  
        pv[t,z-1,:,:] = pvd.array

MeteoInfo

良辰 发表于 2016-8-18 22:17
王老师想问一下，这个句话mpv[t,z-1,:,:] = pv.array想实现什么意思，如果我想单输出pv1，还需要写成这 ...

嗯，这是个很明显的错误，另外.array去掉也是可以的。

良辰

MeteoInfo 发表于 2016-8-18 22:30
嗯，这是个很明显的错误，另外.array去掉也是可以的。

，学生愚钝！王老师，这个是我位涡追踪结果，最后一列！为什么有nan的出现？湿位涡也是一样出现nan,从第9个时次到25个时次！

MeteoInfo

良辰 发表于 2016-8-18 22:50
，学生愚钝！王老师，这个是我位涡追踪结果，最后一列！为什么有nan的出现？湿位涡 ...

可能是高度（气压）超出了数组高度维范围的缘故。

良辰

MeteoInfo 发表于 2016-8-18 23:03
可能是高度（气压）超出了数组高度维范围的缘故。

有解决办法吗？

MeteoInfo

良辰 发表于 2016-8-18 23:06
有解决办法吗？

可以考虑加一个判断，如果轨迹节点的气压值（p1）高于位涡数组最靠近地面的气压值（p2），可以将轨迹节点的气压值改为小于或等于p2。

良辰

好的，我试一下！早点休息王老师，非常感谢！

良辰

MeteoInfo 发表于 2016-8-18 23:12
可以考虑加一个判断，如果轨迹节点的气压值（p1）高于位涡数组最靠近地面的气压值（p2），可以将轨迹节点 ...

王老师，用grads画了一下湿位涡，和您的进行了对比，基本一致，用的fnl1*1的再分析资料。

MeteoInfo

良辰 发表于 2016-8-22 11:14
王老师，用grads画了一下湿位涡，和您的进行了对比，基本一致，用的fnl1*1的再分析资料。

不错，可以用fnl1*1数据在MeteoInfoLab中也做同样的计算，绘图时和GrADS用同样的色标再对比一下，会更清楚。

良辰

王老师，如果后向轨迹的时效比较长，比如后向10天的，一个数据资料集就不能满足计算需求，该如何添加那？学生这样meteof = addfile(meteofn1,meteofn2,meteofn3,meteofn4)，结果rh = meteof['RELH'][:,:,latlim,lonlim]报错，错误信息为TypeError: 'NoneType' object is unsubscriptable。王老师如有时间，希望指点。