AtmosSci

發信人: mark.bbs@casamia.twbbs.org (Mark) 看板: AtmosSci
標 題: Re: 問一個問題：Q1 Q1 的算法
發信站: casamia (04/07/05 13:53:47 Thu)

我把程式計算的部分貼上來，有興趣的人可以看一看：
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
do k=1,L
do j=n1,n2
do i=m1,m2
s(i,j,k)= (cp*t(i,j,k)) + z(i,j,k)
if (rh(i,j,k) .LE. 0.) rh(i,j,k)=0.0001
es(i,j,k)=6.112*( exp( 17.67*( (t(i,j,k)-273.15)
+ /(t(i,j,k)-273.15+243.5) ) ) )
e(i,j,k)= (es(i,j,k)*rh(i,j,k))/100.
r(i,j,k)=0.622*( (e(i,j,k)/(pr(k)-e(i,j,k))) )
LHC(i,j,k)=1000000.*( 2.501-(0.00237*(t(i,j,k)-273.15)) )
enddo
enddo
enddo
do k=1,L
do j=n1+1,n2-1
do i=m1+1,m2-1
svdiv(i,j,k)=
+ (s(i,j,k)*( ( (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(i,j,k)=
+ (r(i,j,k)*( ( (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)*((r(i+1,j,k)-r(i-1,j,k))/(2.*dx(j))) )
+ +( v(i,j,k)*((r(i,j-1,k)-r(i,j+1,k))/(2.*dy)) )
+ )
enddo
enddo
enddo
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
num=(ns2-ns1+1)*(ms2-ms1+1)
do it=1,time
do k=1,L
ave_s(k,it) = 0.
ave_svdiv(k,it) = 0.
ave_w(k,it) = 0.
ave_r(k,it) = 0.
ave_rvdiv(k,it) = 0.
do j=ns1,ns2
do i=ms1,ms2
ave_s(k,it) =ave_s(k,it) + s(i,j,k,it)
ave_svdiv(k,it)=ave_svdiv(k,it) + svdiv(i,j,k,it)
ave_r(k,it) =ave_r(k,it) + (LHC(i,j,k,it)*r(i,j,k,it))
ave_rvdiv(k,it)=ave_rvdiv(k,it)
+ + (LHC(i,j,k,it)*rvdiv(i,j,k,it))
ave_w(k,it) =ave_w(k,it) + w(i,j,k,it)
enddo
enddo
ave_s(k,it) =ave_s(k,it) /float(num)
ave_svdiv(k,it)=ave_svdiv(k,it)/float(num)
ave_q(k,it) =ave_q(k,it) /float(num)
ave_qvdiv(k,it)=ave_qvdiv(k,it)/float(num)
ave_w(k,it) =ave_w(k,it) /float(num)
enddo
enddo

do it=2,time-1
do k=2,L-1
q11(k,it)=(ave_s(k,it+1)-ave_s(k,it-1)) / (2.*6.*3600.)
q12(k,it)=ave_svdiv(k,it)
a=pr(k)-pr(k+1)
b=pr(k-1)-pr(k)
q13(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)
+ )
q21(k,it)=-1.*((ave_q(k,it+1)-ave_q(k,it-1))/(2.*6.*3600.))
q22(k,it)=-1.*(ave_qvdiv(k,it))
q23(k,it)=-1.*( ( (a*ave_q(k-1,it)*ave_w(k-1,it))
+ -(b*ave_q(k+1,it)*ave_w(k+1,it))
+ +((b-a)*ave_q(k,it)*ave_w(k,it)) )
+ / (2.*a*b)
+ )
enddo
enddo
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊

Mark

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