😀开心，今天终于配置好了评论！！！本来在微信公众号就是为了评论！结果没注意现在改版不能评论了！就很烦！所以继续同步github的内容！
因为老师说TN通量的向下不一定是波活动通量向下传播！可能是向上的减少，其异常是减去了气候带的结果，Kodera这篇文章中也提到了这个；
所以现在计算了一下Plumb的结果，plumb是纬向偏差，包含了时间的总场。但是不知道这两者算出来差异会不会很大！按理说感觉应该不会太大吧！
2008-Kodera., et al.,-GRL-Tropospheric impact of reflected planetary waves from the stratosphere.
结果和摸鱼的对比了一下，大体趋势差不多，但是垂直方向没有验证，有问题联系我。
;///////////////////////////////////////////////////////////////
; Calculate 3-D wave-activity flux:plumb[1985]
;///////////////////////////////////////////////////////////////
; Written by Ran,20240606.
; E-mail: 602023280016@nju.edu.cn
;If there is a problem with the code calculation, please contact me by email.
;Input data：hgt t u v
;Note that the data downloaded by ERA5 is the gravitational potential (unites:m^2/s^2).
;The ncep data is the gravity potential height z, which needs to be multiplied by g.
;Input data must include the entire longitude！！！Because to take the zonal deviation!!
;Return value - Four multi-dimentional arrays containing:
; x-component: Fx
; y-component: Fy
; z-component: Fz
; z_ano
;Usage: [/Fx,Fy,Fz,z_ano/] = plumb_pLLL(hgt,t,u,v)
undef("plumb_pLLL")
function plumb_pLLL(hgt[*][*][*][*]:numeric,t[*][*][*][*]:numeric,u[*][*][*][*]:numeric,v[*][*][*][*]:numeric)
local re,sclhgt,k,pi,dName,level,lat,lon,f,wmg,ga,z,coslat,sin2lat,leveltmp,ftmp,coslattmp,sin2lattmp,t_ano,u_ano,v_ano,t_mean,dt_meandz,\
S,Stmp,dvzdlon,duzdlon,dtzdlon
begin
;基本量
; Radius of the earth
re = 6.37*10^6;6378388
; scale height
sclhgt = 8000.
k=0.286
; pi
pi = atan(1.0)*4.
printVarSummary(hgt)
dName = getvardims(hgt)
if (any(ismissing(dName(1:)))) then
print("fatal: plumb_pLLL: requires that all the rightmost dimensions be named")
exit
end if
level = hgt&$dName(1)$
lat = hgt&$dName(2)$
lon = hgt&$dName(3)$
; Coriolis parameter
f = lat(:)
f = (/2.*2.*pi/(60.*60.*24.)*sin(pi/180.*f)/)
f@_FillValue = hgt@_FillValue
;Rotational angular velocity of the earth
wmg=2.*pi/(60.*60.*24.)
; Gravitational acceleration
ga = 9.80665
z=hgt*ga
copy_VarMeta(hgt, z)
; cosine
coslat = cos(lat(:)*pi/180.)
coslat@_FillValue = default_fillvalue(typeof(coslat))
coslat = (/where(coslat.le.0.,coslat@_FillValue,coslat)/)
; sin2lat
sin2lat = sin(2*(lat(:)*pi/180.))
sin2lat@_FillValue = default_fillvalue(typeof(sin2lat))
sin2lat = (/where(sin2lat.le.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)
;printVarSummary(t_ano)
;printVarSummary(z_ano)
;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)
;printVarSummary(S)
; 3-D -> 4-D
Stmp = conform_dims(dimsizes(hgt),S,(/0,1,2/))
;printVarSummary(Stmp)
; 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)
;printVarSummary(sin2lattmp)
;printVarSummary(dvzdlon)
;printVarSummary(wmg)
;printVarSummary(coslattmp)
; 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)"
return [/Fx,Fy,Fz,z_ano/]
end
