- 积分
- 4282
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2021-3-2
- 最后登录
- 1970-1-1
![未绑定新浪微博用户 新浪微博达人勋](source/plugin/sina_login/img/gray.png)
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 是冉冉升起的冉 于 2024-6-6 11:05 编辑
结果跟摸鱼的对比了一下,但感觉在大西洋地区有一点差异,大体趋势差不多,20°低纬地区有一点奇怪,然后垂直方向没有对比;有问题联系我ing;- ;///////////////////////////////////////////////////////////////
- ; 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
复制代码
|
|