AtmosSci

發信人: mark.bbs@casamia.twbbs.org (Mark) 看板: AtmosSci
標 題: Re: 問一個問題：Q1 Q1 的算法
發信站: casamia (04/06/05 14:54:28 Wed)

首先，以EC網格資料計算乾靜能：
s=(cp*t)+gz
cp=1004., t=絕對溫度, z=重力位高度（EC資料的z場已經是 gz，所以沒有乘 g）

再來，由 RH 換算混合比(r)：
if (rh .LE. 0.) rh=0.0001
es=6.112*( exp( 17.67*( (t-273.15)/(t-273.15+243.5) )
)
)
e =(es*rh)/100.
r =0.622*( ( e/(pr-e) ) )
es為飽和水氣壓，e為水氣壓，pr為資料之高度（單位：hPa）

再算凝結潛熱釋放（LHC）：
LHC=1000000.*( 2.501-(0.00237*(t-273.15)) )

之後算 div(sV) 與 div(rV)：（多餘的括號很多，是為了確保計算正確）
svdiv= (s*( ((u(i+1,j,k)-u(i-1,j,k))/(2.*dx(j)))
+((v(i,j-1,k)-v(i,j+1,k))/(2.*dy))
)
)
+( ( u(i,j,k)*((s(i+1,j,k)-s(i-1,j,k))/(2.*dx(j))) )
+( v(i,j,k)*((s(i,j-1,k)-s(i,j+1,k))/(2.*dy)) )
)
rvdiv 類似 svdiv，就不重複貼了
dy=2.*6370000.*(4.*atan(1.))/(360./dd)
dy是南北向一個網格的距離，單位為公尺，dd是網格解析度
dx(j)=sin((j-1)*pi/(180./dd))*dy
dx是東西向一個網格的距離，此值會隨緯度而變

之後就先輸出欲分析的對流區，水平尺度約八度＊六度的小區域
輸出變數包括：s,svdiv,r,rvdiv,w,LHC
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
第二個程式先算 s、svdiv、r、rvdiv、w 的水平平均
此時已先將 r 與 rvdiv 先乘上 LHC
這些資料的矩陣為：（高度,時間）的二維矩陣
隨後 Q1 的第一項為：(ave_s(k,it+1)-ave_s(k,it-1)) / (2.*6.*3600.)
第二項為：ave_svdiv(k,it)
第三項為：( ( (a*ave_s(k-1,it)*ave_w(k-1,it))
-(b*ave_s(k+1,it)*ave_w(k+1,it))
+((b-a)*ave_s(k,it)*ave_w(k,it)) )
/(2.*a*b)
)
a=pr(k)-pr(k+1)
b=pr(k-1)-pr(k)
（此為垂直非固定間距之一階中插分）
Q2的三項與Q1的三項類似，只是把 s 換成 r ，然後把Q2的各項*(-1.)
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

有興趣的人，我可以把程式、資料傳給你玩玩
因為我真的已經不知道問題在哪了......

有人說是垂直運動要以連續方程計算而不能直些使用EC資料的垂直運動
但是我認為那個因素應該不致於使Q1與Q2的尺度差一個量級

Mark

--
※ 發信站: 卡莎米亞(bbs.as.ntu.edu.tw)
◆ From: 140.112.67.16