|
楼主 |
发表于 2018-4-28 10:01:59
|
显示全部楼层
P = [80.3,33.5,8.2,92.7,8.6,0,0,0.3,2.1,0,0,0,0,0,10.3,0,0,0,12.8,0.7,3.7,50.2,0.4,0,0,0.6,3.7,0,0,0]
EP = [4.4,0.4,3.5,1.8,0.1,3.8,3.9,3.6,2.3,4.6,5.3,4.8,5.6,5.5,4.8,4.6,6.5,6.3,4.2,3.2,1.5,0.4,1.2,3.9,2.4,4.8,5.3,5.9,6.0,6.1]
b=0.3;C = 1.0 / 6.0 ;Fc = 2.5;WUM=20;WLM=60;WDM=40;WM=120;EX=1.5;KSS=0.135;SM=40.3;KG=0.45;S=24;SMM=SM*(1+EX);FR=0.8;SMMF=SMM*(1-(1-FR)**(1/EX));SMF=SMMF/(1+EX);AU=SMMF*(1-(1-S/SMF)**(1/(1+EX)));
EU=[0]*365;EL=[0]*365;ED=[0]*365;E=[0]*365;PE=[0]*365;WU=[0]*366;WL=[0]*366
WD=[0]*366;W=[0]*366;R=[0]*366;K=[0]*366;RG=[0]*366;RS=[0]*366;RSS=[0]*366;S=[0]*366;
WU[0]=WUM;WL[0]=WLM;WD[0]=WDM;WMM=(1+b)*WM
for i in range(0,30):
#蒸散发计算
if (WU[i]+P[i])>=EP[i]:
EU[i]=EP[i]
EL[i]=0
ED[i]=0
elif (WU[i]+P[i]<EP[i]) and (WL[i]>=C*WLM):
EU[i]=WU[i]+P[i]
EL[i]=(EP[i]-EU[i])*WL[i]/WLM
ED[i]=0
elif (WU[i]+P[i]<EP[i]) and (WL[i]<C*WLM) and (WL[i]>=C*(EP[i]-EU[i])):
EU[i]=WU[i]+P[i]
EL[i]=(EP[i]-EU[i])*C
ED[i]=0
else:
EU[i]=WU[i]+P[i]
EL[i]=WL[i]
ED[i]=(EP[i]-EU[i])*C-EL[i]
E[i]=EU[i]+EL[i]+ED[i]
PE[i]=P[i]-E[i]
W[i]=WU[i]+WL[i]+WD[i]
#产流计算:
if PE[i]>0:
a=WMM*(1-(1-W[i]/WM)**(1/(1+b)))
if a+PE[i]<WMM:
R[i]=PE[i]-WM+W[i]+WM*(1-(PE[i]+a)/WMM)**(1+b)
else: R[i]=PE[i]+W[i]-WM
K[i]=R[i]/PE[i]
else:
R[i]=0
K[i]=0
#三层蓄水量计算:
if WU[i]+P[i]-EU[i]-R[i]<=WUM:
WU[i+1]=WU[i]+P[i]-EU[i]-R[i]
WL[i+1]=WL[i]-EL[i]
WD[i+1]=WD[i]-ED[i]
else:
WU[i+1]=WUM
if (WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-WUM))<=WLM:
WL[i+1]=WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-WUM)
WD[i+1]=WD[i]-ED[i]
else:
WL[i+1]=WLM
if (WD[i]-ED[i]+WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-WUM)-WLM)<=WDM:
WD[i+1]=WD[i]-ED[i]+WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-WUM)-WLM
else:
WD[i+1]=WDM
#二水源划分:
if 0<PE[i]+AU<SMMF:
RS[i]=(PE[i]-SMF+S+SMF*(1-(PE[i]+(AU)/SMMF)**(EX+1))
RSS[i]=(PE[i]+S-RS/FR)*KSS*FR
RG[i]=(PE[i]+S-RS/FR)*KG*FR
S=S+PE[i]-(RS+RSS+RG)/FR
elif PE[I]+AU>=SMMF:
RS[i] = (PE[I]-SMF+S)*FR
RSS[i] = SMF*KSS*FR
RG[i] = SM*KG*FR
S = SMF-(RSS+RG)/FR
else:
RS=0
RSS=0
RG=0
S=0
print(repr(i+1).rjust(2),'EU','=',round(EU[i],1), 'EL','=',round(EL[i],1), 'ED','=',round(ED[i],1), 'E','=',round(E[i],1), end=' ')
print('PE','=',round(PE[i],1),'WU','=',round(WU[i],1), 'WL','=',round(WL[i],1), 'WD','=',round(WD[i],1), 'W','=',round(W[i],1),end=' ')
print('R','=',round(R[i],1), 'R/PE','=',round(K[i],1), 'RG','=',round(RG[i],1), 'RS','=',round(RS[i],1))
|
|