又来问关于wrfout数据画温度平流的问题探讨了

supersupergirl

以下是我用两种方法画的温度平流 ，感觉都有些勉强，因为数值过大 且冷暖平流夹杂，想请大家帮忙看下哪种计算方法是正确的 或是有其他更好的方法求告知advect_variable函数法：

pressure = 700.
    tc   = wrf_user_getvar(a,"tc",it)      ; T in C
    rh   = wrf_user_getvar(a,"rh",it)      ; relative humidity
    u    = wrf_user_getvar(a,"ua",it)      ; 3D U at mass points
    v    = wrf_user_getvar(a,"va",it)      ; 3D V at mass points


  tc_p = tc+273.166
  copy_VarCoords(tc, tc_p)
  Tadv = advect_var(u, v, tc_p, 0, "advection of temperature", "K/s", 0)
  Tadv = -1*Tadv*100000
  copy_VarCoords(Tadv, u)
  Tadv_p = wrf_user_intrp3d(Tadv, p,"h",pressure,0.,False)


中央差分法：
  tk = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
  z_p = wrf_user_intrp3d(z,p,"h",pressure,0.,False)
  rh_p = wrf_user_intrp3d(rh,p,"h",pressure,0.,False)
  u_p = wrf_user_intrp3d(u,p,"h",pressure,0.,False)
  v_p = wrf_user_intrp3d(v,p,"h",pressure,0.,False)

  tc_p = tk+273.166
  copy_VarCoords(tk, tc_p)
  ;printVarSummary(tc_p)
  tc_p@_FillValue = 9.96921e+36
  dlon = (y(1,2)-y(1,1))*0.0174533
  dTdX = new((/141,153/),typeof(tc_p),tc_p@_FillValue)
  ;printVarSummary(dTdX)
  do nl = 0,141-1
  dX =6378388.*cos(0.0174533*x(nl,1))*dlon
  ;print(dX)
  ;printVarSummary(tc_p(nl:nl,:))
  dTdX(nl:nl,:) = center_finite_diff_n(tc_p(nl:nl,:),dX,True,0,1)
  end do
  dlat = (x(2,1)-x(1,1))*0.0174533
  dTdY = new((/141,153/),typeof(tc_p),tc_p@_FillValue)
  do ml = 0,153-1  
  dY = 6378388.*dlat
  ;print(dY)
  dTdY(:,ml:ml) = center_finite_diff_n(tc_p(:,ml:ml),dY,True,0,0)
  ;printVarSummary(dTdY)
  end do
  advT=-(u_p*dTdX/dX+v_p*dTdY/dY)*100000
  copy_VarCoords(u_p, advT)
  printVarSummary(advT)
  printMinMax(advT, 0)
  advT=smth9(advT,0.5,0.25,False)

不会游泳的fish

楼主你中央差分那里最后错了吧，应该u_p*DTdx就可以了不需要除以dx，因为DTDX已经把dx算进去了

supersupergirl

不会游泳的fish 发表于 2019-10-9 10:16
楼主你中央差分那里最后错了吧，应该u_p*DTdx就可以了不需要除以dx，因为DTDX已经把dx算进去了

不好意思现在才看到 是的 但我去掉dx还是不行 到现在也还没成功 求指教

freekiller

不知道你最后希望你的单位是什么。你应该明白wrf模式坐标， 格点上的经纬度是不相同的，为什么选dlon = (y(1,2)-y(1,1))*0.0174533就代表所有。 为什么不直接用格距算，还要去算两点之间的距离

supersupergirl

freekiller 发表于 2019-11-14 21:31
不知道你最后希望你的单位是什么。你应该明白wrf模式坐标， 格点上的经纬度是不相同的，为什么选dlon = (y( ...

不好意思现在才看到 请问直接用格距那不就都是一样的么 自己一开始转wrf设定的那个dx是吗

freekiller

supersupergirl 发表于 2019-11-15 12:18
不好意思现在才看到 请问直接用格距那不就都是一样的么 自己一开始转wrf设定的那个dx是吗

你去试试看看结果呗。 ncl中算绝对涡度用的就是dx dy。

supersupergirl

freekiller 发表于 2019-11-15 21:55
你去试试看看结果呗。 ncl中算绝对涡度用的就是dx dy。

你好 我试了之后发现和现在的图没有太大区别 还是很杂 请问是什么原因呢