WRF格点下计算涡度散度谜团——操作篇

利用NCL为工具,比较了三种不同方法计算WRF输出风场对应的散度。

介绍

上一篇博客中介绍了WRF格点下计算涡度散度的一些思考,在这一篇博文里我们来看看实际操作中,使用NCL为已有的工具,对比三种思路一致但方法略有不同的操作所得到的结果有什么不同。

三种方法

NCL中对等经纬度格点(再分析资料常用的输出格式之一)的分析计算能力最强,计算散度有成熟的内置函数uv2dv_cfd,但对WRF输出结果的计算支持就相对要弱一些了。

由于uv2dv_cfd要求风场的经纬度是一维单调数组,我们自然想到可以利用NCL的插值函数rcm2rgrid先将WRF输出数据插值到rectilinear grid上,然后再使用NCL针对等经纬度格点的函数接口进行计算。其实WRF输出的任何其他量都可以采用类似的方法降低计算需要考虑的事项,这里我们希望对其准确度进行一个考量。

另一个方法便是采用上一篇博客介绍的公式,直接在WRF原生格点下利用正交曲线坐标系下的散度公式进行求导,并同时考虑Arakawa C交错格点以及地图放大系数,根据上一篇博客的论述,这样求散度时可以利用交错格点的优势使得求导精度高。

在这里还有一个问题——垂直坐标。由于WRF采用了地形跟随坐标系,其模式“水平面”并不平行于等高面或等压面,这一点在地形越为陡峭的地方越突出。而我们在诊断分析中所求的散度往往为针对等高面或等压面而言的水平散度,故模式面中的水平散度与等压面的水平散度往往并不相同。

由此,我们采用三种方法计算WRF输出风场对应的散度:
1.直接在模式“水平”面上采用正交曲线坐标系的公式计算散度,可以利用交错格点优势采用一倍格距精度;

2.先将模式风场插值到等压面上,再利用正交曲线坐标系的公式计算散度,由于NCL的垂直插值要求物理量必须先unstagger,故等压面上的风场已经在质量格点上了,无法利用交错格点优势,只能采用二倍格距精度计算散度;

3.先将模式风场(记得用旋转到经纬方向的风场)插值到等压面上,再将WRF格点数据插值到等经纬度坐标上(格距与模式原格距相当),最后直接采用NCL内置函数计算散度。

数据呈现

我们采用2005年1月10日模拟的一次青藏高原上空过程,wrfout的时间为12UTC,水平格距4公里,垂直层数51层,为了比较的可视化和公平性,我们将计算出的三种散度统一插值到$(86^oE,27.5^oN)–(88^oE,30.8^oN)$的垂直剖面上。

结果对比

compare

上图从左至右依次为方法1~3计算而来的散度(色标相同),散度正值为红色调,实等值线,负值为蓝色调,虚等值线。

可以看到由于地形追随坐标系在地形陡峭(青藏高原南坡)处的倾斜,第一种方法算得的水平散度在高原南坡与其他两种方法有很大的不同:第一种算法得散度在上坡处是强烈深厚的辐散气流,而其他两种先插值到等压面的算法则在南坡表现出了相对一致性,在高原南坡低层为辐合气流,至高层则转为辐散。从我们的需求上来看,等压面水平散度的计算在地形陡峭或者复杂的地区显然不能直接在模式“水平面”上进行再插值到等压面,这里的误差已经到了质变的程度。另外,在高原主体上空,地形相对不那么陡峭的地方,模式面与等压面的水平差异不那么大,可以看到此时的水平散度在pattern上就比较接近了。当然仔细看的话,在复杂小地形的近地面,第一种方法依旧会有相当程度的误差,到了高空误差减少。

对于第二种与第三种方法,可以看到无论是地形陡峭与否,总体的模态是十分相近的,特别是高原主体部分,定性地看两者可视为等同。有趣的是,主要的定量差异依旧在高原南坡。第二种方法在高原南坡为均匀的辐合区,而第三种方法在南坡却有很多小的辐散中心镶嵌于大的辐合区域之中,这里的原因并不明晰,从现在可以想到的理由看,可能是由于NCL从WRF格点到等经纬度双线性插值时引起的误差所造成的,这里还需要咱们聪明的读者来提示我原因。

