使用p_interp.F90程序对wrfout进行气压层插值的一些发现

pandasp1213

关键字：wrfout； p_interp.F90； 气压层插值； Fortran

WRF模式的输出结果wrfout是在混合坐标层上的，不是标准气压层，因此在实际的分析应用时需要对wrfout进行处理，把变量插值到某些固定气压层。

这一问题可以用很多软件解决，之前一直用的ncl，wrf_user_intrp3d函数，最近，因为一些条件限制，需要用Fortran处理。看到WRF手册上提到了p_interp.F90这个程序，在论坛上搜了下，发现用这个程序的人貌似并不多，因为很多软件都可以方便的实现气压层插值。于是自己开始琢磨，下面是此过程中的一些发现。

在WRF官网上找到了p_interp.F90这个程序，下载地址：http://www2.mmm.ucar.edu/wrf/src/。里面有几个不同的P_INTERP压缩包，先下载了2011年的P_INTERP_serial.TAR.gz，开始插值。

先拿500hPa的位势高度场做试验，在namelist.pinterp中设置：process = 'list'；fields = 'GHT'；interp_levels = 850. 。运行程序，出现错误，程序卡着不动了，CTRL+C也切不断，直接关窗口。看wrfout***_PLEV的大小只有2.6K。于是把fields改成默认的fields = 'PRES,TT,GHT,HGT'，再运行，结果一切正常。于是开始对比这两种情况，一步步检查程序。

最后终于找到了错误原因，是由于定义输出文件wrfout***_PLEV的维度和变量时出错。在上面两种情况下，定义输出文件时略有不同，程序在判断输出文件的维度和变量时是分两步来的，第一步循环wrfout自带变量，如果某一变量需要输出，则在wrfout***_PLEV中定义该变量，以及对应的维度，第二步，判断诊断变量，同理添加诊断变量。对于第二种情况，fields中的HGT属于wrfout自带变量，因此在循环完自带变量时，定义的wrfout***_PLEV中已经有时间、层、纬度、经度4个维度，第二步把诊断变量添加进去，一切正常。然而，第一种情况，GHT不是wrfout自带变量，因此第一步循环wrfout自带变量之后，wrfout***_PLEV中的维度只有时间、层，维度和经度是没有的，在第二步时，判断出需要添加诊断变量GHT，在往wrfout***_PLEV中添加新的维度（维度、经度）和变量（GHT）时，程序bug出现了，维度和经度的维度名称搞错了，并且没有打开nf_def_dim，直接导致程序傻掉了。

找到原因就好解决了，顺利解决之后，就准备整理整理思路。突然想到，WRF官网上放了好几个压缩包，不会都是有bug的吧，于是又下了支持并行的P_INTERP.TAR.gz，2010年发布的，设置process = 'list'；fields = 'GHT'，运行成功了！这时看看源程序才发现，并行版的程序和串行版的不一样，判断输出变量时不是分两步来的，这样就不会出现上述问题。至此，前些天都白忙活了。直接用并行版本就行了，并且并行版本还有一些新的实用功能。

时间花了，不能白花，至少留下些记录。

在检查程序时参考的两个网站：http://www.unidata.ucar.edu/soft ... newdocs/netcdf-f77/，查看nf函数，还有
http://www.unidata.ucar.edu/soft ... nc-error-codes.html，查看nf函数返回的status值代表含义。
另外需要注意的，p_interp.F90中的重力加速度g取的9.81，不是9.8或者9.8065，在采用其他软件处理时需要注意g的数值可能不同，出来的位势高度场会不同（图上可以明显反应出来的）。


在编译时需要加上netcdff库，论坛中有人问到过这个问题。

jstzysq

赞一个，学习啦

pandasp1213

jstzysq 发表于 2017-2-2 15:42
赞一个，学习啦

谢谢 谢谢 {:5_257:}{:5_257:}

xyan88

学习了学习了学习了

llw

学习了学习了学习了

竹竿民民

谢谢楼主，学习！

穆斯塔法

赞赞赞 学习了

古魏小

感谢楼主！