MATLAB通量计算实习Word文档格式.docx
- 文档编号:22320439
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:31
- 大小:104.83KB
MATLAB通量计算实习Word文档格式.docx
《MATLAB通量计算实习Word文档格式.docx》由会员分享,可在线阅读,更多相关《MATLAB通量计算实习Word文档格式.docx(31页珍藏版)》请在冰豆网上搜索。
6.8与8.12的比较
净辐射:
在16:
30左右,8.12与6.8相比较有明显的下降,这是降水的影响。
潜热通量:
8.12的潜热通量比6.8明显多,因为降水前大气的含水量及脉动较大。
感热通量:
土壤热通量:
上述4张图的M文件内容如下:
clear
clc
loadd:
\tongliang\ppflux608.dat
\workplace\pflux812.dat
g608=ppflux608(:
3);
g812=pflux812(:
rn608=ppflux608(:
4);
rn812=pflux812(:
le_wpl608=ppflux608(:
6);
le_wpl812=pflux812(:
hs608=ppflux608(:
7);
hs812=pflux812(:
time=0:
0.5:
23.5;
p1=plot(time,g608,'
-b'
);
xlabel('
Ê
±
¼
ä
'
)
ylabel('
Í
Á
È
À
¨
¿
set(p1,'
Color'
'
red'
LineWidth'
2)
holdon
p1=plot(time,g812,'
-r'
blue'
3)
holdoff
p1=plot(time,rn608,'
¾
»
·
ø
É
p1=plot(time,rn812,'
p1=plot(time,le_wpl608,'
Ç
yellow'
p1=plot(time,le_wpl812,'
p1=plot(time,hs608,'
¸
Ð
green'
p1=plot(time,hs812,'
2.订正后的CO2通量的日变化:
A.2006.06.08:
(1.画图picfor5.R。
2.数据资料:
flux812.dat)
在上午8点左右二氧化碳通量有极大值,夜间有正的通量,白天9点到19点有负的通量。
这是为什么?
。
B.2006.08.12:
(1.画图PIC5.R。
pflux812.dat)
二氧化碳日变化很剧烈。
也是天气的影响?
?
3.订正后的感热通量(y轴)与未订正的感热通量(x轴)的线性拟合;
(1.画图pic68.R。
Nuistflux.dat)
B.2006.08.12:
(1.画图pic812.R。
aNuistflux.dat)
4.订正后的潜热通量(y轴)与未订正的潜热通量(x轴)的线性拟合;
A.2006.06.08
Nuistflux.dat)
B.2006.08.12
5.能量平衡闭合情况
fen
(1.画图:
拟合.R。
(1.画图nihe812.R,2.数据资料:
pflux812.dat)
能量不闭合的原因:
(1)通量观测中的采样误差(footprint问题)
(2)测量仪器的系统误差
(3)其他能量源汇项的忽略
主要参考文献:
[1]WebbEK,PearmanGI,LeuningR.Correctionoffluxmeasurementsfordensityeffectsduetoheatandwatervapourtransfer.QuarterlyJournaloftheRoyalMeteorologicalSociety,1980,106:
85-100
[2]SchotanusP,NieuwstadtFTM,DeBruinHAR.Temperaturemeasurementswithasonicanemometeranditsapplicationtoheatandmoisturefluxes.BoundaryLayerMeteorology,1983,26:
81-93
[3]LiuH,PetersG,FokenT.Newequationsforsonictemperaturevarianceandbuoyancyheatfluxwithanomnidirectionalsonicanemometer.Boundary-LayerMeteorology,2001,100:
459-468
[4]于贵瑞,张雷明,孙晓敏,等.亚洲区域陆地生态系统碳通量观测研究进展.中国科学(D辑),2004,34(增Ⅱ):
15-29
[5]于贵瑞,孙晓敏等.陆地生态系统通量观测的原理与方法.高等教育出版社.2006.
主要程序:
1.pic812.R
c()
x<
-read.table("
d:
/tongliang/Nuistflux.dat"
header=FALSE)
hwts<
-x[,1]
hwts_c<
-x[,2]
le_raw<
-x[,3]
le<
-x[,4]
bo_raw<
-x[,5]
f_co2<
-x[,6]
tt<
-seq(from=0,to=23.5,by=.5)
plot(hwts,hwts_c,col="
green"
xlab="
感热通量订正前"
ylab="
感热通量订正后"
lm.h<
-lm(hwts_c~1+hwts)
summary(lm.h)
abline(lm.h,col="
yellow"
text(100,250,expression(hwts_c==-1.120311+0.956567*hwts))
title("
6.8拟合"
plot(le,le_raw,col="
red"
潜热通量订正后"
潜热通量订正前"
lm.le<
-lm(le_raw~1+le)
summary(lm.le)
abline(lm.le,col="
blue"
text(300,50,expression(le_raw==6.66392+0.54679*le))
plot(tt,f_co2,type="
b"
col="
main="
6.8二氧化碳通量"
时间"
二氧化碳通量"
2.picfor5.R
/tongliang/ppflux608.dat"
g<
-x$V3
rn<
-x$V4
fc_wpl<
-x$V5
le_wpl<
-x$V6
hs<
-x$V7
time<
-seq(from=0,to=23.5,by=0.5)
plot(fc_wpl~time,type="
lwd=2,xlab="
col=6)
二氧化碳通量日变化"
text(12,0.5,"
2006.06.08"
add=FALSE
windows(width=10,height=10)
par(mfrow=c(2,2))
plot(rn~time,type="
l"
lwd=2,ylim=c(-200,600),xlab="
净辐射"
col=5)
净辐射日变化"
text(5,400,"
plot(time,g,type="
lwd=2,col=2,ylim=c(-200,600),xlab="
土壤热通量"
土壤热通量日变化"
text(12,400,"
plot(time,le_wpl,type="
lwd=3,col=3,ylim=c(-200,600),xlab="
潜热通量"
潜热通量日变化"
plot(time,hs,type="
lwd=2,col=4,ylim=c(-200,600),xlab="
感热通量"
感热通量日变化"
3.拟合.R
h_add_le<
-x$V1
rn_g<
-x$V2
lm.rh<
-lm(h_add_le~1+rn_g)
summary(lm.rh)
plot(rn_g,h_add_le,xlab="
有效能量"
湍流能量"
abline(lm.rh,col="
text(250,400,expression(湍流能量==0.55403+0.61723*有效能量))
4.shuju1.f90(处理8月份资料的file812.f90与此略有不同)
programmain
implicitnone
integer,external:
:
raw
real,external:
mean
std
integeri,j,k,m,n,l
character*30a(48),ll(4)
character*33ll1,ll2
realtsdata(18000,12),u(18000),v(18000),w(18000),ts(18000),rou_w(18000),press(18000)
realt_fw(18000),q(18000),rou_co2(18000)
realrd,rv,ma,mv,g,alpha,beta,p_ave,ts_ave,t_ave,ta_ave,e_ave,q_ave,rou_dry,rou_vapor
realrou,cpd,cp,lv,uw,wts,wtp,wq_raw,le_raw,hwts,hwts_c,wts_c,bo_raw,le,wq
realF_co2,co2_ave,wco2
rd=287.05
rv=461.5
ma=29
mv=18
g=9.8
open(10,file='
\filename.txt'
status='
old'
open(9,file='
\Nuistflux.dat'
replace'
)!
(打开通量的资料)
read(10,*)(a(i),i=1,48)
close(10)
doi=1,48
tsdata=0
k=raw(a(i))-4!
行数
write(*,*)'
hang..'
k
open(10+i,file=a(i),status='
dol=1,4!
去掉前四行
read(10+i,'
(a30)'
advance='
yes'
)ll(i)
enddo
dom=1,k
read(10+i,*)ll1,ll2,(tsdata(m,n),n=1,12)
dom=1,12!
去野点
if(m==12.or.m==10)then!
unitless
continue
endif
callkill_wp(tsdata(:
m),k)
doj=1,k
u(j)=tsdata(j,1)
v(j)=tsdata(j,2)
w(j)=tsdata(j,3)
ts(j)=tsdata(j,6)
rou_w(j)=tsdata(j,5)
press(j)=tsdata(j,7)
t_fw(j)=tsdata(j,11)
rou_co2(j)=tsdata(j,4)
callrotate_uv(u,v,w,alpha,k)!
坐标旋转uv合成U,新坐标中v为0
callrotate_uw(u,v,w,beta,k)!
同上
ts_ave=mean(ts,k)!
摄氏度求平均
ts_ave=ts_ave+273.16!
转换成绝对温度(K)
p_ave=mean(press,k)
t_ave=mean(t_fw,k)
rou_vapor=mean(rou_w,k)!
水汽密度,g/m3
rou_vapor=rou_vapor/1000!
kg/m3
e_ave=rou_vapor*rv*(t_ave+273.16)/1000!
水汽压
rou_dry=1000*(p_ave-e_ave)/(t_ave+273.16)/rd!
kg/m3干空气密度
rou=1000*(rou_dry+rou_vapor)!
空气密度(g/m3)湿空气
q=rou_w/rou!
比湿g/g水汽密度(观测得)/湿空气密度
q_ave=mean(q,k)
ta_ave=ts_ave/(1+0.51*q_ave)!
绝对温度(K)超声虚温订正
co2_ave=mean(rou_co2,k)
calltrend(u,k)
calltrend(w,k)
calltrend(ts,k)
calltrend(q,k)
calltrend(rou_co2,k)!
!
*************co2
uw=mean(u*w,k)
wts=mean(w*ts,k)
wq_raw=mean(w*q,k)
wco2=mean(w*rou_co2,k)
cpd=1004.67
cp=cpd*(1+0.84*q_ave)
lv=(2.501-0.00237*t_ave)*10**6!
摄氏度
le_raw=rou*lv*wq_raw/1000!
潜热通量
hwts=rou*cp*wts/1000!
感热通量
wts_c=wts-0.51*ta_ave*wq_raw
hwts_c=rou*cp*wts_c/1000!
订正后
bo_raw=hwts_c/le_raw
wq=(1+ma/mv*q_ave)*(1+lv/cp*(q_ave/t_ave)*bo_raw)*wq_raw
le=rou*lv*wq/1000!
潜热通量。
F_co2=wco2+1.6*co2_ave/rou_dry*wq_raw+(1+1.6*q_ave)*co2_ave*wts/ts_ave
write(9,'
(1x,6f)'
)hwts,hwts_c,le_raw,le,bo_raw,F_co2
if(mod(i,20)==0)then
print*,i,"
datashaveprocessed."
close(10+i)
enddo
close(9)
end
subroutinekill_wp(sz,siz)
integer:
siz,i
real:
sz(siz),sig,e
doi=5,siz-4
sig=std(sz(i-4:
i+4),9)
e=mean(sz(i-4:
if(abs(sz(i)-e)>
2.5*sig)then
sz(i)=0.5*sz(i-1)+0.5*sz(i+1)
endsubroutine
subroutinerotate_uv(u,v,w,alph,siz)
u(siz),v(siz),w(siz),alph,pi=3.1415926,alpha,umean,vmean,us,sinphi,cosphi,direc,yaw1,yaw2,yaw3,yaw4,yaw5,yaw6,yaw7,yaw8,yaw9,u_tmp(siz),v_tmp(siz),w_tmp(siz)
umean=mean(u,siz)
vmean=mean(v,siz)
us=(umean**2+vmean**2)**0.5
sinphi=vmean/us
cosphi=umean/us
direc=180*acos(cosphi)/pi
if(sinphi<
0.0)then
direc=360-direc
else
direc=direc
endif
sinphi=sin(pi*direc/180)
cosphi=cos(pi*direc/180)
yaw1=cosphi
yaw2=sinphi
yaw3=0
yaw4=-sinphi
yaw5=cosphi
yaw6=0
yaw7=0
yaw8=0
yaw9=1.0
doi=1,siz
u_tmp(i)=yaw1*u(i)+yaw2*v(i)+yaw3*w(i)
v_tmp(i)=yaw4*u(i)+yaw5*v(i)+yaw6*w(i)
w_tmp(i)=yaw7*u(i)+yaw8*v(i)+yaw9*w(i)
u(i)=u_tmp(i)
v(i)=v_tmp(i)
w(i)=w_tmp(i)
alph=direc
subroutinerotate_uw(u,v,w,bet,siz)
siz,i,beta
u(siz),v(siz),w(siz),bet,pi=3.1415926,umean,wmean,us,sintheta,costheta,direc,pitch1,pitch2,pitch3,pitch4,pitch5,pitch6,pitch7,pitch8,pitch9,u_tmp(siz),v_tmp(siz),w_tmp(siz)
wmean=mean(w,siz)
us=(umean**2+wmean**2)**0.5
sintheta=wmean/us
direc=180*asin(sintheta)/pi
sintheta=sin(pi*direc/180)
costheta=cos(pi*direc/180)
pitch1=costheta
pitch2=0
pitch3=sintheta
pitch4=0
pitch5=1.0
pitch6=0
pitch7=-sintheta
pitch8=0
pitch9=costheta
u_tmp(i)=pitch1*u(i)+pitch2*v(i)+pitch3*w(i)
v_tmp(i)=pitch4*u(i)+pitch5*v(i)+pitch6*w(i)
w_tmp(i)=pitch7*u(i)+pitch8*v(i)+pitch9*w(i)
bet=direc
integerfunctionraw(b)
realz(14)
character*40la
integerj
character*30:
b
raw=0
open(10,file=b,status='
dowhile(.true.)
read(10,*,end=100)(z(j),j=1,14)
read(10,*,end=100)la
raw=raw+1
100continue
subroutinetrend(sz,siz)
i
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 通量 计算 实习