津波波源モデルをJagursで使えるようにする

##目的
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に入力するファイルを作った。