您所在位置：主题：[讨论]一段新安江模型的程序 出错了 高手们帮忙看下啊 非常感谢
bibixiang
[专家分：0] 发布于 2010-12-07 22:01:00
PROGRAM XAJ
!原始降雨蒸发数据
REAL, DIMENSION(24):: E0=(/0.1,0.0,0.1,0.5,0.7,0.9,0.8,0.7,0.5,0.3,0.2,0.1, &
0.0,0.0,0.1,0.6,0.8,1.0,0.9,0.8,0.7,0.5,0.3,0.1/)
REAL, DIMENSION(24):: P=(/10.0,24.1,20.4,18.3,10.1,5.5,0.6,3.1,1.9,4.6,5.0,&
4.8,36.2,29.0,6.0,3.6,0.4,0.0,0.5,3.8,0.0,1.8,0.2,0.3/)
REAL,DIMENSION(3):: UH=(/0.3,0.6,0.1/)
!----------------------------------------------------------------------------------------------
!需要用到的参数
REAL,DIMENSION(24):: EM,PE,R,EU,EL,ED,E,A,FR,AU,RS,RSS,RG,QS,QSS,QG,Q
REAL,DIMENSION(25):: WU,WL,WD,W,S
REAL:: WUM=20.0,WLM=75.0,WDM=80.0,K=0.65,C=0.11,&
B=0.3,IMP=0.0,SM=20.0,EX=1.0,KG=0.3,KSS=0.41,KKG=0.99,KKSS=0.6,&
WM,WWMM,SSM,DT=2.0,QSS0=40.0,QG0=20.0,F=537.0
!----------------------------------------------------------------------------------------------
!临时变量，用来存储降雨后的临时蓄水量
REAL:: DU=0.0,DL=0.0,DD=0.0
REAL:: WUX=0.0,WLX=0.0,WDX=0.0
!-----------------------------------------------------------------------------------------------
S(1)=20.0
!三层蒸发模型计算蒸发量
WM=WUM+WLM+WDM
WU(1)=0.0
WL(1)=70.0
WD(1)=80.0
W(1)=WU(1)+WL(1)+WD(1)
WWMM=(1.0+B)*WM
!------------------------------------------------------------------------------------------------
DO I=1,24
EM(I)=K*E0(I)
!各层蓄水量变化，降水先补充上层，上层蓄满后再补充下层和深层，将降雨分配到各层
DU=WUM-WU(I)
DL=WLM-WL(I)
DD=WDM-WD(I)
IF(P(I)<=DU)THEN
WUX=P(I)+WU(I)
ELSE IF(P(I)>DU.AND.P(I)<=DU+DL)THEN
WUX=WUM
WLX=WL(I)+P(I)-DU
ELSE IF(P(I)>DU+DL.AND.P(I)<=DU+DL+DD)THEN
WUX=WUM
WLX=WLM
WDX=WD(I)+P(I)-DU-DL
ELSE
WUX=WUM
WLX=WLM
WDX=WDM
END IF
!----------------------------------------------------------------------
!三层蒸发模型计算蒸发,先蒸发上层再蒸发下层和深层，将蒸发能力分配到各层
IF(WUX>=EM(I))THEN
EU(I)=EM(I)
EL(I)=0.0
ED(I)=0.0
ELSE
EU(I)=WUX
IF(WLX>=C*WLM)THEN
EL(I)=(EM(I)-EU(I))*WLX/WLM
ED(I)=0.0
ELSE
IF(WLX>=C*(EM(I)-EU(I)))THEN
EL(I)=C*(EM(I)-EU(I))
ED(I)=0.0
ELSE
EL(I)=WLX
ED(I)=C*(EM(I)-EU(I))-EL(I)
END IF
END IF
END IF
WU(I+1)=WUX-EU(I)
WL(I+1)=WLX-EL(I)
WD(I+1)=WDX-ED(I)
E(I)=EU(I)+EL(I)+ED(I)
!---------------------------------------------------------
!产流计算
PE(I)=P(I)-E(I)
A(I)=WWMM*(1.0-(1.0-W(I)/WM)**(1.0/(1.0+B)))
IF(PE(I)<=0.0)THEN
R(I)=0
W(I+1)=W(I)+PE(I)
ELSE
IF(PE(I)+A(I)>WWMM)THEN
R(I)=PE(I)-(WM-W(I))
ELSE
R(I)= PE(I)-((WM-W(I))-WM*(1-(PE(I)+A(I))/WWMM)**(1.0+B)) !落一个*
END IF
W(I+1)=W(I)+P(I)-R(I)
END IF
!----------------------------------------------------------------------
!水源划分
SSM=SM*(1.0+EX)
AU(I)=SSM*(1.0-(1.0-S(I)/SM)**(1.0/(1.0+EX)))
IF(PE(I)>0)THEN
FR(I)=R(I)/PE(I)
IF(PE(I)+AU(I)<SSM)THEN
RS(I)=(PE(I)+S(I)-SM+SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*FR(I)
RSS(I)=(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*KSS*FR(I)
RG(I)=(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*KG*FR(I)
S(I+I)=(1.0-KSS-KG)*(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(1.0+EX))
ELSE
RS(I)=(PE(I)+S(I)-SM)*FR(I)
RSS(I)=SM*KSS*FR(I)
RG(I)=SM*KG*FR(I)
S(I+1)=(1.0-KSS-KG)*SM
END IF
ELSE
FR(I)=1.0-(1.0-W(I)/WM)**(B/(1.0+B))
RS(I)=0.0
RSS(I)=S(I)*KSS*FR(I)
RG(I)=S(I)*KG*FR(I)
S(I+1)=(1-KSS-KG)*S(I)
END IF
QS(I)=0.0;
QSS(I)=0.0;
QG(I)=0.0;
END DO
!------------------------------------------------------------------------------------
!汇流计算，地表径流采用瞬时单位线，壤中流和地下径流采用线型水库
DO I=1,24
DO J=1,3
QS(I+J-1)=QS(I+J-1)+RS(I)*UH(J)*F/(3.6*2.0)
END DO
IF(I==1)THEN
QSS(I)=QSS0*KKSS**(1/12)+RSS(I)*(1.0-KKSS**(1/12))*F/7.2
QG(I)=QGO*KKG**(1/12)+RG(I)*(1.0-KKG**(1/12))*F/7.2
ELSE
QSS(I)=QSS(I-1)*KKSS**(1/12)+RSS(I)*(1.0-KKSS**(1/12))*F/7.2
QG(I)=QG(I-1)*KKG**(1/12)+RG(I)*(1.0-KKG**(1/12))*F/7.2
END IF
Q(I)=QS(I)+QSS(I)+QG(I)
WRITE(*,*) Q(I)
END DO
!--------------------------------------------------------------------------------------
READ(*,*)
STOP
END PROGRAM XAJ
!原始降雨蒸发数据
REAL, DIMENSION(24):: E0=(/0.1,0.0,0.1,0.5,0.7,0.9,0.8,0.7,0.5,0.3,0.2,0.1, &
0.0,0.0,0.1,0.6,0.8,1.0,0.9,0.8,0.7,0.5,0.3,0.1/)
REAL, DIMENSION(24):: P=(/10.0,24.1,20.4,18.3,10.1,5.5,0.6,3.1,1.9,4.6,5.0,&
4.8,36.2,29.0,6.0,3.6,0.4,0.0,0.5,3.8,0.0,1.8,0.2,0.3/)
REAL,DIMENSION(3):: UH=(/0.3,0.6,0.1/)
!----------------------------------------------------------------------------------------------
!需要用到的参数
REAL,DIMENSION(24):: EM,PE,R,EU,EL,ED,E,A,FR,AU,RS,RSS,RG,QS,QSS,QG,Q
REAL,DIMENSION(25):: WU,WL,WD,W,S
REAL:: WUM=20.0,WLM=75.0,WDM=80.0,K=0.65,C=0.11,&
B=0.3,IMP=0.0,SM=20.0,EX=1.0,KG=0.3,KSS=0.41,KKG=0.99,KKSS=0.6,&
WM,WWMM,SSM,DT=2.0,QSS0=40.0,QG0=20.0,F=537.0
!----------------------------------------------------------------------------------------------
!临时变量，用来存储降雨后的临时蓄水量
REAL:: DU=0.0,DL=0.0,DD=0.0
REAL:: WUX=0.0,WLX=0.0,WDX=0.0
!-----------------------------------------------------------------------------------------------
S(1)=20.0
!三层蒸发模型计算蒸发量
WM=WUM+WLM+WDM
WU(1)=0.0
WL(1)=70.0
WD(1)=80.0
W(1)=WU(1)+WL(1)+WD(1)
WWMM=(1.0+B)*WM
!------------------------------------------------------------------------------------------------
DO I=1,24
EM(I)=K*E0(I)
!各层蓄水量变化，降水先补充上层，上层蓄满后再补充下层和深层，将降雨分配到各层
DU=WUM-WU(I)
DL=WLM-WL(I)
DD=WDM-WD(I)
IF(P(I)<=DU)THEN
WUX=P(I)+WU(I)
ELSE IF(P(I)>DU.AND.P(I)<=DU+DL)THEN
WUX=WUM
WLX=WL(I)+P(I)-DU
ELSE IF(P(I)>DU+DL.AND.P(I)<=DU+DL+DD)THEN
WUX=WUM
WLX=WLM
WDX=WD(I)+P(I)-DU-DL
ELSE
WUX=WUM
WLX=WLM
WDX=WDM
END IF
!----------------------------------------------------------------------
!三层蒸发模型计算蒸发,先蒸发上层再蒸发下层和深层，将蒸发能力分配到各层
IF(WUX>=EM(I))THEN
EU(I)=EM(I)
EL(I)=0.0
ED(I)=0.0
ELSE
EU(I)=WUX
IF(WLX>=C*WLM)THEN
EL(I)=(EM(I)-EU(I))*WLX/WLM
ED(I)=0.0
ELSE
IF(WLX>=C*(EM(I)-EU(I)))THEN
EL(I)=C*(EM(I)-EU(I))
ED(I)=0.0
ELSE
EL(I)=WLX
ED(I)=C*(EM(I)-EU(I))-EL(I)
END IF
END IF
END IF
WU(I+1)=WUX-EU(I)
WL(I+1)=WLX-EL(I)
WD(I+1)=WDX-ED(I)
E(I)=EU(I)+EL(I)+ED(I)
!---------------------------------------------------------
!产流计算
PE(I)=P(I)-E(I)
A(I)=WWMM*(1.0-(1.0-W(I)/WM)**(1.0/(1.0+B)))
IF(PE(I)<=0.0)THEN
R(I)=0
W(I+1)=W(I)+PE(I)
ELSE
IF(PE(I)+A(I)>WWMM)THEN
R(I)=PE(I)-(WM-W(I))
ELSE
R(I)= PE(I)-((WM-W(I))-WM*(1-(PE(I)+A(I))/WWMM)**(1.0+B)) !落一个*
END IF
W(I+1)=W(I)+P(I)-R(I)
END IF
!----------------------------------------------------------------------
!水源划分
SSM=SM*(1.0+EX)
AU(I)=SSM*(1.0-(1.0-S(I)/SM)**(1.0/(1.0+EX)))
IF(PE(I)>0)THEN
FR(I)=R(I)/PE(I)
IF(PE(I)+AU(I)<SSM)THEN
RS(I)=(PE(I)+S(I)-SM+SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*FR(I)
RSS(I)=(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*KSS*FR(I)
RG(I)=(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(EX+1.0))*KG*FR(I)
S(I+I)=(1.0-KSS-KG)*(SM-SM*(1.0-(PE(I)+AU(I))/SSM)**(1.0+EX))
ELSE
RS(I)=(PE(I)+S(I)-SM)*FR(I)
RSS(I)=SM*KSS*FR(I)
RG(I)=SM*KG*FR(I)
S(I+1)=(1.0-KSS-KG)*SM
END IF
ELSE
FR(I)=1.0-(1.0-W(I)/WM)**(B/(1.0+B))
RS(I)=0.0
RSS(I)=S(I)*KSS*FR(I)
RG(I)=S(I)*KG*FR(I)
S(I+1)=(1-KSS-KG)*S(I)
END IF
QS(I)=0.0;
QSS(I)=0.0;
QG(I)=0.0;
END DO
!------------------------------------------------------------------------------------
!汇流计算，地表径流采用瞬时单位线，壤中流和地下径流采用线型水库
DO I=1,24
DO J=1,3
QS(I+J-1)=QS(I+J-1)+RS(I)*UH(J)*F/(3.6*2.0)
END DO
IF(I==1)THEN
QSS(I)=QSS0*KKSS**(1/12)+RSS(I)*(1.0-KKSS**(1/12))*F/7.2
QG(I)=QGO*KKG**(1/12)+RG(I)*(1.0-KKG**(1/12))*F/7.2
ELSE
QSS(I)=QSS(I-1)*KKSS**(1/12)+RSS(I)*(1.0-KKSS**(1/12))*F/7.2
QG(I)=QG(I-1)*KKG**(1/12)+RG(I)*(1.0-KKG**(1/12))*F/7.2
END IF
Q(I)=QS(I)+QSS(I)+QG(I)
WRITE(*,*) Q(I)
END DO
!--------------------------------------------------------------------------------------
READ(*,*)
STOP
END PROGRAM XAJ