（４）フォートラン・プログラム・リスト
c*****************************************
c2 段タンク型貯留関数モデルパラメータ同定
c*****************************************cUSING str2t < data‑file > out‑filecstorage‑discharge relationshipc[tank1]
cs1 = k11*q1**p1 + k12*d(q1**p2)/dt
cds1/dt=r‑q1‑f1
c
f1=k13*q1
c[tank2]
cs2 = k21*q2 + k22*dq2/dt
cds2/dt=f1‑q2c
optimization of parameters c11 and c12 and c13
cc11=k11/A**0.24 c12=k12/k11**2/Rave**‑0.2468cc13=1+k13cp1 = 0.6 and p2 = 0.4648
c‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑ccharacter title1*128,title2*128,title3*128,title*128& ,gp*128,gp1*80,gp2*80cdimension e(200),e1(200),e2(200),dpa(9)dimension pas(200,9),bb(200,200),qa(200)ccommon/title/ title1,title2,title3
common/data/ nl,tcommon/data1/ nq,area,r(200),qr(200),q(200)common/data2/qt(200),qq(200),b(200),q1̲1(200),q1̲2(200)cc‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑
ccncase:計算洪水数
read(5,'(i5)') ncase
write(6,'(i5)') ncase
do 9999 kk=1,ncasecタイトル title1：水系名 title2：河川名
ctitle3：観測所名 title：洪水名
read (5,'(a80)') title1
read (5,'(a80)') title2
read (5,'(a80)') title3
read (5,'(a80)') title
carea：流域面積,ipas：計算の実行判別(1:計算する 2:計算しない)ciy,im1,id1,ih1：計算開始年月日時 im2,id2,ih2：計算終了月日時
read (5,'(f8.0,8i5)') area,ipas,iy,im1,id1,ih1,im2,id2,ih2cn:微分方程式の階数 m:モデル定数の個数 nl:分割数
ckount:繰り返し計算数（最大）p:収束条件 rb:平均雨量強度
read(5,4) n,m,nl,kount,p,rb
4 format(4i5,2f8.0)c
計算時間間隔 t179

cccc
c
c
c

t=1.0/float(nl)パラメータ c ,δ,tc 値の入力
read (5,'(5f8.0)') cc1,cc2,c3,delta,tc
nr:雨量データ数 nq:流量データ数
read (5,'(2i5)') nr,nq
実測雨量データ入力
do i=1,nqr(i)=0.0
end doread (5,'(10f8.0)') (r(i),i=1,nr)
実測流量データ入力
read (5,'(10f8.0)') (qr(i),i=1,nq)
地下水分離方法
nal=0 ‑> 地下水が実測を超えたら地下水=実測とする。
nal=1 ‑> 重み係数αを使う。
nal=0

