和歌山サイエンスカフェ「インフィニティ」

働き出して、水理学ってものが
あるのを知っても、なかなか良くわからなかったのですが、
昨年くらいから、腰を据えて勉強をはじめました。

初学者向けに、寒地土木研究所さんの
「現場の水理学」という貴重な資料があるのです。
http://river.ceri.go.jp/contents/tool/suirigaku.html

こちらに、不等流計算など初歩的なところから
のプログラムコードも掲載されています。
言語はfortranです。

でも、それをそのまま打ち込んでも、
動かなかったりしたので、
ほとんどそのままですが、
fortran90の勉強も兼ねて、自分で
ソースコードを作ってみました。

ということで、たまにそのコードを掲載させて
もらいたいと思います。
プログラミングにあたり、参考にした本は、

です。

演習問題３の１階のニュートン法によるもののソース

program newton
implicit none
integer i,j
integer,parameter ::nn=10
double precision q,n,ep,g
double precision aa,bb,cc
double precision fh,f1,dh,hhh
double precision d(nn),dl(nn),z(nn),b(nn),sh(nn),hh(nn)

d=(/ 0d0,500d0,1000d0,1200d0,1800d0,2100d0,2500d0,3000d0,3300d0,3800d0 /)
z=(/ 0d0,.5d0,.9d0,.8d0,2.0d0,2.3d0,3.0d0,3.0d0,3.5d0,4.0d0 /)
b=(/ 300d0,320d0,280d0,250d0,300d0,300d0,320d0,350d0,300d0,250d0 /)


do i=2,nn
dl(i)=d(i)-d(i-1)
enddo

q=1500d0
n=0.025d0
sh(1)=2.5d0
hh(1)=sh(1)+z(1)
ep=1d-3
g=9.8d0

do j=1,nn
if(j==1) then
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
endif

hhh=sh(j-1)

aa=q*q/(2*g*b(j)*b(j))
bb=-q*q*n*n*dl(j)/(2*b(j)*b(j))
cc=z(j)-(z(j-1)+sh(j-1)+q*q/(2*g*b(j-1)*b(j-1)*sh(j-1)*sh(j-1))+&
q*q*n*n*dl(j)/(2d0*b(j-1)*b(j-1)*sh(j-1)**(10d0/3d0)))

10fh=hhh+aa/(hhh**2d0)+bb/(hhh**(10d0/3d0))+cc
f1=1-2d0*aa/(hhh**3d0)-10d0*bb/(3d0*hhh**(13d0/3d0))
dh=-fh/f1
if(abs(fh) < ep ) then
sh(j)=hhh
hh(j)=sh(j)+z(j)
write(*,*) j,d(j),dl(j),b(j),z(j),sh(j),hh(j)
cycle
else
hhh=hhh+dh
goto 10
endif

enddo

end program newton