Fortran批量读取SAC地震数据文件，写SAC文件，并转换成SEGY

zhangzhipeng · 发表于 2016-2-28 20:04:29

[Fortran] 纯文本查看 复制代码

Module SACParameter
    Integer, parameter :: FileUnit = 100
    Character( Len = * ), Parameter :: Directer = '"'//'..\WAV\WAVFormYiliang\yiliang\*.*'//'"'
    Type :: EventPara
        Integer( kind = 4 ) :: StationNumber
        Integer( kind = 4 ) :: FileNumber
        Character( Len = 4 ), allocatable :: StationName(:)        ! 台站名数组   
        Character( Len = 4 ), allocatable :: ComponentName(:,:)        ! 三分量文件名(台站,分量)
        Character( Len = 1 ), allocatable :: stringFileName(:,:)    ! 文件名信息数组二维字符数组
        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), Allocatable :: Header1(:,:)     ! 头端一(台站,分量)
        Real(kind = 4), Allocatable :: Header2(:,:)        ! 头端二(台站,分量)
        Character(Len = 8), Allocatable :: Header3(:,:)    ! 头端三(台站,分量)
        real( kind = 4 ), Allocatable :: Waveform(:,:,:)   ! 波形(台站,分量,记录)
        Integer(kind = 4), Allocatable :: WaveformLong(:,:)             ! 波形长度(台站,分量)
    End Type SACHeader 
End module SACParameter

Program Main
USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  External WriteName, SplitName

  !Character( Len = * ), Intent( In ) :: DirecterPath
  Integer :: i, j, k, m, n
  !Character( Len = * ), Parameter :: Path = '"'//'DirecterPath'//'"'
  Call DoWithWildcard( '..\WAV\WAVFormYiliang\yiliang\*.*' , WriteName , n )
  If ( N >= 0 ) then
    Write(*,*) '共' , N , '个文件！'
  End If
  pause 1
  Call DoWithWildcard( '..\WAV\WAVFormYiliang\yiliang\*.*' , SplitName , n )
  Do i = 1, EVENT1%StationNumber
    Write(*,*) EVENT1%StationName(i)
    Do j = 1, EVENT1%FileNumber
        Write(*,*) EVENT1%ComponentName(i,j)
    End do
  End do
End Program Main


Subroutine DoWithWildcard(cWildcard,CallBack,iTotal)
  Use DFLib,only:GetFileInfoQQ,GetLastErrorQQ,FILE$INFO,FILE$LAST,FILE$ERROR,FILE$FIRST,ERR$NOMEM,ERR$NOENT,FILE$DIR
  USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  Interface 
    Subroutine CallBack( FileName , loop )
      Character(*),Intent(In) :: FileName
      Integer,Intent(In) :: loop
    End Subroutine CallBack
  End Interface
  Character*(*),Intent(In)::cWildcard
  Integer,Intent(Out)::iTotal
  TYPE (FILE$INFO) info
  INTEGER(4)::Wildhandle,length,retInt
  Wildhandle = FILE$FIRST
  iTotal = 0
  DO WHILE (.TRUE.)
      length = GetFileInfoQQ(cWildCard,info,Wildhandle)
      IF ((Wildhandle .EQ. FILE$LAST) .OR.(Wildhandle .EQ. FILE$ERROR)) THEN
        SELECT CASE (GetLastErrorQQ())
        CASE (ERR$NOMEM)  !//内存不足
          iTotal = - 1
          Return
        CASE (ERR$NOENT)  !//碰到通配符序列尾
          Return
        CASE DEFAULT
          iTotal = 0
          Return
        END SELECT
      END IF
      If ((info%permit.AND.FILE$DIR).Eq.0) then
        Call CallBack( Trim(info.Name) , iTotal + 1 )
        iTotal = iTotal + 1
      End If
  END DO
  EVENT1%FileNumber = iTotal
  !Allocate(EVENT1%stringFileName(EVENT1%FileNumber))
End Subroutine DoWithWildcard

Subroutine WriteName( FileName , loop )
USE SACParameter
  Implicit None
    type(EventPara) :: EVENT1
    type(SACHeader) :: SACFILE1
  Character( Len = *),Intent(In) :: FileName
  Integer,Intent(In) :: loop
  Write(*,*) loop , FileName
End Subroutine WriteName

Subroutine SplitName( FileName , loop )          
USE SACParameter                                 
  Implicit None                                  
    type(EventPara) :: EVENT1                  
    type(SACHeader) :: SACFILE1                  
  Character(*),Intent(In) :: FileName             
  Integer,Intent(In) :: loop                    
  Integer( kind = 4) :: i, j, k, m = 0, n = 0, FileNameLong
  Character( Len = 1 ) :: separator = "."
  Character( Len = 4 ) :: Station
  FileNameLong = len_trim(FileName)
  j = Index( FileName, separator )
  k = Index( FileName, separator, .True. )
  If( Station /= FileName(1:j - 1)) then
      n = n + 1
      m = 1
      Station = FileName(1:j - 1)
      write(*,*) Station
      pause 2
      EVENT1%StationName(n) = FileName(1:j - 1)
      EVENT1%ComponentName(n,m) = FileName( k + 1:FileNameLong)
  Else
      EVENT1%ComponentName(n,m) = FileName( k + 1:FileNameLong)
      m = m + 1
  End if
  EVENT1%StationNumber = n
End Subroutine SplitName