cif(ipas.eq.0) go to 9999
cwrite(6,'(i5)') kkwrite(6,'(a80)') title1
write(6,'(a80)') title2
write(6,'(a80)') title3
write(6,'(a80)') titlecdo i=1,80
if(title2(i:i).eq.' ') go to 77
end do77 iz1=ido i=1,80
if(title3(i:i).eq.' ') go to 78
end do78 iz2=iwrite(gp1,'(a2,i2,a4,i2,a4)') '(a',iz1,',a,a',iz2,',a2)'
write(gp2,gp1) title2,'(',title3,' ,'iz3=iz1+iz2+3
write(gp1,'(a2,i2,a5,a)') '(a',iz3,',a128',')'
cwrite(gp,'(''流域面積：'',f8.2,''km2)'',& i4,''年'',i2,''月'',i2,''日〜'',i2,''月'',i2,''日洪水'')')& area,iy,im1,id1,im2,id2
write(6,gp1) gp2,gp
cwrite(6,'(i10)') nqwrite(6,210)
210 format('流域面積',' 雨量強度',' c11 初期値',' c12 初期値',& ' c13 初期値','
δ','Tc')write(6,'(2f10.2,3f10.3,f10.2,f10.4)') area,rb,cc1,cc2,c3,delta,Tccc

m1=m+1
緩和係数

180fac=0.0
p1,p2 値
p1=0.6p2=0.4648
cz1,z2 基準化パラメータ値
z1=1.0z2=1.0c‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑
c00=(delta/tc)**2
c01=delta**2/tc
c‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑c
地下水流出分離call chika(c00,c01,nal)
ccc00,c01:初期値
cnal:地下水分離方法選択
c‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑
c
流量‑‑>流出高
nnn=nqdo i=1,nq
qt(i)=3.6*qt(i)/area
if ( qt(i).le.0.00001 ) then! 中間流出が 0 になる時間
nnn=i‑1
go to 123
end ifend do123 continuecdo 300 kkll=1,kountc
wk1=cc1*z1*area**(0.24)
wk2=cc2*z2*(wk1**2)*rb**(‑0.2648)cc1=(wk1/wk2)*(p1/p2)
c2=(1./wk2)
p3=p1/p2‑1
p4=1./p2cu1=0.0u2=0.0c31=0.0
u3=0.0u4=0.0c32=0.0cc
初期流出高
q1=qt(1)**p2
qb=0.0cq2=0.0
sumuu=0.0c181

sumu1=0.0
sumu2=0.0
sume=0.0
sume1=0.0
sume2=0.0
sotais=0.0c
do 100 i=1,nq
tt=ic
do 200 j=1,nlcif ( q1.eq.0 ) thena1=0.0else
a1=‑c1*p3*q1**(p3‑1)*q2‑c2*c3*p4*q1**(p4‑1)
end ifa2=‑c1*q1**p3
a3=a1+a2**2
a4=a1+a3x=c2*c3*q1**p4*(p4‑1)+c1*p3*q1**p3*q2+c2*r(i)ccall phigam ( a1,a2,a3,a4,phi1,phi2,phi3,phi4,&
gam1,gam2,gam3,gam4 )cd1=‑(1./wk2)*(p1/p2)*q1**p3*q2
d2=(wk1/(wk2**2.))*(p1/p2)*q1**p3*q2+&(1./wk2**2.)*q1**p4*c3‑r(i)/(wk2**2.)d3=(‑1./wk2)*q1**p4c
u1=phi1*u1+phi2*u3+gam2*d1
u2=phi1*u2+phi2*u4+gam2*d2
c31=phi1*c31+phi2*c32+gam2*d3
u3=phi3*u1+phi4*u3+gam4*d1
u4=phi3*u2+phi4*u4+gam4*d2
c32=phi3*c31+phi4*c32+gam4*d3
cq1=phi1*q1+phi2*q2+gam2*xq2=phi3*q1+phi4*q2+gam4*x
cif ( q1.lt.0 ) q1=0
cbb(i,j)=(c3‑1)*q1**(1./p2)
c200 continuecv1=u1*((1./p2)*q1**(1./p2‑1.))*cc1*area**(0.24)
v2=u2*((1./p2)*q1**(1./p2‑1.))*cc2*(wk1**2.)&*rb**(‑0.2648)
v3=c31*((1./p2)*q1**(1./p2‑1.))
qq(i)=q1**(1./p2)

182

sumuu=sumuu+v1*v2
sumu1=sumu1+v1**2.
sumu2=sumu2+v2**2.
e(i)=qt(i)‑qq(i)sume=sume+e(i)e1(i)=e(i)*v1
e2(i)=e(i)*v2
cc

成分回帰分析パラメータ
if(i.le.nnn) thenpas(i,1)=v1
pas(i,2)=v2
pas(i,3)=v3
pas(i,4)=e(i)end ifc
sume1=sume1+e1(i)sume2=sume2+e2(i)c
100 continuecc

成分回帰分析
call momreg(9,m1,200,nnn,pas,dpa)cdz1=dpa(1)dz2=dpa(2)dc3=dpa(3)cg1=abs(dz1/z1)
g2=abs(dz2/z2)
g3=abs(dc3/c3)cif ( g1.lt.0.001.and.g2.lt.0.001.and.g3.lt.0.001 ) go to 999cfac=(1+fac)*0.5
dz1=dz1*fac
dz2=dz2*fac
dc3=dc3*fac
c11=cc1*z1
c12=cc2*z2
z1=z1+dz1
z2=z2+dz2
c3=c3+dc3
cif ( z1.lt.0.001 ) z1=z1‑dz1
if ( z2.lt.0.001 ) z2=z2‑dz2if ( c3.lt.0.001 ) c3=c3‑dc3c
300 continue
999 continue
c183

c

表面・中間流出が 0 になった時間以降を流量 0 とする。
if ( nnn.ne.nq ) thendo i=nnn+1,nqqt(i)=0.0
end doend ifcdo 333 i=1,nq
qq(i)=qq(i)*area/3.6
qt(i)=qt(i)*area/3.6
do j=1,nlbb(i,j)=bb(i,j)*area/3.6
sigb=sigb+bb(i,j)
end doif ( qq(i).lt.0 ) qq(i)=0
sigr=sigr+r(i)sigqt=sigqt+qt(i)sigqq=sigqq+qq(i)333 continuecwk4=(c3‑1.0)/c00
wk3=c01*wk4
cc

地下水流出再計算
c0=1.0/wk4
c1=c01ccall cal2(c0,c1,nq,bb,qa)cdo i=1,nqq1̲2(i)=qa(i)q(i)=qq(i)+qa(i)end doccc
ck3,k4:パラメータbb:損失流量
do i=1,nqb(i)=bb(i,nl)*3.6/area
end docrtot=0.0
qrtot=0.0
qtot=0.0
btot=0.0
qrmax=‑999.99
qmax=‑999.99
sotai=0.0
do i=1,nq
if(qrmax.lt.qr(i)) qrmax=qr(i)if(qmax.lt.q(i)) qmax=q(i)184

rtot=rtot+r(i)qrtot=qrtot+qr(i)qtot=qtot+q(i)btot=btot+b(i)sotai=sotai+abs(qr(i)‑q(i))/qr(i)end doqrtot=qrtot*3.6/area
qtot=qtot*3.6/area
hiryu=qrmax/areasotai=sotai/float(nq)peakg=abs(qrmax‑qmax)/qrmax
cwrite(6,211)
211 format(7x,'c11',7x,'c12',7x,'c13',7x,'k21',7x,'k22',& 7x,'Jpe',7x,'Jre',' 収束回数')
write(6,'(7f10.3,i10)') c11,c12,c3,wk3,wk4,peakg,sotai,kkll
write(6,212)
212 format('比流量
総雨量 観測総流出 計算総流出',&' 観測ピーク 計算ピーク 総浸透')
write(6,'(7f10.3)') hiryu,rtot,qrtot,qtot,qrmax,qmax,btot
write(6,'('' 時間'',4x,''雨量'',2x''実測流量'',2x,''計算流量'',&2x,''直接流出'',2x,''計算直接流出'',2x,''浸透供給量'',&2x,''地下水流出'',2x,''計算地下水流出'')')
do i=1,nq
write(6,'(i5,f8.4,3f10.4,f14.4,2f12.4,f16.4)')&i,r(i),qr(i),q(i),qt(i),qq(i),b(i),q1̲1(i),q1̲2(i)end doc9999 continue
cstop
endccccsubroutine chika(c0,c1,n)** 地下水流出成分の分離 **dimension qrr(200,200),qa(200)common/data/ nl,tcommon/data1/ nq,area,r(200),qr(200),q(200)common/data2/qt(200),qq(200),b(200),q1̲1(200),q1̲2(200)cxnl=float(nl)cdo i=1,nq
do j=1,nl
x=x+1.0
if ( i.eq.1) thenqrr(i,j)=qr(i)elseqrr(i,j)=qr(i‑1)+x*(qr(i)‑qr(i‑1))/xnl

185

end ifend dox=0.0
end doccall cal2(c0,c1,nq,qrr,qa)
cif ( n.eq.0 ) thendo i=1,nq
if ( qr(i).le.qa(i) ) qa(i)=qr(i)q1̲1(i)=qa(i) ! output データ
qt(i)=qr(i)‑qa(i)end doelse
a=1.0
do i=1,nq
aa=qr(i)/qa(i)if ( aa.le.a ) a=aa
end dowrite(6,*)
write(6,'(''重み係数α='',f12.8)') a
do i=1,nqq1̲1(i)=a*qa(i)qt(i)=qr(i)‑a*qa(i)end doend ifreturnendcc
subroutine cal2(c0,c1,nq,bb,q1)cdimension bb(200,200),q1(200)common/data/ nl,tc
a1=‑c0a2=‑c1
a3=a1+a2**2
a4=a1+a3
call phigam(a1,a2,a3,a4,phi1,phi2,phi3,phi4,&
gam1,gam2,gam3,gam4 )g1=0.0g2=0.0do 10 i=1,nq
do 20 j=1,nl
x=c0*bb(i,j)
g1=phi1*g1+phi2*g2+gam2*x
g2=phi3*g1+phi4*g2+gam4*x
20 continue
q1(i)=g1
10 continue

186return
endcc&subroutine phigam(a1,a2,a3,a4,phi1,phi2,phi3,phi4,gam1,gam2,gam3,gam4 )common/data/ nl,tc
phi1=1+(a1*t**2)/2+(a1*a2*t**3)/6+(a1*a3*t**4)/24
phi2=t+(a2*t**2)/2+(a3*t**3)/6+(a2*a4*t**4)/24
phi3=a1*t+(a1*a2*t**2)/2+(a1*a3*t**3)/6+(a1*a2*a4*t**4)/24
phi4=1+a2*t+(a3*t**2)/2+(a2*a4*t**3)/6&+((a1*a3+a2**2*a4)*t**4)/24
gam1=t*(1+(a1*t**2)/6+(a1*a2*t**3)/24)
gam2=t**2*(0.5+(a2*t)/6+(a3*t**2)/24)
gam3=t**2*a1*(0.5+(a2*t)/6+(a3*t**2)/24)
gam4=t*(1.0+(a2*t)/2+(a3*t**2)/6+(a2*a4*t**3)/24)
returnendcc
subroutine momreg(n1,n,m1,md,x,dpa)
component regression method
compute the correction terms of parameters (dpa)
dimension x(m1,n1),dpa(9),cov(9,9),g(200,9),y(200,9)
dimension binv(9,9),coe(9),st(9)na=n‑1c**compute covariance matrix
call sqcov(n1,na,m1,md,x,cov)
do 50 i=1,na
50 st(i)=sqrt(cov(i,i))
do 52 i=1,na
s=st(i)do 52 j=1,i
s1=st(j)cov(i,j)=cov(i,j)/(s*s1)
52 cov(j,i)=cov(i,j)
c**factorization of cov(i,j) by lower triangularc**cholesky method (cov = l * u)c**l = lower triangular u = upper triangularc**compute the inverse of u(i,j)
c**call lowtri(n1,na,cov,binv)
do 54 j=1,na
s=st(j)do 54 i=1,md
54 y(i,j)=x(i,j)/s
do 20 i=1,md
do 20 j=1,na
s=0.
do 22 k=1,jc**c**18722 s=s+y(i,k)*binv(k,j)
20 g(i,j)=s
do 24 i=1,na
s=0.0
do 26 j=1,md
26 s=s+g(j,i)*x(j,n)
24 coe(i)=s
do 30 i=1,na
s=0.
do 29 j=i,na
29 s=s+binv(i,j)*coe(j)30 dpa(i)=s/st(i)return
end
cc
c**
c**
c**
c**
c**
c**
c**
c**

subroutine lowtri(n1,n,p,binv)
lower triangular cholesky factorizationp = u*bp = symmetric matrix
u = lower triangular matrixb = upper triangular matrix (b = ut)
binv = inverse matrix of b
= upper triangular matrix

dimension p(n1,n1),binv(n1,n1)
dimension u(9,9),b(9,9)
c**compute lower triangular u(i,j)
do 5 j=1,n‑1u(j,j)=abs(p(j,j))u(j,j)=sqrt(u(j,j))al=1./u(j,j)
do 5 k=n,j+1,‑1u(k,j)=al*p(k,j)
be=u(k,j)
do 5 i=k,n5 p(i,k)=p(i,k)‑u(i,j)*beu(n,n)=abs(p(n,n))u(n,n)=sqrt(u(n,n))c**c**b = transpose of u
do 40 i=1,ndo 40 j=i,n40 b(i,j)=u(j,i)
c**c**compute inverse if b(i,j)
binv(1,1)=1./b(1,1)
do 50 j=2,nbinv(j,j)=1./b(j,j)
jm1=j‑1
do 50 k=1,jm1

188sum=0.0
do 52 i=k,jm1
52 sum = sum ‑ binv(k,i)*b(i,j)
50 binv(k,j)=sum*binv(j,j)
returnend
ccc**

subroutine sqcov(n1,n,m1,md,x,cov)
compute covariance matrix
dimension x(m1,n1),cov(n1,n1)
do 10 i=1,ndo 10 j=1,is=0.
do 12 k=1,md
12 s=s+x(k,i)*x(k,j)
10 cov(i,j)=s
returnend

189