24小时热门版块排行榜

返回列表

jack小闵

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 9
帖子: 6
在线: 52分钟
虫号: 1145421
注册: 2010-11-13
性别: GG
专业: 化学环境污染与健康

[求助] 请高手指教

function nitrogen
%水田氮素动态模拟代码，假设灌溉模式为持续灌溉。间歇灌溉模式尚须考虑水分动态模拟问题，有待进一步开发。气象数据采用南京2010年5.1-8.14。
%假定生育期148天。
clc
clear all
data=xlsread('climate.xls'); %引用外部气象数据
x=data(:,1); %第几天
Tmean=data(:,2);
Tmax=data(:,3);
Tmin=data(:,4);
prec=data(:,5);%降雨量
wind=data(:,6);
produ=8100;%产量
orgN=3600; % 土壤有机氮初始量,g/kg
NH40=80; % 土壤氨氮初始残留量,g/kg
NO30=40; % 土壤硝氮初始残留量,g/kg
urea=300;%尿素施用量,kg/hm2
Kh=0.5; %尿素水解速率
Kv=0.003; %氨挥发速率。文献中这个数值约为0.05，但根据该值模拟结果，挥发量远大于氮肥使用量，不合理，可能与本文设定的植株吸氮量模拟方法有关。
Kn=0.025; %氨态硝化速率
Km=0.01; %有机氮矿化速率
Ko=0.05; %铵根固定速率
Kd=0.25; %反硝化速度，1/d
lati=31/360*2*pi;%纬度
for i=[1:1:148]
%以下计算总蒸散量
J(i)=119+i;
seta(i)=0.049*sin(0.0172*J(i)-1.39);%日倾角
dr(i)= 1+0.033*cos(2*pi/365*J(i));%日地相对距离
Ws(i)= acos(-tan(lati)*tan(seta(i)));%lati是纬度
Ra(i)= 37.6* dr(i)* (Ws(i)*sin(seta(i))*sin(lati)+cos(seta(i))*cos(lati)*sin(Ws(i)));
ET(i)=0.0023*Ra(i)*(Tmean(i)+17.8)*(Tmax(i)-Tmin(i))^0.5/(2.501-0.002361*Tmean(i));%日蒸散量，mm
end
ETC=sum(ET);
%以上计算总蒸散量
%以下计算土壤氮动态各指标
for i=[1:1:148]
NH4(1)=NH40+urea*(1- exp(-Kh));
NO3(1)=NO30+urea*(1- exp(-Kh))*(1-exp(Kn));
NH41(i)=urea*(1- exp(-Kh));%尿素水解产生的氨氮
urea=urea*exp(-Kh)
NH42(i)=orgN*(1-exp(-Km));%有机质矿化产生的氨氮
orgN=orgN*exp(-Km);%有机氮实时量。
NH43(i)=ET(i)/ETC*produ*0.016*NH4(i)/(NH4(i)+NO3(i));%植株吸氨量。两种方法估计，一是根据蒸散量*土壤中氮素含量；二是蒸散量*单位蒸散量吸氮量；经模拟，第一种方法生育前期吸氮量巨大，
%本研究采用第二者方法确定，吸收量取决于铵态和硝态的比例。其合理性尚待验证，日吸氮量应取决于日生育量，它取决于温度（蒸散量）还是土壤氮素含量？
Kdian=(-0.09018-2729.2/(Tmean(i)+273))^10;%电解常数，假定土壤
Khen=(2.39*10^5/(Tmean(i)+273))*exp(-4151/(Tmean(i)+273));%享利系数，气液相间分配系数。
if prec(i)<50
TAN(i)=NH4(i)/(2250*0.286+500+500+prec(i)*10)%液铵+气铵总浓度。2250为表层20cm土壤干重量，前一500为饱和土壤中水重量，后一500为田面水5cm重量。应用于间歇灌溉灌或旱地农作系统时需计算水分动态。
end
if prec(i)>50
TAN(i)=NH4(i)/(2250*0.286+1000)%液铵+气铵总浓度
end
NH3Y(i)=TAN(i)/(1+10^(-6/Kdian));%液铵浓度，6为土壤pH的默认值，需根据土壤类型确定。
NH3Q(i)=NH3Y(i)*Khen;%可挥发气态铵浓度
NH44(i)=NH3Q(i)*48.4*wind(i)^0.8*(Tmean(i)+273)^(-1.4)*10000*60*60*24/1000;%氨挥发量。
%NH44(i)=NH4(i)*(1-exp(-Kv*(Tmean(i)/20))^1.4);%另一种氨挥量的常见方法，假定铵挥发呈一级动力学，假定20℃时挥发速率系数为0.005，该系数与温度呈幂指数增长关系，即Kv(t)=(t/20)^1.4。
%缺陷是没有考虑水分的影响，实际上，氨挥发量主要取决于田面水铵的浓度、温度和风速。田面水铵浓度与灌溉和降水密切相关
NO31(i)=NH4(i)*(1-exp(-Kn));%铵态硝化产生的硝态氮
NO32(i)= NO3(i)*(1-exp(-Kd));%反硝化损失的硝态氮
NO320(i)=300*0.0164/365;%N2O排放量。300为施氮量。假定N2O年排放率为土壤氮素的0.0164，平均到每天。尚须进行温度和湿度的校正，它们是N2O排放率的决定因素。
NO33(i)=ET(i)/ETC*produ*0.016*NO3(i)/(NH4(i)+NO3(i));%植株吸氮量，计算方法同铵态氮。
NO34(i)=NO3(i)/2250*0.5*10;%硝态氮淋失量。
NO35(i)=0
NH45(i)=0
if prec(i)>50
NO35(i)=(prec(i)-50)*10*NO3(i)/(2250+(data(i,5)+50)*10);%硝态氮径流量。
NH45(i)=(prec(i)-50)*10*NH4(i)/(2250+(data(i,5)+50)*10);%铵态氮径流量。
end
N451(i)=sum(NH45);N351(i)=sum(NO35);%土壤铵态和硝态累计径流损失量
NH4(i+1)=NH40+sum(NH41(i))+sum(NH42(i))-sum(NH43(i))-sum(NH44(i))-N451(i); %土壤氨氮实时量
NO3(i+1)=NO30+sum(NO31(i))-sum(NO32(i))-sum(NO33(i))-sum(NO34(i))-N351(i); %土壤硝氮实时量
N1(i)=NO3(i)/2250;
N2(i)=NH4(i)/2250;
N3(i)=sum(NH43+NO33);
N4(i)=sum(NH44);
N5(i)=sum(NO35+NH45);
N6(i)=sum(NO34);
N7(i)=sum(NO320);
i=i+1;
end
N8=sum(ET);
figure(1); plot(x,N1,'-',x,N2,'-');legend('土壤硝态氮含量,mg/L','土壤铵态氮含量,mg/L')
figure(2); plot(x,ET,'-');legend('日蒸散量/mm')
figure(3); plot(x,NO31,'-');legend('硝氮日产生量,kg/hm2')
figure(4); plot(x,(NH43+NO33),'-',x,NH44,'-');legend('植株日吸氮量kg/hm2','氨日挥发量kg/hm2')
figure(5); plot(x,NO35,'-',x,NH45,'-');legend('硝氮径流量,kg/hm2','铵态径流量kg/hm2')
figure(6); plot(x,NO34);legend('硝氮日淋失量')
figure(7); plot(x,N3,'-',x,N4,'-',x,N5,'-',x,N6,'-',x,N7,'-');legend('植株累计吸氮量,kg/hm2','氨累计挥发量kg/hm2','氮素累计径流损失量','硝氮累计淋溶损失量')
figure(8); plot(x,N7,'-');legend('N2O累计排放量')
问题有以下几点：第一，肥料为尿素，如何加上硝酸铵、硫酸铵和碳酸氢铵的计算代码、第二，基施中如何加上追肥的计算过程

[ Last edited by zongyu_1986 on 2011-12-13 at 18:47 ]

回复此楼