计算波活动通量出错

隽今、

救救孩子，如果顺利毕业了一定记得您的大恩大德！

能帮我看看计算波活动通量的代码有哪里出错了吗？参考了气象家园某位贴主的计算，而且比对了公式感觉没有问题，但是老师说图画的不对，到底是为什么哇

Omega= systemfunc("ls /mnt/d/ncl/aaa/pressure/omega/*")  

Uwnd = systemfunc("ls /mnt/d/ncl/aaa/pressure/uwnd/*")

Vwnd = systemfunc("ls /mnt/d/ncl/aaa/pressure/vwnd/*")

Hgt  = systemfunc("ls /mnt/d/ncl/aaa/pressure/hgt/*")

Air  = systemfunc("ls /mnt/d/ncl/aaa/pressure/air/*")

do i = 44,44

   ;;读取涡度

   Ome  = addfile(Omega(i), "r")

   omega= Ome->omega(273:313,:,:,:)

   ;;读取纬向风

   Uw   = addfile(Uwnd(i), "r")

   u    = Uw->uwnd(273:313,:,:,:)

   ;;读取经向风

   Vw   = addfile(Vwnd(i), "r")

   v    = Vw->vwnd(273:313,:,:,:)

   ;;读取位势高度

   H    = addfile(Hgt(i), "r")

   hgt  = H->hgt(273:313,:,:,:)

   ;;读取温度

   A    = addfile(Air(i), "r")

   lat  = A->lat

   lon  = A->lon

   level= A->level  

   time = A->time(273:313)

   t    = A->air({time},:,:,:)

   ;;———————————计算三维波活动通量———————————

   ntime  = dimsizes(time)

   nlat   = dimsizes(lat)

   nlon   = dimsizes(lon)

   nlevel = dimsizes(level)

   ;;基本量

   ; Gas constant

   gc     = 290

   ; Gravitational acceleration

   ga     = 9.80665

   ; Radius of the earth

   re     = 6378388;6.37*10^6

   ; scale height

   sclhgt = 8000.

   k      = 0.286

   ; pi

   pi     = atan(1.0)*4.

   ; Coriolis parameter

   f      = 2.*2.*pi/(60.*60.*24.)*sin(pi/180.*lat(:))

   f!0    = "lat"

   f&lat  = lat

   f@_FillValue = hgt@_FillValue

   ;Rotational angular velocity of the earth

   wmg     = 2.*pi/(60.*60.*24.)

   z       = hgt*ga

   copy_VarMeta(hgt,z)

   ; cosine

   coslat  = cos(lat(:)*pi/180.)

   coslat@_FillValue = default_fillvalue(typeof(coslat))

   coslat  = (/where(coslat.eq.0.,coslat@_FillValue,coslat)/)

   ; sin2lat

   sin2lat = sin(2*(lat(:)*pi/180.))

   sin2lat@_FillValue= default_fillvalue(typeof(sin2lat))

   sin2lat = (/where(sin2lat.eq.0.,sin2lat@_FillValue,sin2lat)/)

   ; 1-D -> 4-D

   leveltmp   = conform_dims(dimsizes(hgt),level,1)

   ftmp       = conform_dims(dimsizes(hgt),f,2)

   coslattmp  = conform_dims(dimsizes(hgt),coslat,2)

   sin2lattmp = conform_dims(dimsizes(hgt),sin2lat,2)

   ;Zonal deviation

   t_ano   = dim_rmvmean_n_Wrap(t, 3)

   u_ano   = dim_rmvmean_n_Wrap(u, 3)

   v_ano   = dim_rmvmean_n_Wrap(v, 3)

   z_ano   = dim_rmvmean_n_Wrap(z, 3)

   ;Zonal mean

   t_mean    = dim_avg_n_Wrap(t, 3)

   dt_meandz = center_finite_diff_n(t_mean,-sclhgt*log(level/1000.),False,0,1)

   ;

   S = dt_meandz+k*t_mean/sclhgt

   copy_VarCoords_1(hgt, S)

   ; 3-D -> 4-D

   Stmp = conform_dims(dimsizes(hgt),S,(/0,1,2/))

   ; dvz/dlon

   dvzdlon      =  center_finite_diff_n(v_ano*z_ano,lon*pi/180.,True,0,3)

   ; duz/dlon

   duzdlon      =  center_finite_diff_n(u_ano*z_ano,lon*pi/180.,True,0,3)

   ; dtz/dlon

   dtzdlon      =  center_finite_diff_n(t_ano*z_ano,lon*pi/180.,True,0,3)

   ; x-component

   Fx      = leveltmp/1000.*coslattmp*(v_ano*v_ano-1./(2.*wmg*re*sin2lattmp)*dvzdlon)

   ;printVarSummary(Fx)

   ; y-component

   Fy      = leveltmp/1000.*coslattmp*(-u_ano*v_ano+1./(2.*wmg*re*sin2lattmp)*duzdlon)

   ; z-component

   Fz      = leveltmp/1000.*coslattmp*(ftmp/Stmp*(v_ano*t_ano-1./(2.*wmg*re*sin2lattmp)*dtzdlon))

   copy_VarCoords(hgt,Fx)

   copy_VarCoords(Fx,Fy)

   copy_VarCoords(Fx,Fz)

   Fx@units    = "m^2/s^2"

   Fx@var_desc = "wave-activity flux"

   Fx@long_name= "x-component of wave-activity flux (Plumb1985)"

   Fy@units    = "m^2/s^2"

   Fy@var_desc = "wave-activity flux"

   Fy@long_name= "y-component of wave-activity flux (Plumb1985)"

   Fz@units    = "Pa m/s^2"

   Fz@var_desc = "wave-activity flux"

   Fz@long_name= "z-component of wave-activity flux (Plumb1985)"

end do

stepdance

由于是非线性量，波活动量严格来说需要用日数据计算而非月数据。如果使用月数据做回归或合成分析，则对顺序有要求：先算异常，再算活动量（但也只是减小了非线性量的误差而非完全消除）

隽今、

stepdance 发表于 2025-3-18 15:53
由于是非线性量，波活动量严格来说需要用日数据计算而非月数据。如果使用月数据做回归或合成分析，则对顺序 ...

谢谢您的留言，这里使用的是日数据，如果是日数据的话有什么顺序要求吗

stepdance

隽今、 发表于 2025-3-18 16:04
谢谢您的留言，这里使用的是日数据，如果是日数据的话有什么顺序要求吗

一般不会有要求，除非你分析的对象具有明显的日变化

是冉冉升起的冉

波活动通量画的是850 hPa的么？很少见画这个高度的通量，换250之类的看看会不会有什么不一样呢？输入数据和计算过程没有问题的话，也可能结果就是这样啦，老师觉得不对的地方在哪里呢？plumb波纬向分量确实会大一些，可以试试TN通量

隽今、

是冉冉升起的冉 发表于 2025-3-19 20:26
波活动通量画的是850 hPa的么？很少见画这个高度的通量，换250之类的看看会不会有什么不一样呢？输入数据和 ...

200和500hPa的图和850hPa的图类似，我看您画的波活动通量有各个方向的，可能我画的都是向东的就不太对。

是冉冉升起的冉

隽今、 发表于 2025-3-19 22:26
200和500hPa的图和850hPa的图类似，我看您画的波活动通量有各个方向的，可能我画的都是向东的就不太对。

对于plumb波这是正常的，因为纬向分量比较大，我建议你画tn看看。

隽今、

是冉冉升起的冉 发表于 2025-3-20 09:38
对于plumb波这是正常的，因为纬向分量比较大，我建议你画tn看看。

好的谢谢您