目的
JAGURSで津波のシミュレーションをするために、誰かが提案した津波波源モデルをJAGURSで使えるようにする。
前回の続き
波源モデル
独立行政法人 建築研究所 国際地震工学センターが公表している、藤井雄士郎（建築研究所 国際地震工学センター），佐竹健治（東京大学 地震研究所）らによる2011年3月11日東北地方太平洋沖地震の津波波源モデル
(Ver. 6.2, Ver. 7.0 and Ver. 8.0) のVar８を使います。
（そういや、どこぞの論文にここの断層モデルをがっつり使ったけど、コピー送ってない著者の存在を知ってる）
dispCDF55_ver8.0.in 動いたっぽい断層？のパラメータ
snap_55_ver8.0.xls 分割された断層が５分間でどのように動いたか
断層のモデルの中身はさっぱり意味不明ですが、かろうじて緯度経度だけわかります。
測地系について詳しいことは書いてないっぽい気がしますが、旧測地だろうが世界測地だろうが、誤差は1kmも無いので多分大丈夫でしょう！（適当）
JAGURSの入力形式
断層モデルの入力には二通りの方法があります。
動いだ断層を地形ファイルとして入力する方法と、位置や方向などを数字で入力する方法です。
今回は、２つ目のすうじで入力する方法をとります。
JAGURSの英語版のマニュアルを見ると、
だそうです。
というわけで、３０秒ごとの断層の位置などなどをが書き込まれたファイルを作らなければなりません。
変換プログラム
安定のFortran
緯度経度を謎UTM座標に変換して並べ替えと、秒ごとのファイルにする。
fault001.f
implicit none
real(kind=8),dimension(55) ::length,width,depth,strike
real(kind=8),dimension(55) ::dip,rake,slip,lat,lon
real(kind=8),dimension(55,11):: fault
integer i,j,k
integer iosnum
character(len=200) temp,fname
real(kind=8) :: x,y
!lengthwidthdepth strikediprakesliplatlong
open(21,file="dispCDF55_ver8.0.txt",status='old')
read(21,*)temp
do i=1,55
read(21,*)length(i),width(i),depth(i),
c strike(i),dip(i),rake(i),slip(i),lat(i),lon(i),temp
end do
close(21)
open(22,file='snap_55_ver8.0-2.txt',status='old')
do i=1,55
read(22,*) fault(i,:)
end do
close(22)
do k=1,11
write(fname,'("fault_",i3.3,".txt")')(k-1)*30
write(*,*)trim(fname)
do i=1,55
open(31,file=trim(fname),status='unknown')
!x,y,b,l
call wgs2utm_nazo(y,x,lat(i),lon(i))
write(31,'(9F15.7)') y,x,depth(k),length(k),width(k),dip(k),
c strike(k),rake(k),fault(i,k)
end do
end do
close(31)
stop
end
!**************************************************
subroutine utm2wgs_nazo(b,l,x,y,zn)
implicit none
real(kind=8),intent(out):: b,l
real(kind=8),intent(in)::x,y
integer,intent(in) :: zn
real(kind=8) :: a,f,e2,ed2,e4,e6,e8
real(kind=8) :: cad,cbd,ccd,cdd,ced
real(kind=8) :: l0,eps,m0,cb,eta2,rn
real(kind=8) :: xd,yd,tx,ty
real(kind=8) :: db,dl
real(kind=8) :: rd
rd =3.14159265358979323846264338327950288419716939937d0/180.d0
a = 6378137.d0! # 赤道半径 (WGS 84)
f = 1.d0 / 298.257223d0 ! # 扁平率 (WGS 84)
e2 = 2.d0*f - f*f !; # 第1離心率
ed2 = e2 / (1.d0-e2)!; # 第2離心率
e4 = e2 * e2
e6 = e4 * e2
e8 = e6 * e2
cad=a*(1.d0-e2)*(1.d0+3.d0/4.d0*e2+45.d0/64.d0*e4
& +175.d0/256.d0*e6+11025.d0/16384.d0*e8)
cbd=a*(1.d0-e2)*( 3.d0/4.d0*e2+15.d0/16.d0*e4
& +525.d0/512.d0*e6+ 2205.d0/2048.d0*e8)
ccd=a*(1.d0-e2)*( 15.d0/64.d0*e4
& +105.d0/256.d0*e6+ 2205.d0/4096.d0*e8)
cdd=a*(1.d0-e2)*(
& 35.d0/512.d0*e6+ 315.d0/2048.d0*e8)
ced=a*(1.d0-e2)*(
& 315.d0/16384.d0*e8)
l0 =143.d0 !dble(zn)*6.d0-183.d0
eps = 0.0000001d0/3600.d0*rd
m0 = 0.9996d0
ty = y -500000.d0
tx = x/m0
ty = ty/m0
db = 1.d0
dl = 1.d0
b = tx/cad
cb = dcos(b)
eta2 = ed2*cb*cb
rn = a/(1.d0-f)/dsqrt(1.d0+eta2)
l=ty/rn/cb
do while(db>eps .or. dl>eps)
call ll2xy(xd,yd,b,l,a,f,e2)
db = (tx-xd)/cad
cb = dcos(b)
eta2=ed2*cb*cb
rn = a/(1.d0-f)/dsqrt(1.d0+eta2)
dl = (ty-yd)/rn/cb
b = b+db
l = l+dl
end do
b = b/rd
l = l/rd+l0
end subroutine utm2wgs_nazo
!****************************************************
!***********************************************
subroutine wgs2utm_nazo(x,y,b,l)
implicit none
real(kind=8) b,l
real(kind=8) x,y
real(kind=8) rd
integer zn
real(kind=8) l0
real(kind=8) a,f,e2
rd =3.14159265358979323846264338327950288419716939937d0/180.d0
a=6378137.d0
f=1.0d0/298.257223d0
e2=2.0d0*f-f*f
zn = mod(int((l+180.d0)/6.d0),60)+1
l0= 143.d0!dble(zn)*6.d0-183.d0 !原点だそうだ。
!write(*,*) '中央罫線：東経',l0
call ll2xy(x,y,b*rd,(l-l0)*rd,a,f,e2)
x = x*0.9996d0
y = y*0.9996d0
y = y+500000.d0
return
end
!***********************************************
!***********************************************
subroutine ll2xy(x,y,b,dl,a,f,e2)
implicit none
real(kind=8) x,y
real(kind=8) b,dl
real(kind=8) a,f,e2
real(kind=8) dl2,dl4
real(kind=8) sb,cb,cb3,cb5,tb2,tb4,ed2,eta2,rn
real(kind=8) lom
x=0.d0
y=0.d0
dl2=dl**2
dl4=dl**4
sb=dsin(b)
cb=dcos(b)
cb3 = cb*cb*cb
cb5 = cb3*cb*cb
tb2 = sb*sb/cb/cb
tb4 = tb2*tb2
ed2=e2/(1.d0-e2)
eta2 = ed2*cb*cb
rn=a/(1.d0-f)/dsqrt(1.d0+eta2)
x=lom(b,a,f,e2)
!write(*,*) '赤道からの距離:',nint(x)/1000.d0
x =x+ rn/2.d0*sb*cb*dl2
x =x+ rn/24.d0*sb*cb3*(5.d0-tb2+9.d0*eta2+4.d0*eta2*eta2)*dl4
x =x+ rn/720.d0*sb*cb5*(61.d0-58.d0*tb2+tb4)*dl4*dl2
y = rn*cb*dl
y =y+ rn/6.d0*cb3*(1.d0-tb2+eta2)*dl*dl2
y =y+ rn/120.d0*cb5*(5.d0-18.d0*tb2+tb4)*dl*dl4
return
end
!*****************************************************
!*****************************************************
function lom(b,a,f,e2)
implicit none
real(kind=8) lom
real(kind=8) b,a,f,e2
real(kind=8) e4,e6,e8
real(kind=8) cad,cbd,ccd,cdd,ced
e4=e2*e2
e6=e4*e2
e8=e6*e2
cad=a*(1.d0-e2)*(1.d0+3.d0/4.d0*e2+45.d0/64.d0*e4
& +175.d0/256.d0*e6+11025.d0/16384.d0*e8)
cbd=a*(1.d0-e2)*( 3.d0/4.d0*e2+15.d0/16.d0*e4
& +525.d0/512.d0*e6+ 2205.d0/2048.d0*e8)
ccd=a*(1.d0-e2)*( 15.d0/64.d0*e4
& +105.d0/256.d0*e6+ 2205.d0/4096.d0*e8)
cdd=a*(1.d0-e2)*(
& 35.d0/512.d0*e6+ 315.d0/2048.d0*e8)
ced=a*(1.d0-e2)*(
& 315.d0/16384.d0*e8)
lom= cad*b-cbd/2.d0*dsin(2.d0*b)+ccd/4.d0*dsin(4.d0*b)
& -cdd/6.d0*dsin(6.d0*b)+ced/8.d0*dsin(8.d0*b)
return
end
!***************************************************
!*****************************************************
subroutine llh2xyz(x,y,z,b,l,h,a,e2)
implicit none
real(kind=8) x,y,z
real(kind=8) b,l,h,a,e2
real(kind=8) br,lr,sb,cb,rn
real(kind=8) rd
rd =3.14159265358979323846264338327950288419716939937d0/180.d0
br = b*rd
lr = l*rd
sb = dsin(br)
cb = dcos(br)
rn = a/dsqrt(1.d0-e2*sb*sb)
x = (rn+h) * cb * dcos(lr)
y = (rn+h) * cb * dsin(lr)
z = (rn*(1.d0-e2)+h) * sb
!write(*,*)e2
end
!***************************************************
!*************************************************
subroutine xyz2llh(b,l,h,x,y,z,a,e2)
implicit none
real(kind=8) x,y,z,a,e2
real(kind=8) b,l,h
real(kind=8) bda,p,t,st,ct,sb,rn
real(kind=8) rd
rd =3.14159265358979323846264338327950288419716939937d0/180.d0
bda = dsqrt(1.d0-e2)
p = dsqrt(x**2+y**2)
t = datan2(z, p*bda)
st = dsin(t)
ct = dcos(t)
b = datan2(z+e2*a/bda*st**3, p-e2*a*ct**3)
l = datan2(y, x)
sb = dsin(b)
rn = a / dsqrt(1.d0-e2*sb*sb)
h = p/dcos(b) - rn
b = b/rd
l = l/rd
end
!******************************************************
このプログラムで、fault_000.txt ~ fault_300.txt が作られます。
dispCDF55_ver8.0.txtのslipは全移動量なので、snap_55_ver8.0-2.txtの各秒ごとの値に入れ替える必要があります。（これに気づくのに、かつての僕は半日悩んだ）
このファイルをまとめたリストをかっこよくコマンドで作ります。
find ./ | grep fault_ | sort > faultlist.txt
まとめ
3.11の公開されている波源モデルを使ってJAGURSに入力するファイルを作った。