从散度中心的强度上来看,不出所料,第一种方法由于利用了交错格点的优势,导数格距小,精度高,所对应的中心强度也是最强的,所求得场棱角分明;第二种方法则由于NCL垂直插值的限制,使用了两倍格距,其中心强度则较第一种弱一些,也较平滑一些;第三种方法由于还经过了格点的再插值,然后又使用了中央差分进行散度计算,相当于进行了两次平滑,故平滑程度程度最高,对应的中心强度也最弱。这也是为何高精度模拟和再分析资料所画的图会比较棱角分明,而我们用中央差分求得的场会比较平滑,这是与导数的空间精度密不可分的。即导数的量值大小与导数的空间精度、所用资料的格距直接相关,不同格距、精度的导出量没有可比性。

一些结论

1.在地形陡峭与复杂的地区,必须先将风场插值到等压面,再求等压面上的诊断量;

2.将WRF格点插值到等经纬度格点上,再用标准方法求一些诊断量,从定性、模态观察、作图角度是可以接受的,但要将经纬度格距调小一些;

3.从定量角度,在WRF原生格点下计算可以保留更多模式原信息,提高计算精度;在地形陡峭处也似乎可以得到更为准确的值;

4.格距越小、求导空间精度越高,所求的导数数值越大、越不平滑(本质上其包含的尺度、描述的现象越小);

5.如果所求量需要进一步代入其他等式进行计算,务必要确保所有参与计算的量具有相同的空间精度,否则没有可比性。

一些途中遇到的NCL题外话

函数wrf_user_intrp3d不允许被插值的函数有缺测(不会做缺测值检查与处理),故如果用自己算出来的值再去用这个函数插值时要小心(用一下where保证没有奇异值)。

NCL的array的coordinate subscripts并不做插值,也就是说,如果coordinate是非整数,而你的subscript是整数,给出的值并不会将其插值到整数坐标位置,而是找一个已有坐标里离这个整数最近且在subscript范围内(如果是{a:b}的话),或者是直接绝对值最近的坐标(如果是单单一个{a}的话),所以数据对应不同坐标时你看上去用了同样的subscript,数据给的值所在位置也许并不相同,上面三种方法,横坐标的实际位置其实是不一样的,包括山体(缺测值)的形状不同,这都与subscript不会插值、不同方法所用空间格点不同有关。

关于为何第二种方法地形在高原平台上比较准而在南坡缺测过多的原因:可以看到在高原的形状(数据缺测的形状位置)上,第二种方法在高原主体比较准确,但是在南坡缺测值过多,第三种方法则正好相反。直接在模式面上算(第一种方法)一定是缺测最少的(原数据就没有缺测);后两种方法,插值到等压面大家是一样的,第二种方法最后一步是将凡是用到的u,v有缺测的div就赋予缺测(Fortran中),而我怀疑uv2dv_cfd是尽量保留不要缺测,所以它在地形交界处的地形比用第二种方法细节更多。


关于在NCL中使用Fortran90子函数

在NCL中使用Fortran90子函数(procedure与function均可)其实很简单,Fortran中按照往常习惯编写就好,不用加什么特别的注释。然后需要写一个stub文件,把Fortran子函数的接口和传入传出的变量声明部分写好,头尾加上C NCLFORSTART与C NCLEND,然后使用WRAPIT命令行做一个.so的库文件

WRAPIT stub文件 Fortran90文件

这样便可以在NCL中使用这个Fortran90子函数了。只需要在NCL脚本开头用external给这个库文件起个别名,用的时候采用

别名::Fortran90函数名

就可以了!

和f2py一样,NCL(Python)中数组最右侧(row-major)循环最快,Fortran数组则是最左侧循环最快(column-major),但是在编写的时候不用考虑这个问题,Fortran中将最左侧理解为经度,NCL中将最右侧理解为经度,WRAPIT会自动给你互换这些数字在数组中的对应位置。

×

碰并增长

扫码支持
扫码打赏,你说多少就多少

打开支付宝扫一扫,即可进行扫码打赏哦

文章目录
  1. 1. 介绍
  2. 2. 三种方法
  3. 3. 数据呈现
  4. 4. 结果对比
  5. 5. 一些结论
  6. 6. 一些途中遇到的NCL题外话
  7. 7. 关于在NCL中使用Fortran90子函数
,