[Fortran] 纯文本查看 复制代码
Module SACParameter
Integer, parameter :: FileUnit = 100
Character( Len = * ), Parameter :: FileName = "..\WAV\WAVFormYiliang\yiliang\AST.GZ.2012251031817.00.BHE"
Type :: EventPara
Integer( kind = 4 ) :: StationNumber = 1
Integer( kind = 4 ) :: FileNumber
Character( Len = 4 ) :: StationName
Character( Len = 4 ) :: FileName
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Real(kind = 4), allocatable :: StationDELTA(:)
End Type EventPara
Type :: SACHeader
Real( Kind = 4 ), Allocatable :: DELTA, DEPMIN, DEPMAX, SCALE, ODELTA, B, E, &
O, A, INTERNAL1, T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, F, &
RESP0,RESP1,RESP2,RESP3,RESP4,RESP5,RESP6,RESP7,RESP8,RESP9, &
STLA,STLO,STEL,STDP,EVLA,EVLO,EVEL,EVDP,MAG, &
USER0, USER1, USER2, USER3, USER4, USER5, USER6, USER7, USER8, USER9, &
DIST, AZ, BAZ, GCARC, INTERNAL2, INTERNAL3, &
DEPMEN, CMPAZ, CMPINC, XMINIMUM, XMAXIMUM, YMINIMUM, YMAXIMUM, &
ADJTM, UNUSED1, UNUSED2, UNUSED3, UNUSED4, UNUSED5, UNUSED6
Integer( Kind = 4 ), Allocatable :: NZYEAR, NZJDAY, NZHOUR, NZMIN, NZSEC, NZMSEC, NVHDR, &
NORID, NEVID, NPTS, NSPTS, NWFID, NXSIZE, NYSIZE, UNUSED7, &
IFTYPE, IDEP, IZTYPE, UNUSED8, IINST, ISTREG, IEVREG, IEVTYP, IQUAL, ISYNTH, IMAGTYP, IMAGSRC, &
UNUSED9, UNUSED10, UNUSED11, UNUSED12, UNUSED13, UNUSED14, UNUSED15, UNUSED16, &
LEVEN, LPSPOL, LOVROK, LCALDA, UNUSED17
Character( Len = 16 ), Allocatable :: KEVNM
Character( Len = 8 ), Allocatable :: KSTNM, KHOLE, KO, KA, KT0, KT1, KT2, KT3, KT4, KT5, KT6, KT7, KT8, KT9, &
KF, KUSER0, KUSER1, KUSER2, KCMPNM, KNETWK, KWaveformRD, KINST
Integer(kind = 4),Dimension(70) :: Header1
Real(kind = 4),Dimension(40) :: Header2
Character(Len = 8), Dimension(22) :: Header3
real( kind = 4 ), Allocatable :: Waveform(:)
Integer(kind = 4) :: WaveformLong
End Type SACHeader
End module SACParameter
Program TestReadSAC
USE SACParameter
Implicit none
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
integer :: i, j, k, DELTAdt
real( kind = 4 ), Allocatable :: WaveformOne(:)
Open( Unit = FileUnit, File = 'FileName', Form = 'Unformatted', Status = 'OLD', Action = 'READ' )
Call GetFileN( FileUnit , FileName )
close( FileUnit )
Call RecordSACHeader( SACFILE1%Header1, SACFILE1%Header2, SACFILE1%Header3, SACFILE1%WaveformLong)
Do i = 1, EVENT1%StationNumber
Do j = 1, SACFILE1%WaveformLong
WaveformOne(j) = SACFILE1%Waveform(j)
End do
End do
Write(*,*)"please input the resample delta T!"
read(*,*)DELTAdt
Call Resample(WaveformOne, SACFILE1%WaveformLong, SACFILE1%DELTA, DELTAdt)
Call Toonewaveform( EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, EVENT1%StationWaveformLong, EVENT1%StationMaximumAmplitude)
Call WavePlots(EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, SACFILE1%DIST, SACFILE1%BAZ, SACFILE1%A, dtss, ntss, SACFILE1%KSTNM, nstr_sta, tlong, nnk)
End Program TestReadSAC
!---------------------------------------------
Subroutine GetFileN( iFileUnit )
USE SACParameter
USE, INTRINSIC :: ISO_FORTRAN_ENV
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Logical, Parameter :: FileIsNotEnd = .True.
Integer, Intent( INOUT ) :: iFileUnit
Integer :: i = 1, j = 1, k = 1
!Character( len = 1 ) ::cDummy
!GetFileN = 0
Rewind( iFileUnit )
Do i = 1,70
Read( Unit = iFileUnit, iostat = k, FMT = "(70G15.7)" ) SACFILE1%Header1(i)
!GetFileN = GetFileN + 1
End do
Do i = 1,40
Read( Unit = iFileUnit, iostat = k, FMT = "(40I10)" ) SACFILE1%Header2(i)
!GetFileN = GetFileN + 1
End do
Read( Unit = iFileUnit, iostat = k, FMT = "(A8)" ) SACFILE1%Header3(1)
!GetFileN = GetFileN + 1
Read( Unit = iFileUnit, iostat = k, FMT = "(A16)" ) SACFILE1%KEVNM
!GetFileN = GetFileN + 1
Do i = 2,22
Read( Unit = iFileUnit, iostat = k, FMT = "(21A8)" )SACFILE1%Header3(i)
End do
SACFILE1%WaveformLong = 0
Do while( FileIsNotEnd )
Read( Unit = iFileUnit, iostat = k) SACFILE1%Waveform( SACFILE1%WaveformLong )
!GetFileN = GetFileN + 1
SACFILE1%WaveformLong = SACFILE1%WaveformLong + 1
If( k == IOSTAT_END ) Exit
End do
Rewind( iFileUnit )
End Subroutine GetFileN
!-------------------------------------------
Subroutine RecordSACHeader(Header1, Header2, Header3, WaveformLong)
USE SACPArameter
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ), Dimension(70) :: Header1
Real(kind = 4), Intent( IN ), Dimension(40) :: Header2
Character(Len = 8), Intent( IN ), Dimension(22) :: Header3
Integer(kind = 4), Intent( IN ) :: WaveformLong
!此处为头段数据的记录，见下一篇帖子
End Subroutine RecordSACHeader
!----------------------------------
Subroutine toonewaveform(stationnumber, waveformlong, waveform, StationWaveformLong, StationMaximumAmplitude)
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Real(kind = 4), allocatable :: waveform(:,:)
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Integer(kind = 4), Intent( IN ) :: WaveformLong
Integer(kind = 4), Intent( IN ) :: stationnumber
allocate(waveform(waveformlong, stationnumber))
allocate(StationWaveformLong(stationnumber))
allocate(StationMaximumAmplitude(stationnumber))
StationWaveformLong(stationnumber) = WaveformLong
Do i = 1, stationnumber
StationMaximumAmplitude(i) = 1.E-15
Do j = 1, ntss(i)
If (StationMaximumAmplitude(i)<abs(waveform(j,i))) StationMaximumAmplitude(i) = abs(waveform(j,i))
End Do
Do j = 1, StationWaveformLong(i)
waveform(j, i) = waveform(j, i) / StationMaximumAmplitude(i)
End Do
End Do
Return
End Subroutine toonewaveform
!------------------------------------
Subroutine Resample(Waveform, WaveformLong, Delta, DeltaTdt)
!resampling for one dimension signal
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ) :: WaveformLong
real( kind = 4 ), Allocatable :: Waveform(:), WaveformTemp(:)
Integer(kind = 4) :: tLong
!Allocate(WaveformTemp(SACFILE1%WaveformLong))
!Allocate(Waveform(SACFILE1%WaveformLong))
tLong = ( SACFILE1%WaveformLong - 1 ) * SACFILE1%Delta
WaveformLongTemp = int( tLong / Deltadt + 0.5) + 1
Do i = 1, WaveformLongTemp
t = ( i - 1 ) * DELTAdt
Do j = 1, WaveformLong
t1 = ( j - 1 ) * Deltadt
t2 = j * Deltadt
If (t >= t1 .And. t < t2) Then
WaveformTemp( j ) = SACFILE1%Waveform( j ) + ( SACFILE1%Waveform( j + 1 ) - SACFILE1%Waveform( j )) / SACFILE1%Delta * ( t - t1 )
Goto 301
End If !j
End Do
301 Continue
End Do
Do i = 1, WaveformLongTemp
SACFILE1%Waveform(i) = WaveformTemp(i)
End Do
SACFILE1%WaveformLong = WaveformLongTemp
SACFILE1%Delta = DeltaTdt
Return
End Subroutine Resample
!----------------------------------------------
Subroutine BubbleSort(x, ny, n)
Dimension x(n), ny(n)
nsig = 1
m = n - 1
Do While ((m>0) .And. (nsig==1))
nsig = 2.
Do j = 1, m
If (x(j)>x(j+1)) Goto 100
nsig = 1
temp = x(j)
x(j) = x(j+1)
x(j+1) = temp
temp = ny(j)
ny(j) = ny(j+1)
ny(j+1) = temp
100 End Do
m = m - 1
End Do
Return
End Subroutine BubbleSort
F 币
贡献
权杖
窥视卡
发表于 2016-2-28 23:58:30
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
楼主
发表于 2016-2-29 16:39:04