| 查看: 837 | 回复: 0 | |||
[求助]
请高手指教
|
|
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 ] |
回复此楼» 猜你喜欢
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有5人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有6人回复
-
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有5人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有9人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有5人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有6人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有5人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有5人回复
-
售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急
已经有7人回复
找到一些相关的精华帖子，希望有用哦~
关于N-苯基马来酰亚胺合成后点板出现2个点。怎么处理？求高手指教下啊
已经有11人回复
- 关于液相的一个奇怪问题，请各位高手指教~~~~
已经有30人回复
- 【讨论】大家来交流一下二氯亚砜卤化 的经验嘛，遇到了一些问题，求高手指教啊
已经有10人回复
- 【求助/交流】半定量PCR啊，做不出结果啊。。。高手指教啊...
已经有5人回复
- 【求助/交流】PCR没有结果，希望高手指教
已经有29人回复
- 【求助】本人首次做DSC，不知如何分析，求高手指教。
已经有8人回复
- 研二怎么上，我该干点什么，请高手指教，我很焦虑
已经有11人回复
- 【求助】高手请进~~重谢
已经有8人回复
- 请高手指教，讲师一年最多能参加几个国家基金啊
已经有28人回复
- 【求助】请高手指教,有图求真相——Pt/C催化剂稳定性和CV曲线问题
已经有8人回复
- 【求助】我的玻碳电极出来什么问题？请高手指教
已经有33人回复
- 【交流】请各位高手指教中药饮片如何保藏过夏季
已经有7人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
