PythonのWRF出力から3D Frontogenesisを計算する

Question

現在、私はWRF出力データの3次元フロントジェネシスを計算するためにPythonスクリプトを実行して、断面解析を実行しようとしています。私は既に、2次元のバージョンが何であるかの非常にマイナーな側面をキャプチャしているにもかかわらず、比較するためにペターソンの公式の二次元のバージョンを動作させています。

私が3D frontogenesisを計算するために使用している式は次のとおりです。勾配項がある

： ：

ここに私の2次元コード（925 hPaの）が生産するものの一例のイメージです

：

そして、ここでは、同じ2次元（925 HPA）の表面にまで補間された三次元コードです

少なくともいくつかの読みが画像上の同じ地理的領域に現れるという事実は、私が少なくとも正しい軌道にいること、そしておそらくどこかで軽度のエラーを起こしていることを示している。私が見たところでは、私のPythonコードは、ここでは、少なくともnp.gradient機能がどのように動作するかの私の理解で、正しい表示されますfrontogenesis計算コードである：参照、DXと同様

# Fetch the fields we need 
p = getvar(ncFile, "pressure") * 100 # Convert to Pa 
z = getvar(ncFile, "z") 
ua = getvar(ncFile, "ua") 
va = getvar(ncFile, "va")   
theta = getvar(ncFile, "theta") 
omega = getvar(ncFile, "omega") 

dz = np.gradient(z, axis=0) 
dp = np.gradient(p, axis=0) 

theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (np.gradient(theta, axis=0)/dz)**2) 

zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1))) 
meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1))) 
vertical_gradient = (-1 * (np.gradient(theta, axis=0)/dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1))) 

F3D = 1.08e9 * (1/theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient) 
return F3D 

をdyの用語は、NetCDFファイル自体から直接取得されます（属性DXとDYは両方とも4000mと定義されています）

私はnetCDFをインポートするgetvarを使用してデータをインポートするためにwrf-ファイル。参考として、netCDFファイルは、標準numpy配列と同様の配列構造を使用しています。 [bottom-top、south-north、east-west]

したがって、axis引数の順序は正しいはずです（z = 0、y = 1、x = 2）。

私の教員アドバイザーの一人は、グラデーションのエッジが計算の内部で問題を引き起こしていると考えていますが、各点はエッジとは独立して計算されるため、違いはありませんこれは私が100％確信していないものです。

上記の画像に示されているように計算結果が正しくない理由が分かっている人はいますか？

Answer 1

もう少し頭打ちをした後、私は問題を解決することができました。あなたの完全なユニット分析を検討することを確認するためにショーに行く！

潜在的な温度の勾配の項は、垂直軸をzと考えていましたが、オメガの垂直勾配の項はpを考えていました。オメガは圧力の垂直単位なので、潜在的な温度勾配のための私の用語の選択は間違っていた。 zからpへの導関数を交換すると、問題の前半が修正されました。

第2に、垂直方向の導関数を計算する場合、そのnumpyを考慮する必要があります。勾配は、2つ以上のデータポイントを旅していると仮定しているので、それは間違って2で結果を分割し、私は分析の内側と外側のポイントを処理するための関数ラッパーを作成しました。その後

def calc_center_difference(A, ax): 
    gradient = np.gradient(A, axis=ax) 
    gradient *= 2.0 
    if ax == 0: 
     gradient[0,:,:] /= 2.0 
     gradient[-1,:,:] /= 2.0 
    elif ax == 1: 
     gradient[:,0,:] /= 2.0 
     gradient[:,-1,:] /= 2.0 
    elif ax == 2: 
     gradient[:,:,0] /= 2.0 
     gradient[:,:,-1] /= 2.0 
    else: 
     return ValueError("Invalid axis passed to calc_center_difference") 
    return gradient 

、私はスワップラッパー関数インスタンスでdzまたはdpのいずれかを考慮して、条件が正しく計算されていることを確認した派生条件のインスタンス、およびvoila！すべてが動作します！

ここでのpythonでフル3D frontogenesis計算することができます私の以前の機能の修正形です：ちょうどリマインダーとして

# Input netcdf 
# - [bottom_top, north_south, west_east] 
def three_d_fronto(ncFile): 
    # Fetch the fields we need 
    p = to_np(getvar(ncFile, "pressure") * 100) 
    z = to_np(getvar(ncFile, "z")) 
    ua = to_np(getvar(ncFile, "ua")) 
    va = to_np(getvar(ncFile, "va"))  
    theta = to_np(getvar(ncFile, "theta")) 
    omega = to_np(getvar(ncFile, "omega")) 

    dz = calc_center_difference(z, 0) 
    dp = calc_center_difference(p, 0) 

    theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (calc_center_difference(theta, 0)/dp)**2) 
    zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1))) 
    meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1))) 
    vertical_gradient = (-1 * (calc_center_difference(theta, 0)/dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1))) 

    F3D = 1.08e9 * (1/theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient) 
    return F3D 

は、あなたがこの機能を使用するためにインストールWRF-pythonのライブラリを必要としています。