http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html
ある場所（緯度経度）での変位を出力しようと思ったら、結構頭がこんがらがった。というわけで、wrapperを作ったので、プログラムをメモ代わりにアップしておく。断層パラメータの記載は一般的な定義のままで、位置は断層上端左端にした。ああ、ややこしかった。
注意としては、深さは正の値を入れること、長さLの場合はAL1=0, AL2=Lとすること、幅Wの場合は、AW1=-W, AW2=0とすること。
一般的な断層パラメータの定義が数式と少し違い（走向は北から反時計回りとか）、注意が必要ですね。すぐ忘れそうなので、やはりwrapperを作っておいた方が良いな。もう少しモジュールがたまってきたら、arでライブラリにしてしまいたい。
ちなみに、こんな感じで呼び出す。
test.f90:
program test
use okada92_util
implicit none
type(disp) :: du
type(flt) :: f
real(8) :: lat, lon, h
integer :: i, j
f%lambda = 1.d0 ! ずるしている
f%mu = 1.d0
f%lat = 32.7725d0
f%lon = 136.0282d0
f%depth = 5000.d0
f%str = 261.766d0
f%dip = 8.903d0
f%rake = 90.d0
!f%rake = 0.d0
f%length = 100000.d0
f%width = 25263.533d0
f%slip = 5.d0
do i = -50, 50
do j = -50, 50
lat = f%lat + dble(i)*0.1d0
lon = f%lon + dble(j)*0.1d0
h = 0.d0
du = flt_disp(f, lat, lon, h)
write(*,'(F7.3,",",F7.3,3(",",F15.7))') lon, lat, du%e, du%n, du%u
enddo
enddo
stop
end program test
モジュールは下記のようなもの。
okada92_util.f90:
module okada92_util
type disp
real(8) :: e, n, u
end type disp
type flt
!real(8) :: alpha ! (lambda+mu)/(lambda+2*mu) !αをインプットにした方が良かったかも
real(8) :: lambda, mu ! (lambda+mu)/(lambda+2*mu)
real(8) :: lat, lon ! latitude, longitude (deg) at the corner of upper edge
real(8) :: depth ! fault depth (positive is ...)
real(8) :: str ! strike angle (deg)
real(8) :: dip ! dip angle (deg)
real(8) :: rake ! rake angle (deg)
real(8) :: length, width ! fault length, width (m)
real(8) :: slip ! slip (m)
end type flt
! ellipsoid
type ellip
real(8) :: a, b, e2, f, fi
end type ellip
! parameters
real(8), parameter :: PI = 3.1415926535897932d0
real(8), parameter :: D2R = PI/180.d0
real(8), parameter :: R2D = 180.d0/PI
contains
! ----------------------------------------------------------
type(disp) function flt_disp(f, lat, lon, h)
! ----------------------------------------------------------
! Calculates displacement at (lat, lon, h) due to a
! rectangular fault. The equations are given by
! Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040].
!
! This function is a wrapper of 'DC3D' on the website:
! http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html
!
! Input:
! f : type(flt), includes
! f%lambda, mu ! (lambda+mu)/(lambda+2*mu)
! f%lat, lon ! latitude, longitude (deg) at the corner of upper edge.
! f%depth ! fault depth (positive is downward)
! f%str ! strike angle (deg)
! f%dip ! dip angle (deg)
! f%rake ! rake angle (deg)
! f%length, width ! fault length, width (m)
! f%slip ! slip (m)
! lat, lon, h : latitude(deg), longitude(deg), h(m)
! positive h is upward.
! Output:
! disp : type(disp),
! disp%e ! eastward disp. (m)
! disp%n ! northward disp. (m)
! disp%u ! vertical disp. (m)
!
! Created:
! 2013/04/29 Sako
! ----------------------------------------------------------
type(flt), intent(in) :: f
real(8), intent(in) :: lat, lon, h
! inputs of DC3D
real(8) :: a, x, y, z, d, dip
real(8) :: al1, al2, aw1, aw2
real(8) :: disl1, disl2, disl3
! output of DC3D
real(8) :: ux, uy, uz
real(8) :: uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz
integer :: iret
! tmp
real(8) :: dx, dy, t
!real(8) :: alpha ! (lambda+mu)/(lambda+2*mu)
a = (f%lambda+f%mu)/(f%lambda+2.d0*f%mu)
d = f%depth ! 断層上端の深さ。深さは正の値にする。
dip = f%dip
al1 = 0.d0
al2 = f%length ! 断層の長さ。
aw1 = -f%width ! 断層の幅。図を良く見るとAW1はマイナスということで結構紛らわしい。
aw2 = 0.d0
disl1 = f%slip * dcos(f%rake*D2R) ! すべり量を各軸に分解。
disl2 = f%slip * dsin(f%rake*D2R)
disl3 = 0.d0
! 目的の緯度経度をローカル座標系のx, yに変換する。走向の定義が北からなので、ここもややこしい。
call bl2xy(f%lat, f%lon, lat, lon, dx, dy)
t = (-f%str+90.d0)*D2R; ! rotate to local coordinate
x = dx*dcos(t) + dy*dsin(t) ! t = -str+90
y =-dx*dsin(t) + dy*dcos(t)
z =-h
! SUBROUTINE DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
!* AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
!* UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET)
call DC3D(a, x, y, z, d, dip, al1, al2, aw1, aw2, disl1, disl2, disl3, &
ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz, iret)
if (iret == 0) then
flt_disp%e = ux*dsin(f%str*D2R) - uy*dcos(f%str*D2R) ! 忘れずに座標変換。走向の向きは北からだからこうなる。
flt_disp%n = ux*dcos(f%str*D2R) + uy*dsin(f%str*D2R)
flt_disp%u = uz
else
! singular location
flt_disp%e = 0.d0
flt_disp%n = 0.d0
flt_disp%u = 0.d0
endif
return
end function flt_disp
! ----------------------------------------------------------
subroutine bl2xy(lat0, lon0, lat, lon, dx, dy)
! ----------------------------------------------------------
! Calculates distance between reference station(lat0, lon0) and
! station2(lat, lon).
!
! Reference: Ohno S.(1987), Geodetic Methods (in Japanese).
! see equations [1.59] and [1.62].
!
! Input:
! lat0, lon0 : position of station1 (reference)
! lat, lon : position of station2
!
! Output:
! dx, dy : distance between stations (sta1 -> sta2).
!
! Created:
! 2013/04/29 Sako
! ----------------------------------------------------------
implicit none
real(8), intent(in) :: lat0, lon0, lat, lon
real(8), intent(out) :: dx, dy
type(ellip) :: ell
real(8) :: dlat, dlon, Bm
real(8) :: a, fi, f, b, e2
! WGS84
ell = ellipsoid('WGS84')
a = ell%a
fi= ell%fi
f = ell%f
b = ell%b
e2= ell%e2
dlat = (lat - lat0)*D2R
dlon = (lon - lon0)*D2R
Bm = dsin(dlat)
dx = a*dcos(lat0*D2R) / dsqrt(1.d0-e2*dsin(lat0*D2R)**2.d0)*dlon
dy = a*(1.d0-e2)*dlat*(1.d0+3.d0/4.d0*e2-3.d0/4.d0*e2*dcos(2.d0*Bm)+1.d0/8.d0*e2*(dlat)**2.d0*dcos(2.d0*Bm))
return
end subroutine bl2xy
! ----------------------------------------------------------
type(ellip) function ellipsoid(ellip_name)
! ----------------------------------------------------------
implicit none
character(len=*), intent(in), optional :: ellip_name
real(8) :: a, fi, f
! default is GRS80
a = 6378137.d0
fi= 298.257222101d0
if (present(ellip_name)) then
if (ellip_name == 'GRS80') then
a = 6378137.d0
fi= 298.257222101d0
elseif (ellip_name == 'WGS84') then
a = 6378137.d0;
fi= 298.257223563d0;
else
! ...
endif
endif
f = 1.d0/fi
ellipsoid%a = a
ellipsoid%fi= fi
ellipsoid%f = f
ellipsoid%b = a*(1.d0-f)
ellipsoid%e2= 1.d0 - (1.d0-f)**2.d0
return
end function ellipsoid
end module okada92_util
0 件のコメント:
コメントを投稿