EP通量和TN波通量
EP通量
局地 E-P 通量（Trenberth, 1986）能有效诊断瞬变波对时间平均流的作用。局地 E-P 通量的计算如下
其中上标“－”和“′”分别表示时间平均（根据研究需要来设定时间平均范围）和时间扰动。
Plumb 波作用通量（Plumb, 1985）在球坐标下三维波活动通量的诊断方程为
其中上标“－”和“′”分别表示纬圈平均和纬向偏差。
(1) 式和 (2) 式中， 表示纬度、经度、位势，， 分别表示纬度、经度、位势，科氏参数，地球半径和地球自转速率。, , , 其中 是气体常数与定压比热之比， 是标高。（1）中 下标 表示偏导。（2）中 、 是地转风。
TN波通量
T-N 波作用通量 (Takaya and Nakamura, 2001) 的三维表达式如下：
其中：
是准地转流函数相对于气候场的扰动，基本流场 表示气候场，其他参数同 (2) 式。NCL代码
Takaya 和 Nakamura (2001, JAS: TN01)推导出的3D波活动通量，这里使用月平均数据进行计算。NCL代码如下：
;-------------------------------------------------------
; 3-D wave-activity flux derived by Takaya and Nakamura (1999, 2001)
; See equation (38) in Takaya and Nakamura (2001, JAS)
;-------------------------------------------------------
;
; Used data:
; - Monthly-mean data of NCEP/NCAR reanalysis 1
; - Geopotential height (hgt: m) : hgt.mon.mean.nc
;
; - Monthly climatology of NCEP/NCAR reanalysis 1
; - Geopotential height (hgt: m) : hgt.mon.mean.nc
; - Air temperature (air: degC) : air.mon.mean.nc
; - Zonal wind (uwnd: m/s) : uwnd.mon.mean.nc
; - Meridional wind (vwnd: m/s) : vwnd.mon.mean.nc
;
; Data source:
; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
;
; Calculation period: January 1990 to December 1990
; (To modify, change fyear, fmon, lyear, lmon)
;
; Pressure unit: [hPa]
; Basic state: Monthly climatology
; Perturbation: Deviation from climatology
;
; Outputs:
; x-component: TN2001-Fx.monthly.1990.nc
; y-component: TN2001-Fy.monthly.1990.nc
; z-component: TN2001-Fz.monthly.1990.nc
; QG stream function anomaly: psidev.monthly.1990.nc
; Brunt-Vaisala frequency: NN.monthly.1990.nc
;-------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; The first and last dates of calculation
fyear = 1990
fmon = 1
lyear = 1990
lmon = 12
; Directory paths
diri = "/e3b/ncep/plev/monthly/nc/"
dirci = "/e3b/ncep/plev/monthly-climate/nc/"
; Monthly mean data
filename = systemfunc("ls "+diri+"hgt.mon.mean.nc")
zfile = addfile(filename, "r")
; Monthly climatology data
filename = systemfunc("ls "+dirci+"air.mon.ltm.nc")
btfile = addfile(filename, "r")
filename = systemfunc("ls "+dirci+"hgt.mon.ltm.nc")
bzfile = addfile(filename, "r")
filename = systemfunc("ls "+dirci+"uwnd.mon.ltm.1968-1996.nc")
bufile = addfile(filename, "r")
filename = systemfunc("ls "+dirci+"vwnd.mon.ltm.1968-1996.nc")
bvfile = addfile(filename, "r")
; Test variable type and read data accordingly
if (getfilevartypes(zfile, "hgt").eq."short") then
zvarorg = short2flt(zfile->hgt)
else
zvarorg = zfile->hgt
end if
if (getfilevartypes(btfile, "air").eq."short") then
btvar = short2flt(btfile->air) + 273.15 ; Convert air to Kelvin
else
btvar = btfile->air + 273.15
end if
; Similarly read other variables (bzvar, buvar, bvvar)
if (getfilevartypes(bzfile, "hgt").eq."short") then
bzvar = short2flt(bzfile->hgt)
else
bzvar = bzfile->hgt
end if
if (getfilevartypes(bufile, "uwnd").eq."short") then
buvar = short2flt(bufile->uwnd)
else
buvar = bufile->uwnd
end if
if (getfilevartypes(bvfile, "vwnd").eq."short") then
bvvar = short2flt(bvfile->vwnd)
else
bvvar = bvfile->vwnd
end if
; Extract dimensions
time = zfile->time
lat = zfile->lat
lon = zfile->lon
level = zfile->level
ntime = dimsizes(time)
nlat = dimsizes(lat)
nlon = dimsizes(lon)
nlevel = dimsizes(level)
; Calendar setup for input file
time@calendar = "standard"
option = 0
option@calendar = time@calendar
utc_date = cd_calendar(time, option)
syear = tointeger(utc_date(0,0))
smon = tointeger(utc_date(0,1))
; Step calculations for range of data to process
fstep = (fyear - syear)*12 + fmon - smon
lstep = (lyear - syear)*12 + lmon - smon
; Extract the height anomaly
zvar = zvarorg(fstep:lstep,:,:,:)
; Allocate arrays for climatology
czvar = new((/ntime,nlevel,nlat,nlon/), float, zvarorg@_FillValue)
ctvar = new((/ntime,nlevel,nlat,nlon/), float, zvarorg@_FillValue)
cuvar = new((/ntime,nlevel,nlat,nlon/), float, zvarorg@_FillValue)
cvvar = new((/ntime,nlevel,nlat,nlon/), float, zvarorg@_FillValue)
; Loop through each step to compute values
do istep = 0, ntime-1
iyear = tointeger(utc_date(fstep + istep,0))
imon = tointeger(utc_date(fstep + istep,1))
; Find corresponding climatological step
jstep = imon - bsmon
if (jstep.lt.0) then
jstep = jstep + 12
end if
czvar(istep,:,:,:) = bzvar(jstep,:,:,:)
ctvar(istep,:,:,:) = btvar(jstep,:,:,:)
cuvar(istep,:,:,:) = buvar(jstep,:,:,:)
cvvar(istep,:,:,:) = bvvar(jstep,:,:,:)
end do
; Calculate height anomaly and other relevant variables
zavar = zvar - czvar
delete(czvar)
; Constants
gc = 290 ; Gas constant
ga = 9.80665 ; Gravitational acceleration
re = 6378388 ; Earth radius
sclhgt = 8000. ; Scale height
pi = 4.*atan(1.) ; Pi value
; Calculate Coriolis parameter and other necessary variables
f = 2. * 2. * pi / (60. * 60. * 24.) * sin(pi/180. * lat(:))
f!0 = "lat"
f&lat = lat
f@_FillValue = zvarorg@_FillValue
; Other variable calculations, including stream function, gradient computations...
; Output results
ncFx = addfile("TN2001-Fx.monthly.1990.nc", "c")
ncFy = addfile("TN2001-Fy.monthly.1990.nc", "c")
ncFz = addfile("TN2001-Fz.monthly.1990.nc", "c")
ncpsidev = addfile("psidev.monthly.1990.nc", "c")
ncNN = addfile("NN.monthly.1990.nc", "c")
; Save outputs
ncFx->Fx = Fx
ncFy->Fy = Fy
ncFz->Fz = Fz
ncpsidev->psidev = psidev
ncNN->NN = NN
end
Python代码
import numpy as np
import xarray as xr
### Constants
gc = 290.0 # Gas constant
g = 9.80665 # Gravitational acceleration
re = 6378388.0 # Earth radius
sclhgt = 8000.0 # Atmospheric scale height
omega = 7.292e-5 # Earth's angular velocity
def TN_Wave_Activity_Flux(T_Climo: xr.DataArray,
U_Climo: xr.DataArray,
V_Climo: xr.DataArray,
Geo_Climo: xr.DataArray,
Geo_ano: xr.DataArray) -> xr.DataArray:
"""
Compute T-N flux. The input is 3D DataArray data (level, lat, lon).
Output are flux components.
Parameters:
----------
T_Climo : xarray.DataArray (level, lat, lon)
Climatological temperature field in Kelvin.
U_Climo : xarray.DataArray (level, lat, lon)
Climatological zonal wind (east-west component) in m/s.
V_Climo : xarray.DataArray (level, lat, lon)
Climatological meridional wind (north-south component) in m/s.
Geo_Climo : xarray.DataArray (level, lat, lon)
Climatological geopotential.
Geo_ano : xarray.DataArray (level, lat, lon)
Geopotential anomaly.
"""
# Get the data coordinates
data_coords = T_Climo.coords
lon, lat, level = T_Climo.lat, T_Climo.lon, T_Climo.level
# Calculate total wind speed
UV_Climo = np.sqrt(U_Climo ** 2 + V_Climo ** 2)
# Get longitude, latitude, and pressure level data
lon = np.deg2rad(lon)[np.newaxis, np.newaxis, :]
lat = np.deg2rad(lat)[np.newaxis, :, np.newaxis]
level = np.array(level)[:, np.newaxis, np.newaxis]
# Convert DataArrays to numpy arrays
T_Climo, U_Climo, V_Climo, UV_Climo, Geo_Climo, Geo_ano = map(np.array, [T_Climo, U_Climo, V_Climo, UV_Climo, Geo_Climo, Geo_ano])
# Retain only eastward wind, set westward wind to NaN
positive_flow_mask = U_Climo >= 0
T_Climo, U_Climo, V_Climo, UV_Climo, Geo_Climo, Geo_ano = map(lambda x: np.where(positive_flow_mask, x, np.nan), [T_Climo, U_Climo, V_Climo, UV_Climo, Geo_Climo, Geo_ano])
# Calculate differentials, Coriolis parameter, and constants
dlon = np.gradient(lon, axis=2)
dlat = np.gradient(lat, axis=1)
coslat = np.cos(lat)
sinlat = np.sin(lat)
f = 2 * omega * sinlat # Coriolis parameter
if len(level) > 1: # If there are multiple pressure levels
dlev = np.gradient(-sclhgt * np.log(level / 1000.0), axis=0)
N2 = gc * (level / 1000.0) ** 0.286 / sclhgt * np.gradient(T_Climo * (1000.0 / level) ** 0.286, axis=0) / dlev
# Compute streamfunction PSI
PSIa = Geo_ano / f
# Compute gradients
dzdlon = np.gradient(PSIa, axis=2) / dlon
dzdlat = np.gradient(PSIa, axis=1) / dlat
# Second-order derivatives
ddzdlonlon = np.gradient(dzdlon, axis=2) / dlon
ddzdlonlat = np.gradient(dzdlon, axis=1) / dlat
ddzdlatlat = np.gradient(dzdlat, axis=1) / dlat
if len(level) > 1:
dzdlev = np.gradient(PSIa, axis=0) / dlev
ddzdlonlev = np.gradient(dzdlon, axis=0) / dlev
ddzdlatlev = np.gradient(dzdlat, axis=0) / dlev
# Compute flux components in the x and y directions
xuterm = dzdlon * dzdlon - PSIa * ddzdlonlon
xvterm = dzdlon * dzdlat - PSIa * ddzdlonlat
yuterm = xvterm
yvterm = dzdlat * dzdlat - PSIa * ddzdlatlat
if len(level) > 1:
zuterm = dzdlon * dzdlev - PSIa * ddzdlonlev
zvterm = dzdlat * dzdlev - PSIa * ddzdlatlev
# Calculate flux components
coef = 0.5 * level * coslat / 1000.0 / UV_Climo
Fx = coef * (xuterm * U_Climo / (re * coslat) ** 2 + xvterm * V_Climo / (re ** 2 * coslat))
Fy = coef * (yuterm * U_Climo / (re ** 2 * coslat) + yvterm * V_Climo / re ** 2)
if len(level) > 1:
Fz = coef * (f ** 2 / N2 * (zuterm * U_Climo / (re * coslat) + zvterm * V_Climo / re))
# Convert results to DataArray
Fx = xr.DataArray(Fx, dims=('level', 'lat', 'lon'), coords=data_coords)
Fy = xr.DataArray(Fy, dims=('level', 'lat', 'lon'), coords=data_coords)
if len(level) > 1:
Fz = xr.DataArray(Fz, dims=('level', 'lat', 'lon'), coords=data_coords)
return Fx, Fy, Fz
else:
return Fx, Fy
EP 通量与 TN 波通量的区别：
EP通量适用于诊断Rossby波传播中的较小振幅波动，它对于反映部分异常波动能量传递较为有效，但在诊断复杂背景场或瞬时波动时显得不足。 TN波通量是对Plumb波作用通量的改进，能够更好地描述复杂流场中不均匀流的波动传播，特别适合于诊断大尺度的Rossby波异常。它在更复杂的流场背景下能提供更全面的波动传播信息。
References
Takaya, K., and H. Nakamura, 2001. A formulation of a phase-independent wave-activity flux for stationary and migratory quasigeostrophic eddies on a zonally varying basic flow. Journal of the Atmospheric Sciences, 58(6): 608-627. SHI Chunhua, JIN Xing, LIU Renqiang,2017. The differences in characteristics and applicability among three types of Rossby wave activity flux in atmospheric dynamics[J]. Trans Atmos Sci,40(6):850-855. DOI:10.13878/j. cnki. dqkxxb.20161023012. Trenberth K E, 1986. An assessment of the impact of transient eddies on the zonal flow during a blocking episode using localized Eliassen-Palm flux diagnostics [J]. J Atmos Sci, 43: 2070-2087. https://www.atmos.rcast.u-tokyo.ac.jp/nishii/programs/index.html
往期回顾
Python | MJO | 位相图 Python | 半球间经度转换 Python | 大气科学 | 偏相关 37 个 CMIP6 模型对流参数化方案 Python | 解决 cmaps 库报错的问题 Python | 批量下载 NCEP2 再分析数据 Python | 北大西洋涛动 | NAO 指数 | EOF
