WRF格点计算散度和涡度

happymgj

自己突发奇想，尝试求取了wrfout数据的涡度散度，也请大家指点指正一下。
wrfout数据在水平方向上是等距格点，而非等经纬度格点。这种等距格点也是“正交曲线坐标系”。
根据场论，正交曲线坐标系计算散度、涡度时，应当考虑拉梅系数（Lamé coefficient）。
首先看正交曲线坐标系的散度：

其中，H1，H2为拉梅系数。

而拉梅系数H与地图放大系数m互为倒数。

则也可以写为

其中，m1和m2分别为“x”与“y”方向的地图放大系数（拉梅系数倒数）。

注意“x”和“y”，并非是局地直角坐标系的经向和纬向，而是正交曲线坐标系下的坐标方向。

同样，涡度也有如下写法：

由此，编写出相应求涡度、散度程序的代码

    u00 = wrf_user_getvar(f,"ua",0)      ;mass格点的u速度，3d变量

    v00 = wrf_user_getvar(f,"va",0)       ;mass格点的v速度，3d变量

    dx=30000.                    

    dy=30000.

    m_x0=f->MAPFAC_MX(0,:,:)       ;mass格点的“x”方向的地图放大系数，但是其为2d变量

    m_y0=f->MAPFAC_MY(0,:,:)       ;mass格点的“y”方向的地图放大系数，但是其为2d变量

    m_x=conform(u00, m_x0, (/1,2/))  ;把“x”方向地图放大系数维度变成与u00 一致

    m_y=conform(u00, m_y0, (/1,2/))  ;把“y”方向地图放大系数维度变成与u00 一致

   div_mx=center_finite_diff_n(u00/m_y, dx, False, 1, 2)  ;求取"du/dx”

   div_my=center_finite_diff_n(v00/m_x, dy, False, 1, 1)  ;求取"du/dy”

   div_m=m_x*m_y*(div_mx+div_my)*1.e5

   vort_mx=m_x*center_finite_diff_n(v00/m_y, dx, False, 1, 2);求取"dv/dx”

   vort_my=m_y*center_finite_diff_n(u00/m_x, dy, False, 1, 1)；求取"dv/dy”

   vort_m=(m_y*vort_mx-m_x*vort_my)*1.e5

变量下标带了m，就是考虑了“x”和“y”方向的地图放大因子（拉梅系数的倒数）。

对此也编写了程序进行对比

右图为：用上述公式计算的相对涡度

左图为：用rcm2rgrid插值成等经纬度格点计算的相对涡度

受rcm2rgrid插值经度及绘图前数据平滑影响，其误差应该是可以接受的。

freekiller

{:eb502:}{:eb502:}{:eb502:},之前ncl论坛有个fortran程序，用ncl调用计算涡度，散度，考虑了地图放大因子

smy@202010

哇，学习了学习了，感谢楼主！

Ivy

太棒了，感谢楼主提供的计算方法，很有用！

HD_chan

我看了wrf公式中绝对涡度的计算，还考虑了质点格点的地图放大系数和科式力，不知道在相对涡度中是否可以忽略？