MATLAB通量计算实习.docx
《MATLAB通量计算实习.docx》由会员分享，可在线阅读，更多相关《MATLAB通量计算实习.docx（31页珍藏版）》请在冰豆网上搜索。
MATLAB通量计算实习
通量计算实习
资料说明：
文件1：
湍流资料
文件2：
辐射、土壤温度、土壤热通量、土壤含水量以及降水记录
研究对象：
6月和8月分别选取一天资料。
利用湍流资料中的超声风速、温度、水汽、二氧化碳观测计算湍流感热通量、潜热通量，并结合辐射观测和土壤热通量观测资料分析能量平衡闭合特征。
资料处理步骤：
1.首先修改datapartition文件，将原始的湍流资料分割为每半个小时一个文件的资料；
2.修改湍流处理程序，处理资料。
资料处理具体要求：
1.野点剔除采用九点平滑法，
视为野点；
2.坐标旋转采用二维坐标旋转（uv,uw）；
3.考虑对感热通量的湿度订正（schotanus订正）和对潜热通量、CO2通量的Webb订正（WPL订正）。
4.时间平均方案：
半小时。
研究内容：
（1）订正后的感热通量H、潜热通量LE、净辐射Rn、的日变化图；
（2）订正后的CO2通量的日变化；
（3）订正后的感热通量（y轴）与未订正的感热通量（x轴）的线性拟合；给出拟合方程。
（4）订正后的潜热通量（y轴）与未订正的潜热通量（x轴）的线性拟合；给出拟合方程。
对上述图表给出简要分析说明，并将6月份选取的某天和8月份的某天进行对比分析。
作业提交内容：
1.订正后的感热通量H、潜热通量LE、净辐射Rn、土壤热通量G的日变化图；
A.2006.6.8:
分析：
（1.数据处理pflux68.f90程序。
2.画图picfor5.R程序。
3.数据资料：
flux608.dat，pflux608.dat）
净辐射、潜热通量、感热通量在14:
00左右有极大值，而土壤热通量在14:
00左右有极小值。
它们都有明显的日变化。
B.2006.8.12:
分析：
（1.数据处理flux2.f90。
2.画图PIC5.R。
3.数据资料：
flux812.dat，pflux812.dat）
该天的16:
30左右，20:
00左右，21:
00-22:
00存在降雨。
净辐射、潜热通量、感热通量、土壤热通量日变化明显，且都存在极大值。
6.8与8.12的比较
净辐射：
在16:
30左右，8.12与6.8相比较有明显的下降，这是降水的影响。
潜热通量：
8．12的潜热通量比6.8明显多，因为降水前大气的含水量及脉动较大。
感热通量：
土壤热通量：
上述4张图的M文件内容如下：
clear
clc
loadd:
\tongliang\ppflux608.dat
loadd:
\workplace\pflux812.dat
g608=ppflux608（:
3）;
g812=pflux812（:
3）;
rn608=ppflux608（:
4）;
rn812=pflux812（:
4）;
le_wpl608=ppflux608（:
6）;
le_wpl812=pflux812（:
6）;
hs608=ppflux608（:
7）;
hs812=pflux812（:
7）;
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'）;
set（p1,'Color','blue','LineWidth',3）
holdoff
p1=plot（time,rn608,'-b'）;
xlabel（'ʱ¼ä'）
ylabel（'¾»·øÉä'）
set（p1,'Color','red','LineWidth',2）
holdon
p1=plot（time,rn812,'-r'）;
set（p1,'Color','blue','LineWidth',3）
holdoff
p1=plot（time,le_wpl608,'-b'）;
xlabel（'ʱ¼ä'）
ylabel（'DZÈÈͨÁ¿'）
set（p1,'Color','yellow','LineWidth',2）
holdon
p1=plot（time,le_wpl812,'-r'）;
set（p1,'Color','blue','LineWidth',3）
holdoff
p1=plot（time,hs608,'-b'）;
xlabel（'ʱ¼ä'）
ylabel（'¸ÐÈÈͨÁ¿'）
set（p1,'Color','green','LineWidth',2）
holdon
p1=plot（time,hs812,'-r'）;
set（p1,'Color','blue','LineWidth',3）
2.订正后的CO2通量的日变化：
A.2006.06.08：
分析：
（1.画图picfor5.R。
2.数据资料：
flux812.dat）
在上午8点左右二氧化碳通量有极大值，夜间有正的通量，白天9点到19点有负的通量。
这是为什么？
。
。
。
B.2006.08.12：
分析：
（1.画图PIC5.R。
2.数据资料：
pflux812.dat）
二氧化碳日变化很剧烈。
也是天气的影响？
？
？
？
3.订正后的感热通量（y轴）与未订正的感热通量（x轴）的线性拟合；给出拟合方程。
A.2006.06.08：
分析：
（1.画图pic68.R。
2.数据资料：
Nuistflux.dat）
B.2006.08.12:
分析：
（1.画图pic812.R。
2.数据资料：
aNuistflux.dat）
4.订正后的潜热通量（y轴）与未订正的潜热通量（x轴）的线性拟合；给出拟合方程。
A.2006.06.08
分析：
（1.画图pic812.R。
2.数据资料：
Nuistflux.dat）
B.2006.08.12
分析：
（1.画图pic812.R。
2.数据资料：
aNuistflux.dat）
5.能量平衡闭合情况
A.2006.06.08：
fen
分析：
（1.画图：
拟合.R。
2.数据资料：
aNuistflux.dat）
B.2006.08.12：
分析：
（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",xlab="潜热通量订正后",ylab="潜热通量订正前"）
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））
title（"6.8拟合"）
plot（tt,f_co2,type="b",col="green",main="6.8二氧化碳通量",xlab="时间",ylab="二氧化碳通量"）
2.picfor5.R
c（）
x<-read.table（"d:
/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="b",lwd=2,xlab="时间",ylab="二氧化碳通量",col=6）
title（"二氧化碳通量日变化"）
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="时间",ylab="净辐射",col=5）
title（"净辐射日变化"）
text（5,400,"2006.06.08"）
plot（time,g,type="l",lwd=2,col=2,ylim=c（-200,600）,xlab="时间",ylab="土壤热通量"）
title（"土壤热通量日变化"）
text（12,400,"2006.06.08"）
plot（time,le_wpl,type="l",lwd=3,col=3,ylim=c（-200,600）,xlab="时间",ylab="潜热通量"）
title（"潜热通量日变化"）
text（12,400,"2006.06.08"）
plot（time,hs,type="l",lwd=2,col=4,ylim=c（-200,600）,xlab="时间",ylab="感热通量"）
title（"感热通量日变化"）
text（12,400,"2006.06.08"）
3.拟合.R
c（）
x<-read.table（"d:
/tongliang/ppflux608.dat"）
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="有效能量",ylab="湍流能量",col="blue",main="6.8拟合"）
abline（lm.rh,col="red"）
text（250,400,expression（湍流能量==0.55403+0.61723*有效能量））
4.shuju1.f90（处理8月份资料的file812.f90与此略有不同）
programmain
implicitnone
integer,external:
:
raw
real,external:
:
mean
real,external:
:
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='d:
\filename.txt',status='old'）
open（9,file='d:
\Nuistflux.dat',status='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='old'）
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）
enddo
dom=1,12!
去野点
if（m==12.or.m==10）then!
unitless
continue
endif
callkill_wp（tsdata（:
m）,k）
enddo
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）
enddo
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."
endif
close（10+i）
enddo
close（9）
end
subroutinekill_wp（sz,siz）
implicitnone
real,external:
:
mean
real,external:
:
std
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:
i+4）,9）
if（abs（sz（i）-e）>2.5*sig）then
sz（i）=0.5*sz（i-1）+0.5*sz（i+1）
endif
enddo
endsubroutine
subroutinerotate_uv（u,v,w,alph,siz）
implicitnone
real,external:
:
mean
integer:
:
siz,i
real:
:
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）
enddo
alph=direc
end
subroutinerotate_uw（u,v,w,bet,siz）
implicitnone
real,external:
:
mean
integer:
:
siz,i,beta
real:
:
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）
umean=mean（u,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
doi=1,siz
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）
u（i）=u_tmp（i）
v（i）=v_tmp（i）
w（i）=w_tmp（i）
enddo
bet=direc
endsubroutine
integerfunctionraw（b）
implicitnone
!
realz（14）
character*40la
integerj
character*30:
:
b
raw=0
open（10,file=b,status='old'）
dowhile（.true.）
!
read（10,*,end=100）（z（j）,j=1,14）
read（10,*,end=100）la
raw=raw+1
enddo
100continue
close（10）
end
subroutinetrend（sz,siz）
implicitnone
real,external:
:
mean
i
升级会员 