diff --git a/.gitignore b/.gitignore index 259148f..3628328 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,14 @@ *.exe *.out *.app +*.x + +# Data files +*.win +*.win32 +*.cnt +*.sac +*.dat +test/testdata/* +*tmp* +*lst \ No newline at end of file diff --git a/README.md b/README.md index c542427..21206f3 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,45 @@ # fwin -A module for seismic waveform data in WIN and WIN32 formats +A module for reading seismic waveform data in WIN and WIN32 formats + +## Important Notice +This library is under being development. Significant changes in designs of module libraries might occur in future versions. + +## Compile + +```bash +$ cd src + (edit makefile if necessary) +$ make +``` + +To fully use the asynchronous I/O in Fortran2008, gfortran version 9 or later is required, while the program can be successfully compiled in more earlier versions. + +## Core modules + +Please read the block comments in the code for the detail of usage. Source codes of the utility programs (see below) may be useful for understanding the usage. + +### module m_win (in m_win.f90) +- `subroutine win__read_files`: Asynchronously read a set of WIN/WIN32 files +- `subroutine win__read_file`: Read a win file + +## module m_winch (in m_winch.f90) +- `type winch__hdr`: defines channel information. Usually use as an array +- `subroutine winch__read_tbl`: Read a channel table file +- `subroutine winch__get_all_chid`: Obtain all channel ID from the channel table data array +- `subroutine winch__get_all_stnm`: Obatin all station names contained in the channel table data array +- `subroutine winch__get_all_cmpnm`: Obtain all component names contained in the channel table data array + + +## Utility Programs + +### fchinf.x +Display selected channel information from a given channel table file + +### fdewin_s.x +Convert WIN/WIN32 files to ascii text: Synchronous version + +### fdewin_a.x +Convert WIN/WIN32 files to ascii text: Asynchronous version + +### fwin2sac.x +Convert WIN/WIN32 files to SAC-formatted datafiles diff --git a/src/fchinf.f90 b/src/fchinf.f90 new file mode 100644 index 0000000..3b505cb --- /dev/null +++ b/src/fchinf.f90 @@ -0,0 +1,88 @@ +!-------------------------------------------------------------------------------------------------! +!> Display channel information +!-- +program fchinf + + use iso_fortran_env + use m_win + use m_winch + use m_getopt + use m_util + implicit none + !-- + type(winch__hdr), allocatable :: ch(:) + character(256) :: fn_chtbl, chbuf, stbuf, cmpbuf + character(4), allocatable :: chid(:) + character(16), allocatable :: stnm(:) + character(16), allocatable :: cmpnm(:) + logical :: is_opt, is_opt_ch, is_opt_st, is_opt_cmp + integer :: nch, nst, ncmp + integer :: i, j, k + character(256) :: fmt + logical :: chid_mode + integer :: ichid + logical :: is_all_chid, is_all_st, is_all_cmp + character(4) :: chid_tmp + !---- + + call getopt('f', is_opt, fn_chtbl) + if(.not. is_opt) call usage_stop() + + call getopt('c', is_opt_ch, chbuf) + call getopt('s', is_opt_st, stbuf) + call getopt('p', is_opt_cmp, cmpbuf) + + ! priority is on specified channel ID than station&components + chid_mode = is_opt_ch + + if( (.not. is_opt_ch) .and. ( (.not. is_opt_st) .or. (.not. is_opt_cmp) ) ) call usage_stop() + + call winch__read_tbl(fn_chtbl, ch) + + if( chid_mode ) then + call util__read_arglst(chbuf, nch, is_all_chid, chid) + if( is_all_chid ) call winch__get_all_chid(ch, chid) + else + call util__read_arglst(stbuf, nst, is_all_st, stnm) + if( is_all_st ) call winch__get_all_stnm(ch, stnm) + call util__read_arglst(cmpbuf, ncmp, is_all_cmp, cmpnm) + if( is_all_cmp ) call winch__get_all_cmpnm(ch, cmpnm) + end if + + fmt = '(A, ": ", A, " (", A, ")", " unit: ", A, " T0=", F6.3, ' & + // '" h=", F6.3, " location: ", F9.5, " E,", F9.5, " N,", F9.2, " m")' + + if( chid_mode ) then + do j=1, size(ch) + do i=1, size(chid) + ichid = win__ach2ich(chid(i)) + if( ichid == ch(j)%ichid ) then + !! export channel information + write(output_unit, fmt) ch(i)%achid, trim(ch(i)%stnm), trim(ch(i)%cmpnm), & + trim(ch(i)%unit), ch(i)%period, ch(i)%damp, ch(i)%lon, ch(i)%lat, ch(i)%elev + exit + end if + end do + end do + else + do j=1, size(stnm) + do k=1, size(cmpnm) + call winch__st2chid(ch, stnm(j), cmpnm(k), chid_tmp, i) + if( i > 0 ) then + write(output_unit, fmt) ch(i)%achid, trim(ch(i)%stnm), trim(ch(i)%cmpnm), & + trim(ch(i)%unit), ch(i)%period, ch(i)%damp, ch(i)%lon, ch(i)%lat, ch(i)%elev + end if + end do + end do + end if + + contains + subroutine usage_stop() + + write(error_unit,'(A)') 'usage: fchinf.x <-f chtbl> [-c chid] [-s stnm] [-p cmpnm]' + stop + + end subroutine usage_stop + +end program fchinf +!-------------------------------------------------------------------------------------------------! diff --git a/src/fdewin_a.f90 b/src/fdewin_a.f90 new file mode 100644 index 0000000..b7eadb8 --- /dev/null +++ b/src/fdewin_a.f90 @@ -0,0 +1,111 @@ +!-------------------------------------------------------------------------------------------------! +!> Read Win/Win32-formatted seismograph data and export to ascii files: asynchronous version +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!-- +program fdewin_a + + use iso_fortran_env + use m_win + use m_util + use m_getopt + + implicit none + + integer, parameter :: fsmax = 200 + character(256), allocatable :: fn_win(:) + integer :: nch, nw, nsec + character(4), allocatable :: chid(:) + integer, allocatable :: dat(:,:) + integer, allocatable :: npts(:,:), sfreq(:) + character(80) :: d_out + logical :: is_test_mode + !! ---- + + !-----------------------------------------------------------------------------------------------! + !> command-line option processing + !-- + block + + character(80) :: fn_winlst + character(80) :: fn_chlst + logical :: is_opt, is_all + + call getopt('l', is_opt, fn_winlst, '' ) + + if( is_opt ) then + call util__readlst( fn_winlst, nw, fn_win ) + else + nw = 1 + allocate( fn_win(1) ) + call getopt('w', is_opt, fn_win(1), '') + if(.not. is_opt) call usage_stop() + end if + + call getopt('c', is_opt, fn_chlst, '' ) + + if( is_opt ) then + call util__read_arglst( fn_chlst, nch, is_all, chid ) + else + call usage_stop() + end if + + call getopt('d', is_opt, d_out, '.' ) + + call getopt('Z', is_test_mode) + + end block + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Read the data + !-- + block + integer :: tim + !---- + + allocate( sfreq(nch) ) + allocate( dat(fsmax*60*nw,nch) ) !! initial size + call win__read_files(fn_win, chid, sfreq, nsec, tim, dat, npts) + end block + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Export + !-- + block + integer :: i, j, io + character(80) :: fn_asc + !---- + + do i=1, nch + if( sfreq(i) > 0 ) then + fn_asc = trim(d_out) //'/'//trim(chid(i))//'.dat' + open(newunit=io, file=fn_asc, action='write', status='unknown') + do j=1, sfreq(i) * nsec + write(io,'(I0)') dat(j,i) + end do + close(io) + if( is_test_mode ) exit + else + write(error_unit,'(A)') 'Channel ' // chid(i) // ' : no data in the file' + end if + end do + end block + !-----------------------------------------------------------------------------------------------! + +contains + +!-----------------------------------------------------------------------------------------------! + subroutine usage_stop() + + write(error_unit,'(A)') 'usage: fdewin_a.x <-l winlst|-w winfile> <-c chid|lst> [-d dir]' + stop + end subroutine usage_stop + +end program fdewin_a +!-------------------------------------------------------------------------------------------------! diff --git a/src/fdewin_s.f90 b/src/fdewin_s.f90 new file mode 100644 index 0000000..9198bdd --- /dev/null +++ b/src/fdewin_s.f90 @@ -0,0 +1,122 @@ +!-------------------------------------------------------------------------------------------------! +!> Read Win/Win32-formatted seismograph data and export to ascii files: synchronous version +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!-- +program fdewin_s + + use iso_fortran_env + use m_win + use m_util + use m_getopt + + implicit none + !-- + + integer :: nch, nw, nsec + character(256), allocatable :: fn_win(:) + character(4), allocatable :: chid(:) + integer, allocatable :: dat(:,:), dat0(:,:) + integer, allocatable :: npts(:,:), sfreq(:) + character(80) :: d_out + integer, parameter :: fsmax = 200 + logical :: is_test_mode + !---- + + !-----------------------------------------------------------------------------------------------! + !> command-line option processing + !-- + block + + character(80) :: fn_winlst + character(80) :: fn_chlst + logical :: is_opt, is_all + + call getopt('l', is_opt, fn_winlst, '' ) + + if( is_opt ) then + call util__readlst( fn_winlst, nw, fn_win ) + else + nw = 1 + allocate( fn_win(1) ) + call getopt('w', is_opt, fn_win(1), '') + if( .not. is_opt ) call usage_stop() + end if + + call getopt('c', is_opt, fn_chlst, '' ) + + if( is_opt ) then + call util__read_arglst( fn_chlst, nch, is_all, chid ) + else + if( .not. is_opt ) call usage_stop() + end if + + call getopt('d', is_opt, d_out, '.' ) + + call getopt('Z', is_test_mode) + + end block + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Read the data + !-- + block + integer :: i, j + integer :: tim + !---- + + allocate( sfreq(nch) ) + allocate( dat(fsmax*60*nw,nch) ) !! initial size + dat(:,:) = 0 + + do i=1, nw + call win__read_file(fn_win(i), chid, sfreq, nsec, tim, dat0, npts) + do j=1, nch + if( sfreq(j) > 0 ) then + dat( (i-1)*sfreq(j)*nsec+1:i*sfreq(j)*nsec, j) = dat0(1:sfreq(j)*nsec,j) + end if + end do + + end do + end block + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Export + !-- + block + integer :: i, j, io + character(80) :: fn_asc + !---- + + do i=1, nch + if( sfreq(i) > 0 ) then + fn_asc = trim(d_out) //'/'//trim(chid(i))//'.dat' + open(newunit=io, file=fn_asc, action='write', status='unknown') + do j=1, sfreq(i)*nsec*nw + write(io,'(I0)') dat(j,i) + end do + close(io) + if( is_test_mode ) exit + else + write(error_unit,'(A)') 'Channel ' // chid(i) // ' : no data in the file' + end if + end do + end block + !-----------------------------------------------------------------------------------------------! + +contains + !-----------------------------------------------------------------------------------------------! + subroutine usage_stop() + + write(error_unit,'(A)') 'usage: fdewin_s.x <-l winlst|-w winfile> <-c chid|lst> [-d dir]' + stop + end subroutine usage_stop + +end program fdewin_s +!-------------------------------------------------------------------------------------------------! diff --git a/src/fwck.f90 b/src/fwck.f90 new file mode 100644 index 0000000..7f97b27 --- /dev/null +++ b/src/fwck.f90 @@ -0,0 +1,54 @@ +!-------------------------------------------------------------------------------------------------! +!> Display win/win32 file information +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!-- +program fwck + + use iso_fortran_env + use m_win + implicit none + !-- + character(256) :: fn + logical :: is_exist + character, allocatable :: buf(:) + character(4), allocatable :: chid(:) + integer :: nch + integer :: i + type(win__hdr) :: wh + integer :: io + !---- + + if( command_argument_count() /= 1 ) error stop 'usage: fwck.x winfile' + + call get_command_argument(1, fn) + inquire( file=fn, exist=is_exist ) + if( .not. is_exist ) error stop 'File not found: ' // trim(fn) + + call win__init() + call win__start_file_buf( fn, wh, buf, io ) + call win__finish_file_buf(io) + call win__scan_buf(wh, buf) + call win__get_chid( wh, buf, 1, nch, chid ) + + write(output_unit,'(A)') " FILE: " // trim( fn ) + write(output_unit,'(A,I8,A)') " SIZE: ", wh%sz, " bytes" + write(output_unit,'(A,I3)') " TYPE: ", wh%fmt + write(output_unit,*) + + do i=1, nch + write(output_unit,'(A,I5,A,A)') " Channel #", i, ": ", chid(i) + end do + + do i=1, wh%nb + call win__get_chid(wh, buf, i, nch, chid ) + write(*,'(A,I3,A,I4,2("-",I2.2)," ",2(I2.2,":"),I2.2)') & + "Block", i, ": ", wh%yr(i), wh%mo(i), wh%dy(i), wh%hr(i), wh%mi(i), wh%sc(i) + end do + +end program fwck +!-------------------------------------------------------------------------------------------------! diff --git a/src/fwin2sac.f90 b/src/fwin2sac.f90 new file mode 100644 index 0000000..620791b --- /dev/null +++ b/src/fwin2sac.f90 @@ -0,0 +1,245 @@ +!-------------------------------------------------------------------------------------------------! +!> win2sac fortran version +!-- +program fwin2sac + + use iso_fortran_env + use m_win + use m_winch + use m_getopt + use m_wsac + use m_util + implicit none + !-- + + character(4), allocatable :: chid(:) + integer :: nch + character(256), allocatable :: fn_win(:) + integer :: nwin + logical :: is_chtbl + character(256) :: dn + type(winch__hdr), allocatable :: ch(:), ch_tbl(:) + !---- + + !-----------------------------------------------------------------------------------------------! + ! win files + !-- + block + integer :: io, ierr, i + character(256) :: fn_winlst + logical :: is_opt + + call getopt('l', is_opt, fn_winlst) + if(.not. is_opt) call usage_stop() + open(newunit=io, file=fn_winlst, iostat=ierr, action='read', status='old') + if( ierr /= 0 ) error stop 'file not found: ' // trim(fn_winlst) + call util__countline(io, nwin) + allocate( fn_win(nwin) ) + do i=1, nwin + read(io, '(A)') fn_win(i) + end do + close(io) + + call getopt('d', is_opt, dn, '.') + end block + + !-----------------------------------------------------------------------------------------------! + ! channel table + !-- + block + character(256) :: fn_chtbl + call getopt('t', is_chtbl, fn_chtbl) + call winch__read_tbl(fn_chtbl, ch_tbl) + end block + + !-----------------------------------------------------------------------------------------------! + ! select channel IDs from command-line arguments + !-- + block + integer :: i, j, k, l + integer :: nst, ncmp + character(4) :: chid0 + character(16), allocatable :: stnm(:) + character(16), allocatable :: cmpnm(:) + character(256) :: chbuf, stbuf, cmpbuf + logical :: is_opt_ch, is_opt_st, is_opt_cmp + logical :: is_all_chid, is_all_st, is_all_cmp + integer :: ichid + !-- + + call getopt('c', is_opt_ch, chbuf) + call getopt('s', is_opt_st, stbuf) + call getopt('p', is_opt_cmp, cmpbuf) + + if( (.not. is_opt_ch) .and. ( (.not. is_opt_st) .or. (.not. is_opt_cmp) ) ) call usage_stop() + nch = 0 + ! priority is on specified channel ID than station & components + if( is_opt_ch ) then + call util__read_arglst(chbuf, nch, is_all_chid, chid) + if( is_all_chid ) then + call winch__get_all_chid(ch_tbl, chid) + nch = size(chid) + end if + + else + call util__read_arglst(stbuf, nst, is_all_st, stnm) + if( is_all_st ) call winch__get_all_stnm(ch_tbl, stnm) + call util__read_arglst(cmpbuf, ncmp, is_all_cmp, cmpnm) + if( is_all_cmp ) call winch__get_all_cmpnm(ch_tbl, cmpnm) + + nch = 0 + do i=1, size(stnm) + do j=1, size(cmpnm) + !try + call winch__st2chid(ch_tbl, stnm(i), cmpnm(j), chid0, k) + if( k > 0 ) nch = nch + 1 + end do + end do + allocate(chid(nch)) + l = 0 + do i=1, size(stnm) + do j=1, size(cmpnm) + !try + call winch__st2chid(ch_tbl, stnm(i), cmpnm(j), chid0, k) + if( k > 0 ) then + l = l + 1 + chid(l) = chid0 + end if + end do + end do + end if + !! prepare channel type data + if( is_chtbl ) then + allocate(ch(0)) + do i=1, nch + ichid = win__ach2ich(chid(i)) + do j=1, size(ch_tbl) + if( ichid == ch_tbl(j)%ichid ) then + ch_tbl(j)%conv = ch_tbl(j)%conv * 1d9 + ch = [ch, ch_tbl(j)] + exit + end if + end do + end do + else + allocate(ch(nch)) + do i=1, nch + call winch__init(ch(i)) + ch(i)%achid = chid(i) + ch(i)%stnm = ch(i)%achid + ch(i)%cmpnm = '' + ch(i)%conv = 1.0_real64 + end do + end if + end block + + !-----------------------------------------------------------------------------------------------! + ! Read win files and export to SAC file + !-- + block + integer, allocatable :: dat(:,:), npts(:,:), sfreq(:) + character(256) :: fn_sac + integer :: tim, nsec + integer :: i + type(sac__hdr) :: sh + character(8) :: ymd + character(6) :: hms + character(6) :: clen + + allocate(sfreq(nch)) + + call win__read_files(fn_win, ch(:)%achid, sfreq, nsec, tim, dat, npts) + call sac__init(sh) + call util__localtime(tim, & + sh%nzyear, sh%nzmonth, sh%nzday, sh%nzhour, sh%nzmin, sh%nzsec, sh%nzjday) + + sh%nzmsec = 0 + sh%b = 0 + + write(ymd,'(I4.4,I2.2,I2.2)') sh%nzyear, sh%nzmonth, sh%nzday + write(hms,'(3I2.2)') sh%nzhour, sh%nzmin, sh%nzsec + write(clen,'(I6.6)') nsec + + do i=1, nch + call ch2sh(ch(i), sh) + sh%npts = sfreq(i) * nsec + sh%delta = 1/dble(sfreq(i)) + sh%e = (sh%npts - 1) * sh%delta + + fn_sac = trim(dn) // '/' // ymd // '__' // hms // '__' // clen // '__' // & + trim(sh%kstnm) // '__' // trim(adjustl(sh%kcmpnm)) // '__.sac' + + if(sum(npts(:,i))>0) then + call sac__write(fn_sac, sh, dat(:,i)*ch(i)%conv, .true.) + end if + end do + end block + + contains + + subroutine usage_stop() + + character(18) :: sp1 = ' ' + write(error_unit,'(A)') 'usage: fwin2sac.x <-l listfile> ' + write(error_unit,'(A)') sp1//'[-t chtbl] [-c chids|chlist|all] [-s stnms|stlist|all] [-p cmpnm|cmplist|all]' + stop + end subroutine usage_stop + + subroutine ch2sh(ch, sh) + + type(winch__hdr), intent(in) :: ch + type(sac__hdr), intent(out) :: sh + !-- + integer :: iscan + !---- + + + sh%stla = ch%lat + sh%stlo = ch%lon + sh%stdp = - ch%elev + sh%stel = ch%elev + if( len_trim( ch%stnm ) <= 8 ) then + sh%kstnm = trim(ch%stnm) + else + sh%kstnm = ch%stnm(1:8) + end if + sh%kcmpnm = trim(ch%cmpnm) + + !! cmpaz, cmpinc + !! n + iscan = scan( ch%cmpnm, 'NX' ) + if( iscan > 0 ) then + sh%cmpaz = 0 + sh%cmpinc = 90 + end if + + !! e + iscan = scan( ch%cmpnm, 'EY' ) + if( iscan > 0 ) then + sh%cmpaz = 90 + sh%cmpinc = 90 + end if + + !! v + iscan = scan( ch%cmpnm, 'UZ' ) + if( iscan > 0 ) then + sh%cmpaz = 0 + sh%cmpinc = 0 + end if + + !! unit + select case ( trim(ch%unit) ) + case( 'm' ) + sh%idep = 6 + case( 'm/s' ) + sh%idep = 7 + case( 'm/s/s' ) + sh%idep = 8 + case default + sh%idep = 5 + end select + + end subroutine ch2sh + +end program fwin2sac +!-------------------------------------------------------------------------------------------------! diff --git a/src/m_getopt.f90 b/src/m_getopt.f90 new file mode 100644 index 0000000..bdb5c85 --- /dev/null +++ b/src/m_getopt.f90 @@ -0,0 +1,195 @@ +!-------------------------------------------------------------------------------------------------! +!> Processig command line option +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!< +!-- +module m_getopt + + use iso_fortran_env + use m_system + implicit none + private + + public :: getopt + + !-----------------------------------------------------------------------------------------------! + !> Obtain command-line option + !! + !! @par Usage + !! getopt( OPTNAME, exist, var, default_bar ) + !! - var, default_var are optional. + !! - types of var can be double, single integer or character + !-- + interface getopt + + module procedure getopt_d, getopt_s, getopt_i, getopt_c + + end interface getopt + !-----------------------------------------------------------------------------------------------! + +contains + + !-----------------------------------------------------------------------------------------------! + !> Obtain command line as characters + !! + !! @par Example + !! For the command + !! hoge -A -T 01 + !! call getopt( 'A', isA, Aval ) -> isA=.true., Aval='' + !! call getopt( 'T', isT, Tval ) -> isT=.true., Tval='01' + !! call getopt( 'C', isC, Cval ) -> isC=.false., Cval='' + !! call getopt( 'A', isA ) -> isA=.true. + !-- + subroutine getopt_c(opt, isExist, val, default) + character(*), intent(in) :: opt + logical, intent(out) :: isExist + character(*), optional, intent(out) :: val + character(*), optional, intent(in) :: default + !-- + integer :: narg + character(256), allocatable :: argv(:) + integer :: i + character(256) :: optkey + !---- + + narg = system__iargc() + allocate( argv(1:narg) ) + + do i=1, narg + call system__getarg( i, argv(i)(:) ) + end do + + if( present( val ) ) then + val = '' + end if + + isExist = .false. + optkey = '-'//trim(adjustl( opt ) ) + + do i=1, narg + + if( trim(optkey) == trim(argv(i)) ) then + + if( isExist ) then + write(error_unit,*) 'getopt: ', & + 'option '//trim(optkey)//' is multiplly defined. ' + end if + + isExist = .true. + + if( present( val ) ) then + + if( i==narg ) then + val = '' + else + val = argv(i+1) + end if + + end if + end if + end do + + deallocate( argv ) + + if( .not. isExist .and. present( default ) ) then + val = default + end if + + end subroutine getopt_c + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Obtain command line option for double precision number + !-- + subroutine getopt_d( opt, isExist, val, default ) + + character(*), intent(in) :: opt + logical, intent(out) :: isExist + real(real64), intent(out) :: val + real(real64), optional, intent(in) :: default + !-- + character(1024) :: aval + !---- + + call getopt_c( opt, isExist, aval ) + + if( .not. isExist ) then + if( present( default ) ) then + val = default + else + val = -99999.9_real64 + end if + return + end if + + read( aval, * ) val + + end subroutine getopt_d + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Obtain command line option for single precision number + !-- + subroutine getopt_s( opt, isExist, val, default ) + + character(*), intent(in) :: opt + logical, intent(out) :: isExist + real, intent(out) :: val + real, optional, intent(in) :: default + !-- + character(1024) :: aval + !---- + + call getopt_c( opt, isExist, aval ) + + if( .not. isExist ) then + if( present( default ) ) then + val = default + else + val = -99999.9 + end if + return + end if + + read(aval, *) val + + end subroutine getopt_s + !-----------------------------------------------------------------------------------------------! + + + !-----------------------------------------------------------------------------------------------! + !> Obtain command line option for integer number + !-- + subroutine getopt_i( opt, isExist, val, default ) + + character(*), intent(in) :: opt + logical, intent(out) :: isExist + integer, intent(out) :: val + integer, optional, intent(in) :: default + !-- + character(1024) :: aval + !---- + + call getopt_c( opt, isExist, aval ) + + if( .not. isExist ) then + if( present( default ) ) then + val = default + else + val = -99999 + end if + return + end if + + read(aval, *) val + + end subroutine getopt_i + !-----------------------------------------------------------------------------------------------! + +end module m_getopt +!-------------------------------------------------------------------------------------------------! diff --git a/src/m_system.f90 b/src/m_system.f90 new file mode 100644 index 0000000..ead6cc3 --- /dev/null +++ b/src/m_system.f90 @@ -0,0 +1,158 @@ +!-------------------------------------------------------------------------------------------------! +!> Processig command line argument, environment variables, and system call +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!< +!-- +module m_system + + use iso_fortran_env + implicit none + private + + public :: system__getarg + public :: system__getenv + public :: system__iargc + public :: system__expenv + + !!----------------------------------------------------------------------------------------------! + !> Obtain command line arguments for character, integer, single and double precision data + !-- + interface system__getarg + + module procedure getarg_a, getarg_i, getarg_f, getarg_d + + end interface + !-----------------------------------------------------------------------------------------------! + +contains + + !-----------------------------------------------------------------------------------------------! + !> Returns a number of arguments. Fortran2003 wrapper function + !-- + integer function system__iargc() + + system__iargc = command_argument_count() + + end function system__iargc + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Obtain environmental variable + !-- + subroutine system__getenv( name, value ) + + character(*), intent(in) :: name + character(*), intent(out) :: value + !---- + + call get_environment_variable( name, value ) + + end subroutine system__getenv + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> get i-th command line argument as a character + !-- + subroutine getarg_a (i, arg) + + integer, intent(in) :: i ! order of the arguments + character(*), intent(out) :: arg ! argument + !---- + + call get_command_argument( i, arg ) + + end subroutine getarg_a + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> get i-th command line argument as an integer + !-- + subroutine getarg_i (i, arg) + + integer, intent(in) :: i + integer, intent(out) :: arg + !-- + character(256) :: carg + !---- + + call getarg_a( i, carg ) + read(carg,*) arg + + end subroutine getarg_i + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> get i-th command line argument as a floating variable + !-- + subroutine getarg_f (i, arg) + + integer, intent(in) :: i + real, intent(out) :: arg + !-- + character(256) :: carg + !---- + + call getarg_a( i, carg ) + read(carg,*) arg + + end subroutine getarg_f + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> get i-th command line argument as a double precision variable + !-- + subroutine getarg_d (i, arg) + + integer, intent(in) :: i + real(real64), intent(out) :: arg + !-- + character(256) :: carg + !---- + + call getarg_a( i, carg ) + read(carg,*) arg + + end subroutine getarg_d + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Expand shell environmental variables. Variables must be wrapped in '${' and '}' + !-- + subroutine system__expenv( str ) + + character(*), intent(inout) :: str + !-- + character(256) :: str2 + integer :: ikey1, ikey2 + integer :: iptr + character(256) :: str3 + !---- + + iptr = 1 + str2 = '' + do + ikey1 = scan( str(iptr:), "${" ) + iptr - 1 + if( ikey1==iptr-1 ) exit + + ikey2 = scan( str(iptr:), "}" ) + iptr -1 + str2=trim(str2) // str(iptr:ikey1-1) + call system__getenv( str(ikey1+2:ikey2-1), str3 ) + str2 = trim(str2) // trim(str3) + iptr = ikey2+1 + + end do + str2 = trim(str2) // trim(str(iptr:)) + + str = trim(str2) + + end subroutine system__expenv + !-----------------------------------------------------------------------------------------------! + + +end module m_system +!-------------------------------------------------------------------------------------------------! diff --git a/src/m_util.f90 b/src/m_util.f90 new file mode 100644 index 0000000..a48eb1f --- /dev/null +++ b/src/m_util.f90 @@ -0,0 +1,230 @@ +!-------------------------------------------------------------------------------------------------! +!> Utility routines +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!< +!-- +module m_util + + use iso_fortran_env + implicit none + public + +contains + + !-----------------------------------------------------------------------------------------------! + !> Read command line buffer + !! given buffer (command-line option) may be one of the following: + !! - a filename contains data + !! - 'all' or 'ALL' + !! - comma-separated data + !-- + subroutine util__read_arglst( buf, n, is_all, lst ) + + character(*), intent(in) :: buf + integer, intent(out) :: n + logical, intent(out) :: is_all + character(*), allocatable, intent(out) :: lst(:) + + !-- + logical :: is_exist + integer :: io + integer :: i + !---- + + is_all = .false. + + !! first assume buf is a filename + inquire( file=buf, exist=is_exist ) + + if( is_exist ) then + + !! if it is filename, read list from it + open(newunit=io, file=buf, action='read', status='old') + + call util__countline(io, n) + + allocate( lst(n) ) + do i=1, n + read(io,'(A)') lst(i) + end do + close(io) + + is_all = .false. + + else if( trim(buf) == 'ALL' .or. trim(buf) == 'all' ) then + + is_all = .true. + n = 0 + + else + + !! buf should be a comma-separated channel number list + call util__split( buf, ',', n, lst ) + is_all = .false. + + end if + + end subroutine util__read_arglst + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + subroutine util__split(buf, sep, n, nm) + + character(*), intent(in) :: buf !< input buffer + character(*), intent(in) :: sep !< separator + integer, intent(out) :: n !< number of separated fields + !-- + character(*), allocatable, intent(out) :: nm(:) !< separated variables + integer :: i, j, k + !---- + + ! count separator + n = 0 + i = 1 + j = 1 + do + i = index(buf(j:), SEP) + if( i == 0 ) exit + j = j + i + n = n + 1 + end do + n = n + 1 !! #words = #separator + 1 + allocate(nm(n)) + + i = 1 + j = 1 + do k=1, n-1 + i = index(buf(j:), SEP) + nm(k) = trim(adjustl(buf(j:i+j-2))) + j = j + i + end do + nm(n) = trim(adjustl(buf(j:))) + + end subroutine util__split + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Count the line number except for the comment line + !-- + subroutine util__countline(io, n, comment) + + integer, intent(in) :: io + integer, intent(out) :: n + character, intent(in), optional :: comment + !-- + integer :: stat + character (256) :: line + !---- + + n = 0 + rewind(io) + do + read( io, '(A256)', iostat=stat) line + if( stat /= 0 ) exit + + if( present( comment ) ) then + if( line(1:1) == comment ) cycle + end if + + if( trim(line) /= '' ) then + n = n + 1 + end if + + end do + rewind( io ) + + end subroutine util__countline + !----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + subroutine util__readlst(fn, n, lst) + + character(*), intent(in) :: fn !< filename + integer, intent(out) :: n !< size + character(*), intent(out), allocatable :: lst(:) !< lst + !-- + character(256) :: abuf + integer :: io, ierr + integer :: i + !---- + + open(newunit=io, file=fn, action='read', status='old', iostat=ierr) + if( ierr /= 0 ) error stop "file "//trim(fn) //" not found" + call util__countline(io, n) + + allocate(lst(n)) + do i=1, n + read(io,'(A)') abuf + lst(i) = trim(adjustl(abuf)) + end do + + end subroutine util__readlst + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> Number of days since 0001-01-01 by Fairfield's formula + !-- + integer function fairfield(y, m, d) result(nd) + integer, intent(in) :: y, m, d + integer :: y0, m0 + + if( m==1 .or. m==2 ) then + y0 = y - 1 + m0 = m + 12 + else + y0 = y + m0 = m + end if + + nd = 365*y0 + floor(y0/4.) - floor(y/100.) + floor(y/400.) + floor(306*(m0+1)/10.) + d - 428 + + end function fairfield + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> calculate UNIX(POSIX) time from date & time + !-- + subroutine util__timelocal(yr, mo, dy, hr, mi, sc, tim) + + integer, intent(in) :: yr, mo, dy, hr, mi, sc + integer, intent(out) :: tim + integer, parameter :: tim0 = 719163 + + tim = (((fairfield(yr, mo, dy) - tim0) * 24 + hr )*60 + mi)*60 + sc + + end subroutine util__timelocal + !-----------------------------------------------------------------------------------------------! + + !-----------------------------------------------------------------------------------------------! + !> calculate date & time & number of days in the year from UNIX(POSIX) time + !-- + subroutine util__localtime(tim, yr, mo, dy, hr, mi, sc, jday) + + integer, intent(in) :: tim + integer, intent(out) :: yr, mo, dy, hr, mi, sc, jday + integer :: dy0, sc0 + dy0 = floor( tim / (24*60*60.) ) + 719163 + sc0 = tim - (dy0-719163) * 24*60*60 + + do yr=1970, 2099 + if( fairfield(yr+1,1,1) > dy0 ) exit + end do + do mo=1, 11 + if( fairfield(yr, mo+1,1) > dy0 ) exit + end do + dy = dy0 - fairfield(yr, mo, 1) + 1 + hr = sc0 / 3600 + mi = (sc0 - hr*3600) / 60 + sc = sc0 - hr*3600 - mi * 60 + jday = fairfield(yr, mo, dy) - fairfield(yr, 1, 1) + 1 + + end subroutine util__localtime + !-----------------------------------------------------------------------------------------------! + +end module m_util +!-------------------------------------------------------------------------------------------------! diff --git a/src/m_win.f90 b/src/m_win.f90 new file mode 100644 index 0000000..8949adb --- /dev/null +++ b/src/m_win.f90 @@ -0,0 +1,968 @@ +!-------------------------------------------------------------------------------------------------! +!> Read Win/Win32-formatted seismograph data +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!-- +module m_win + + use, intrinsic :: iso_fortran_env + use m_util + implicit none + private + + !! main subroutines for win/win32 files + public :: win__read_files + public :: win__read_file + + !! supporting public procedures & header type + public :: win__start_file_buf + public :: win__finish_file_buf + public :: win__decode_buf + public :: win__scan_buf + public :: win__get_chid + public :: win__ach2ich + public :: win__ich2ach + public :: win__init + public :: win__free + public :: win__hdr + + !! win header + type win__hdr + integer :: nb !< number of second blocks + integer, allocatable :: pb(:) !< pointer to the second-blocks + integer, allocatable :: yr(:), mo(:), dy(:) !< date of the second-blocks + integer, allocatable :: hr(:), mi(:), sc(:) !< time of the second-blocks + integer :: fmt !< WIN1(1) or WIN32(32) + integer :: sz !< file size in bytes + logical :: chkbuf = .false. + end type win__hdr + + !! private parameters + integer, private, parameter :: NS_MAX = 200 ! maximum samples per second + integer, private, parameter :: WIN1 = 1 + integer, private, parameter :: WIN32 = 32 + integer, private, parameter :: NB_BUF_INIT = 3600 !< initial size of the buffer memory (3600 s) + integer, private, parameter :: NO_CH = -99999 + + logical, save, private :: initialized = .false. + integer(int16), save, private :: twelve_bit + character(1), save, private :: azero !< binary image of zero in single-bit integer + + + integer, parameter :: YR_MIN = 1913 !! avoid misreading 201912 as 1912 + integer, parameter :: YR_MAX = 2100 + +contains + + !----------------------------------------------------------------------------------------------! + !> read single win data file for given channel ids chid. + !! if specified channel is not available, returns with npts(ich) = 0 + !-- + subroutine win__read_file(fn_win, chid, sfreq, nsec, tim, dat, npts) + + character(*), intent(in) :: fn_win + character(4), intent(in) :: chid(:) + integer, intent(out) :: sfreq(:) + integer, intent(out) :: nsec + integer, intent(out) :: tim + integer, intent(out), allocatable :: dat(:,:) + integer, intent(out), allocatable :: npts(:,:) + !-- + integer :: io_win + integer :: j + integer, allocatable :: dat0(:,:) + integer :: nch + character, allocatable :: buf(:) + type(win__hdr) :: wh + !---- + + nch = size(chid(:)) + + if( .not. initialized ) call win__init() + + write(error_unit,'(A)') '[win__read_file]: ' // trim(fn_win) + + call win__start_file_buf(fn_win, wh, buf, io_win ) + call win__finish_file_buf(io_win) + call win__scan_buf(wh, buf) + if( .not. allocated(npts) ) allocate(npts(wh%nb, nch)) + npts(:,:) = 0 + call win__decode_buf(wh, buf, nch, chid, dat0, npts, sfreq) + nsec = wh%nb + allocate(dat(maxval(sfreq)*nsec,nch)) + + deallocate(buf) + + call util__timelocal(wh%yr(1), wh%mo(1), wh%dy(1), wh%hr(1), wh%mi(1), wh%sc(1), tim) + + do j=1, nch + dat(1:sfreq(j)*nsec,j) = dat0(1:sfreq(j)*nsec,j) + end do + deallocate(dat0) + + call win__free(wh) + + end subroutine win__read_file + + !----------------------------------------------------------------------------------------------! + !> read multiple win data files having same temporal lengths + !-- + subroutine win__read_files(fn_win, chid, sfreq, nsec, tim, dat, npts) + + character(*), intent(in) :: fn_win(:) + character(4), intent(in) :: chid(:) + integer, intent(out) :: sfreq(:) + integer, intent(out) :: nsec !< number of seconds + integer, intent(out) :: tim !< Unix (POSIX) time of the head of the first file + integer, intent(out), allocatable :: dat(:,:) + integer, intent(out), allocatable :: npts(:,:) + !-- + integer :: nw + type(win__hdr) :: wh1, wh2 + integer :: io_win1, io_win2 + integer :: i, j + integer, allocatable :: dat0(:,:), npts0(:,:) + integer :: nb + integer :: nch + character, allocatable :: buf1(:), buf2(:) + logical :: first_touch + integer :: sfreq_max + !---- + + nw = size(fn_win(:)) + nch = size(chid(:)) + if( .not. initialized ) call win__init() + write(error_unit,'(A)') '[win__read_files]: ' // trim(fn_win(1)) + + first_touch = .true. + + call win__start_file_buf(fn_win(1), wh1, buf1, io_win1 ) + call win__finish_file_buf(io_win1) + + + do i=2, nw-1, 2 + write(error_unit,'(A)') '[win__read_files]: ' // trim(fn_win(i)) + call win__start_file_buf(fn_win(i), wh2, buf2, io_win2 ) + call win__scan_buf(wh1, buf1) + if(.not. allocated(npts0)) allocate(npts0(wh1%nb,nch)) + call win__decode_buf(wh1, buf1, nch, chid, dat0, npts0, sfreq) + + if( first_touch ) then + nb = wh1%nb + nsec = nb * nw + sfreq_max = maxval(sfreq) + allocate(dat(sfreq_max*nsec,nch)) + dat(:,:) = 0 + call util__timelocal(wh1%yr(1), wh1%mo(1), wh1%dy(1), wh1%hr(1), wh1%mi(1), wh1%sc(1), tim) + allocate(npts(nsec,nch)) + npts(:,:) = 0 + first_touch = .false. + else + if( maxval(sfreq) > sfreq_max ) then + call expand_dat(dat, nch, sfreq_max*nsec, maxval(sfreq)*nb*nw) + sfreq_max = maxval(sfreq) + end if + end if + + npts((i-2)*nb+1:(i-1)*nb,:) = npts0(1:nb,:) + + do j=1, nch + dat((i-2)*sfreq(j)*nb+1:(i-1)*sfreq(j)*nb, j) = dat0(1:sfreq(j)*nb,j) + end do + call win__free(wh1) + + call win__finish_file_buf(io_win2) + + write(error_unit,'(A)') '[win__read_files]: ' // trim(fn_win(i+1)) + call win__start_file_buf(fn_win(i+1), wh1, buf1, io_win1 ) + + call win__scan_buf(wh2, buf2) + call win__decode_buf(wh2, buf2, nch, chid, dat0, npts0, sfreq) + + deallocate(buf2) + + if( maxval(sfreq) > sfreq_max ) then + call expand_dat(dat, nch, sfreq_max*nsec, maxval(sfreq)*nb*nw) + sfreq_max = maxval(sfreq) + end if + + npts((i-1)*nb+1:i*nb,:) = npts0(1:nb,:) + do j=1, nch + dat((i-1)*sfreq(j)*nb+1:i*sfreq(j)*nb, j) = dat0(1:sfreq(j)*nb,j) + end do + call win__free(wh2) + + call win__finish_file_buf(io_win1) + + end do + + if( mod(nw, 2) == 0 ) then + write(error_unit,'(A)') '[win__read_files]: ' // trim(fn_win(nw)) + call win__start_file_buf(fn_win(nw), wh2, buf2, io_win2 ) + i = nw + else + i = nw + 1 + end if + + if( allocated(dat0) ) deallocate(dat0) + + call win__scan_buf(wh1, buf1) + if(.not. allocated(npts0)) allocate(npts0(wh1%nb,nch)) + call win__decode_buf(wh1, buf1, nch, chid, dat0, npts0, sfreq) + deallocate(buf1) + if( first_touch ) then + nb = wh1%nb + nsec = nb * nw + sfreq_max = maxval(sfreq) + allocate(dat(sfreq_max*nsec,nch)) + dat(:,:) = 0 + call util__timelocal(wh1%yr(1), wh1%mo(1), wh1%dy(1), wh1%hr(1), wh1%mi(1), wh1%sc(1), tim) + allocate(npts(nsec,nch)) + npts(:,:) = 0 + first_touch = .false. + else + if( maxval(sfreq) > sfreq_max ) then + call expand_dat(dat, nch, sfreq_max*nsec, maxval(sfreq)*nb*nw) + sfreq_max = maxval(sfreq) + end if + end if + npts((i-2)*nb+1:(i-1)*nb,:) = npts0(1:nb,:) + do j=1, nch + dat((i-2)*sfreq(j)*nb+1:(i-1)*sfreq(j)*nb, j) = dat0(1:sfreq(j)*nb,j) + end do + call win__free(wh1) + if( mod(nw, 2) == 0 ) then + call win__finish_file_buf(io_win2) + call win__scan_buf(wh2, buf2) + call win__decode_buf(wh2, buf2, nch, chid, dat0, npts0, sfreq) + deallocate(buf2) + + if( maxval(sfreq) > sfreq_max ) then + call expand_dat(dat, nch, sfreq_max*nsec, maxval(sfreq)*nb*nw) + sfreq_max = maxval(sfreq) + end if + + npts((i-1)*nb+1:i*nb,:) = npts0(1:nb,:) + do j=1, nch + dat((i-1)*sfreq(j)*nb+1:i*sfreq(j)*nb, j) = dat0(1:sfreq(j)*nb,j) + end do + call win__free(wh2) + end if + + deallocate(dat0, npts0) + + end subroutine win__read_files + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> initialize module + !-- + subroutine win__init() + + integer :: iptr + integer(int8) :: izero + !---- + + !! this variable will be used for reading 12-bit-width integer data from buffer + twelve_bit = 0 + + do iptr=0, 11 + twelve_bit = ibset(twelve_bit,iptr) + end do + + izero = 0 + azero = transfer( izero, azero ) + + initialized = .true. + + end subroutine win__init + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Open WIN1/WIN32 file and start reading all data into buffer buf(:) + !! The buffer memory will be allocated. + !! If buf(:) are already being used, these data will be destoroied + !-- + subroutine win__start_file_buf(fn, wh, buf, io) + + character(*), intent(in) :: fn !< input filename + type(win__hdr), intent(inout) :: wh !< header data + character, allocatable, intent(inout) :: buf(:) + integer, intent(out) :: io + !-- + integer :: ierr + !---- + + !! initialize module + if( .not. initialized ) call win__init() + + open(newunit=io, file=fn, access='stream', iostat=ierr, action='read', asynchronous='yes') + if(ierr == 0) then + inquire(io, size=wh%sz) + if(allocated(buf)) then + deallocate(buf) + end if + allocate(buf(wh%sz)) + read(io, asynchronous='yes') buf(:) + else + write(error_unit,*) 'file open error: ' // trim(fn) + stop + end if + + end subroutine win__start_file_buf + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Finalize reading process (win__start_file_buf) and close the file + !-- + subroutine win__finish_file_buf(io) + + integer, intent(in) :: io + logical :: is_open + !---- + + inquire(io, opened=is_open) + if( is_open ) close(io) + + end subroutine win__finish_file_buf + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! !> + !> Check the data buffer buf to + !! - investigate file type (WIN1 or WIN32) + !! - allocate necessary memory space based on the number of second-blocks + !! - create channel ID tables contained in the data buffer + !! This routine must be called before start reading data + !-- + subroutine win__scan_buf(wh, buf) + + type(win__hdr), intent(inout) :: wh !< win header + character, intent(in) :: buf(:) !< buffer + !-- + integer(int16) :: idb + integer :: sz_block + integer, allocatable :: pb(:) !! address of second blocks + integer :: nb_buf + integer :: n + integer :: block_header_size + integer :: dt1, dt2 !< location of date&time record + integer :: bs1, bs2 !< location of blocksize record + integer, allocatable :: yr(:), mo(:), dy(:), hr(:), mi(:), sc(:) + integer :: ichk + integer :: yr_32, mo_32, dy_32, hr_32, mi_32, sc_32 + !---- + + !! initialize module + if( .not. initialized ) call win__init() + + !! this routine is called only once per file + if( wh%chkbuf ) return + + !! initial size of buffer memory. Expanded if necessary + nb_buf = NB_BUF_INIT + allocate( pb(nb_buf), yr(nb_buf), mo(nb_buf), dy(nb_buf), hr(nb_buf), mi(nb_buf), sc(nb_buf) ) + + !! file type check: assume win32-type and evaluate if date & time are valid + ichk = int( transfer(buf(1:2), idb) ) + call decode_datetime(32, buf(5:12), yr_32, mo_32, dy_32, hr_32, mi_32, sc_32) + if( YR_MIN <= yr_32 .and. yr_32 <= YR_MAX .and. & + 1 <= mo_32 .and. mo_32 <= 12 .and. & + 1 <= dy_32 .and. dy_32 <= 31 .and. & + 0 <= hr_32 .and. hr_32 <= 23 .and. & + 0 <= mi_32 .and. mi_32 <= 59 .and. ichk == 0 ) then + wh%fmt = WIN32 + pb(1) = 5 + else ! if not, it must be WIN1 file + wh%fmt = WIN1 + pb(1) = 1 + end if + + if( wh%fmt == WIN1 ) then + block_header_size = 0 + dt1 = 4 + dt2 = 9 + bs1 = 0 + bs2 = 3 + else + block_header_size = 16 + dt1 = 0 + dt2 = 7 + bs1 = 12 + bs2 = 15 + end if + + n = 0 + !! obtain block locations & datetime information + do + + n = n + 1 + + call decode_datetime(wh%fmt, buf(pb(n)+dt1:pb(n)+dt2), & + yr(n), mo(n), dy(n), hr(n), mi(n), sc(n)) + + sz_block = transfer(buf(pb(n)+bs2:pb(n)+bs1:-1), sz_block ) + + !! location of the next block + pb(n+1) = pb(n) + block_header_size + sz_block + + !! EOF detection + if( pb(n+1) > wh% sz ) then + wh%nb = n + exit + end if + + !! expand buffer memory if necessary + if( n == nb_buf ) then + call expand_i( pb, nb_buf ) + end if + + end do + + allocate(wh%pb(wh%nb+1) ) + allocate(wh%yr(wh%nb), wh%mo(wh%nb), wh%dy(wh%nb), wh%hr(wh%nb), wh%mi(wh%nb), wh%sc(wh%nb)) + + wh%pb(1:wh%nb+1) = pb(1:wh%nb+1) + wh%yr(1:wh%nb) = yr(1:wh%nb) + wh%mo(1:wh%nb) = mo(1:wh%nb) + wh%dy(1:wh%nb) = dy(1:wh%nb) + wh%hr(1:wh%nb) = hr(1:wh%nb) + wh%mi(1:wh%nb) = mi(1:wh%nb) + wh%sc(1:wh%nb) = sc(1:wh%nb) + + deallocate(pb, yr, mo, dy, hr, mi, sc) + + wh%chkbuf = .true. + + end subroutine win__scan_buf + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Return a list of channel IDs in the ib-th block of the file buffer + !-- + subroutine win__get_chid( wh, buf, ib, nch, chid ) + + type(win__hdr), intent(inout) :: wh + character, intent(in) :: buf(:) + integer, intent(in) :: ib !< block number + integer, intent(out) :: nch + !-- + character(4), allocatable :: chid(:) + character(4), allocatable :: chid_tmp(:) + integer :: p + integer :: nmaxch + integer(int8) :: isb + integer(int16) :: idb + integer :: ss, ns + integer :: p1, p2 + !---- + + !! initialize module + if( .not. initialized ) call win__init() + + if( .not. wh%chkbuf ) call win__scan_buf(wh, buf) + + if( wh%fmt == WIN32 ) then + p1 = 16 !< block header size + p2 = 2 !< institution ID & network ID + else + p1 = 10 + p2 = 0 !< no institution&network IDs in WIN1 + end if + + if( ib > wh%nb ) then + write(error_unit,*) "ERROR [win__chid]: No such block ", ib + return + end if + + nmaxch = 300 + allocate(chid_tmp(nmaxch)) + if( allocated( chid ) ) deallocate(chid) + nch = 0 + p = wh%pb(ib) + p1 + + do while ( p < wh%pb(ib+1) ) + + !! temp memory expand + if( nch == nmaxch ) call expand_a( 4, chid_tmp, nmaxch ) + + ! read channel ID + p = p + p2 + write(chid_tmp(nch+1),'(Z4)') transfer( buf(p+1:p:-1), idb ) + isb = transfer( buf(p+2), isb ) + ss = int( ishft(isb, -4) ) + ns = int( iand( transfer( buf(p+3:p+2:-1), idb), twelve_bit )) + + nch = nch + 1 + if( ss == 0 ) then + p = p + 8 + ns/2 + else + p = p + 8 + (ns-1) * ss + end if + end do + + allocate( chid(nch) ) + chid(1:nch) = chid_tmp(1:nch) + deallocate(chid_tmp) + + !! zero padding + do p1=1, nch + do p2=1, 4 + if( chid(p1)(p2:p2) == ' ' ) chid(p1)(p2:p2) = '0' + end do + end do + + end subroutine win__get_chid + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Decode win data buffer for specified channel ids + !-- + subroutine win__decode_buf(wh, buf, nch, chid, dat, npts, sfreq) + + type(win__hdr), intent(inout) :: wh + character, intent(in) :: buf(:) !< win data buffer + integer, intent(in) :: nch !< #channels to decode + character(4), intent(in) :: chid(nch) !< channel IDs + integer, intent(out), allocatable :: dat(:,:) !< (npts,nch) + integer, intent(out) :: npts(wh%nb, nch) !< #samples. 0 for not-found + integer, intent(out) :: sfreq(nch) !< sampling frequency + !-- + integer :: ib + integer :: ns(nch) + integer :: dbuf(NS_MAX,nch) + integer :: ich + integer(int16) :: ichid(nch) + integer :: chid_tbl(-32768:32767) + integer :: i + integer :: sfreq_max + !---- + + !! initialize module + if( .not. initialized ) call win__init() + if( .not. wh%chkbuf ) call win__scan_buf(wh, buf) + do ich=1, nch + ichid(ich) = win__ach2ich( chid(ich) ) + end do + + ns(:) = -1 + npts(:,:) = 0 + + dbuf = 0 + sfreq(:) = 0 + !! set channel ID inverse table + chid_tbl(:) = NO_CH + do i=1, size(ichid) + chid_tbl(ichid(i)) = i + end do + + do ib=1, wh%nb + call decode_1sec( wh, buf, ib, ichid, chid_tbl, ns, dbuf(:,:)) + if( ib == 1 ) then + sfreq_max = max(maxval(ns(:)),1) + if(allocated(dat)) deallocate(dat) + allocate(dat(sfreq_max*wh%nb, nch)) + dat(:,:) = 0 + else + if( maxval(ns(:)) > sfreq_max ) then + !! expand memory size + call expand_dat(dat, nch, sfreq_max*wh%nb, maxval(ns(:))*wh%nb) + sfreq_max = maxval(ns(:)) + end if + end if + + do ich=1, nch + npts(ib,ich) = ns(ich) + if( ns(ich) > 0 ) then + dat((ib-1)*ns(ich)+1:ib*ns(ich), ich) = dbuf(1:ns(ich),ich) + sfreq(ich) = ns(ich) + end if + end do + end do + + end subroutine win__decode_buf + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Convert ascii channel ID to integer + !-- + function win__ach2ich( achid ) result(chid) + + character(4), intent(in) :: achid + !-- + integer(int16) :: chid + !---- + + read(achid,'(Z4)') chid + + end function win__ach2ich + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Convert integer channel ID to ascii + !-- + function win__ich2ach( chid ) result(achid) + + integer(int16), intent(in) :: chid + !-- + character(4) :: achid + !---- + + write(achid,'(Z4.4)') chid + + end function win__ich2ach + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Free memory + !-- + subroutine win__free( wh ) + + type(win__hdr), intent(inout) :: wh + !---- + + deallocate(wh%pb) + deallocate(wh%yr, wh%mo, wh%dy, wh%hr, wh%mi, wh%sc) + + wh%chkbuf = .false. + wh%nb = 0 + wh%sz = 0 + + end subroutine win__free + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Double the size of the input allocable array x(1:n) to x(1:2*n) + !! Array size n will be overwritten to 2*n + !-- + subroutine expand_i(x, n) + + integer, intent(inout), allocatable :: x(:) + integer, intent(inout) :: n + !-- + integer, allocatable :: tmp(:) + !---- + + if( size(x) /= n ) then + write(error_unit,*) "Error [expand]: array size mismatch" + return + end if + + allocate( tmp(n) ) + tmp(:) = x(:) + deallocate(x) + allocate(x(2*n)) + x(1:n) = tmp(:) + deallocate(tmp) + n = 2 * n + + end subroutine expand_i + !----------------------------------------------------------------------------------------------! + !----------------------------------------------------------------------------------------------! + !> Double the size of the input allocable array x(1:n) to x(1:2*n) + !! Array size n will be overwritten to 2*n + !-- + subroutine expand_a(m, x, n) + + integer, intent(in) :: m + character(m), intent(inout), allocatable :: x(:) + integer, intent(inout) :: n + !-- + character(m), allocatable :: tmp(:) + !---- + + if( size(x) /= n ) then + write(error_unit,*) "Error [expand]: array size mismatch" + return + end if + + allocate( tmp(n) ) + tmp(:) = x(:) + deallocate(x) + allocate(x(2*n)) + x(1:n) = tmp(:) + deallocate(tmp) + n = 2 * n + + end subroutine expand_a + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> memory size change dat(npts0,nch) -> dat(npts1,nch) + !-- + subroutine expand_dat(dat, nch, npts0, npts1) + + integer, intent(inout), allocatable :: dat(:,:) + integer, intent(in) :: nch, npts0, npts1 + !-- + integer, allocatable :: dtmp(:,:) + !---- + + if( size(dat, dim=1) /= npts0 ) then + write(error_unit,*) '[expand_dat]: size mismatch' + return + end if + allocate(dtmp(npts0,nch)) + dtmp(:,:) = dat(:,:) + deallocate(dat) + allocate(dat(npts1,nch)) + dat(1:npts0,:) = dtmp(:,:) + dat(npts0+1:npts1,:) = 0 + deallocate(dtmp) + + end subroutine expand_dat + !----------------------------------------------------------------------------------------------! + !----------------------------------------------------------------------------------------------! + !> memory size change dat(npts0,nch) -> dat(npts1,nch) + !-- + subroutine expand_rdat(dat, nch, npts0, npts1) + + real, intent(inout), allocatable :: dat(:,:) + integer, intent(in) :: nch, npts0, npts1 + !-- + real, allocatable :: dtmp(:,:) + !---- + + if( size(dat, dim=1) /= npts0 ) then + write(error_unit,*) '[expand_dat]: size mismatch' + return + end if + allocate(dtmp(npts0,nch)) + dtmp(:,:) = dat(:,:) + deallocate(dat) + allocate(dat(npts1,nch)) + dat(1:npts0,:) = dtmp(:,:) + dat(npts0+1:npts1,:) = 0 + deallocate(dtmp) + + end subroutine expand_rdat + !----------------------------------------------------------------------------------------------! + !----------------------------------------------------------------------------------------------! + !> Read win/win32-formatted date and time in the second block header + !-- + subroutine decode_datetime( fmt, tbuf, yr, mo, dy, hr, mi, sc ) + + integer, intent(in) :: fmt + character, intent(in) :: tbuf(:) + integer, intent(out) :: yr, mo, dy, hr, mi, sc + !-- + integer(int8) :: hb1, hb2 + integer :: y1, y2 + !---- + + if( fmt == WIN1 ) then + call bcd_int(tbuf(1),hb1,hb2); y1 = hb1*10+hb2 + call bcd_int(tbuf(2),hb1,hb2); mo = hb1*10+hb2 + call bcd_int(tbuf(3),hb1,hb2); dy = hb1*10+hb2 + call bcd_int(tbuf(4),hb1,hb2); hr = hb1*10+hb2 + call bcd_int(tbuf(5),hb1,hb2); mi = hb1*10+hb2 + call bcd_int(tbuf(6),hb1,hb2); sc = hb1*10+hb2 + + !! assume observation period 1981 -- 2080 + if( y1 > 80 ) then + yr = y1 + 1900 + else + yr = y1 + 2000 + end if + + else if ( fmt == WIN32 ) then + call bcd_int(tbuf(1),hb1,hb2); y1 = hb1*10+hb2 + call bcd_int(tbuf(2),hb1,hb2); y2 = hb1*10+hb2 + call bcd_int(tbuf(3),hb1,hb2); mo = hb1*10+hb2 + call bcd_int(tbuf(4),hb1,hb2); dy = hb1*10+hb2 + call bcd_int(tbuf(5),hb1,hb2); hr = hb1*10+hb2 + call bcd_int(tbuf(6),hb1,hb2); mi = hb1*10+hb2 + call bcd_int(tbuf(7),hb1,hb2); sc = hb1*10+hb2 + yr = y1 * 100 + y2 + end if + + end subroutine decode_datetime + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Decode bcd-formatted buffer into two integers + !-- + subroutine bcd_int( b, j, k ) + + character, intent(in) :: b + integer(int8), intent(out) :: j + integer(int8), intent(out) :: k + !-- + integer(int8) :: i + integer :: l + !---- + + i = transfer( b, i ) + j=ishft(i,-4) + k=i + do l=4,7 + k=ibclr(k,l) + end do + + end subroutine bcd_int + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + subroutine hbyte_int( i, j, k ) + + integer(int8), intent(in) :: i + integer(int8), intent(out) :: j ! i(1:4) + integer(int8), intent(out) :: k ! i(5:8) + !-- + integer :: l + !---- + + !! Decompose + j=ishft(i,-4) + k=i + do l=4,7 + k=ibclr(k,l) + end do + + !! Positive/Negative + if( btest(j,3) ) then + j=ibset(j,4) + j=ibset(j,5) + j=ibset(j,6) + j=ibset(j,7) + end if + + if( btest(k,3) ) then + k=ibset(k,4) + k=ibset(k,5) + k=ibset(k,6) + k=ibset(k,7) + end if + + end subroutine hbyte_int + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Read specified second-block and decode it for given channels + !-- + subroutine decode_1sec( wh, buf, ib, chid, chid_tbl, ns, dat ) + + type(win__hdr), intent(in) :: wh + character, intent(in) :: buf(:) + integer, intent(in) :: ib !< block number + integer(int16), intent(in) :: chid(:) + integer, intent(in) :: chid_tbl(-32768:32767) + integer, intent(out) :: ns(:) + integer, intent(out) :: dat(:,:) + !-- + integer(int8) :: isb + integer(int8) :: isb1, isb2 + integer(int16) :: idb + integer(int16) :: chid_buf + integer :: p, q + integer :: ss, ns0 + integer :: p1, p2 + integer :: ich + integer :: i + integer :: nch + !---- + + nch = size(chid) + + if( wh%fmt == WIN32 ) then + p1 = 16 !< block header size + p2 = 2 !< institution ID & network ID + else + p1 = 10 + p2 = 0 !< no institution&network IDs in WIN1 + end if + + if( ib > wh%nb ) then + write(error_unit,*) "ERROR [decode_1sec]: No such block ", ib + return + end if + + p = wh%pb(ib) + p1 + + !! initialize + ns(:) = -1 + + do while ( p < wh%pb(ib+1) ) + + !! skip institution & network ID + p = p + p2 + + !! read channel ID + chid_buf = transfer( buf(p+1:p:-1), idb ) + + !! decompose sample size ss (0.5 byte) & sample length ns (1.5 byte) + isb = transfer( buf(p+2), isb ) + ss = int( ishft(isb, -4) ) + ns0 = int( iand( transfer( buf(p+3:p+2:-1), idb), twelve_bit )) + + ich = chid_tbl(chid_buf) + + if( ich /= NO_CH ) then + + ns(ich) = ns0 + + !! first data + dat(1,ich) = transfer( buf(p+7:p+4:-1), ich ) + + select case( ss ) + case( 0 ) !! 0.5 byte + + do i=2, ns0-2, 2 + q = 8 + p + (i-2) / 2 + call hbyte_int( transfer(buf(q:q), isb), isb1, isb2) + dat(i, ich) = dat(i-1,ich) + int(isb1) + dat(i+1,ich) = dat(i-1,ich) + int(isb1) + int(isb2) + end do + if( ns0 /= 1 ) then + q = 8 + p + (ns0-2) / 2 + call hbyte_int( transfer(buf(q:q), isb), isb1, isb2) + dat(ns0, ich) = dat(ns0-1,ich) + int(isb1) + end if + case( 1 ) + do i=2, ns0 + q = 8 + p + (i-2) + dat(i,ich) = dat(i-1,ich) + int( transfer( buf(q:q), isb) ) + end do + case( 2 ) + do i=2, ns0 + q = 8 + p + (i-2) * 2 + dat(i,ich) = dat(i-1,ich) + int( transfer( buf(q+1:q:-1), idb) ) + end do + case( 3 ) + do i=2, ns0 + q = 8 + p + (i-2) * 3 + dat(i,ich) = dat(i-1,ich) + int( transfer( azero // buf(q+2)//buf(q+1)//buf(q), i ) ) / 256 + end do + case( 4 ) + do i=2, ns0 + q = 8 + p + (i-2) * 4 + dat(i,ich) = dat(i-1,ich) + transfer(buf(q+3:q:-1), i) + end do + case( 5 ) + do i=2, ns0 + q = 8 + p + (i-2) * 4 + dat(i,ich) = transfer(buf(q+3:q:-1), i) + end do + end select + end if + + if( ns0 /= 1 ) then + if ( ss == 0 ) then + p = p + 8 + ns0/2 + else + p = p + 8 + (ns0-1) * ss + end if + else + p = p + 8 + end if + + end do + + end subroutine decode_1sec + !----------------------------------------------------------------------------------------------! + +end module m_win +!------------------------------------------------------------------------------------------------! diff --git a/src/m_winch.f90 b/src/m_winch.f90 new file mode 100644 index 0000000..6869157 --- /dev/null +++ b/src/m_winch.f90 @@ -0,0 +1,331 @@ +!-------------------------------------------------------------------------------------------------! +!> WIN/WIN32 Channel Table +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. +!-- +module m_winch + + use, intrinsic :: iso_fortran_env + use m_win + implicit none + + public :: winch__hdr + + type winch__hdr + + character(4) :: achid !< channel id in 4-digit character + integer(int16) :: ichid !< channel id in integer form + character(16) :: stnm !< station name + character(16) :: cmpnm !< component name + character(16) :: unit !< measurement unit + real(real64) :: period !< natural period of the sensor + real(real64) :: damp !< damping constant + real(real64) :: lon !< station longitude (degrees-east) + real(real64) :: lat !< station latitude (degrees-north) + real(real64) :: elev !< station elevation (m) + real(real64) :: sens !< sensor sensitivity + integer :: ampl !< amplification factor + real(real64) :: step !< 1-digit step in voltage + real(real64) :: conv !< conversion coef. from DC to physical unit + + end type winch__hdr + +contains + + !----------------------------------------------------------------------------------------------! + !> Initialize win channel header + !-- + subroutine winch__init(wch) + + type(winch__hdr), intent(inout) :: wch + !---- + + wch%achid = '-123' + wch%ichid = 0 + wch%stnm = '-12345.0' + wch%cmpnm = '-12345.0' + wch%unit = '-12345.0' + wch%ampl = -12345 + wch%period = -12345.0_real64 + wch%damp = -12345.0_real64 + wch%lon = -12345.0_real64 + wch%lat = -12345.0_real64 + wch%elev = -12345.0_real64 + wch%sens = -12345.0_real64 + wch%step = -12345.0_real64 + wch%conv = -12345.0_real64 + + end subroutine winch__init + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Read the spacified channel table file + !-- + subroutine winch__read_tbl(fn_chtbl, wch) + + character(*), intent(in) :: fn_chtbl + type(winch__hdr), intent(out), allocatable :: wch(:) + !! -- + integer :: io + integer(int16) :: ich + integer :: ierr + character(512) :: line + type(winch__hdr) :: wch0 + !! ---- + + allocate(wch(0)) + + open(newunit=io, file=trim(fn_chtbl), iostat=ierr, status='old', action='read') + if( ierr /= 0 ) then + write(error_unit,'(A)') 'error [winch__readchtbl]: file not found.' + return + end if + + do + read(io, '(A)', iostat=ierr) line + if( ierr/= 0 ) exit + + line = adjustl(line) + + ! skip comment or blank lines + if( line(1:1) .eq. "#" ) cycle + if( len_trim(line) == 0 ) cycle + + call winch__init(wch0) + + ich = win__ach2ich(line(1:4)) + + ! store the channel information + wch0 % ichid = ich + wch0 % achid = line(1:4) + + ! the space (separation) should be searched via index function rather than automated read + ! command, as the channel table info may contain other separator characters such as '/' + + call readline_skip() + call readline_skip() + call readline_skip() + + call readline_ascii (wch0%stnm) ! station name + call readline_ascii (wch0%cmpnm) ! station name + + call readline_skip() + call readline_skip() + + call readline_real64 (wch0%sens, 1.0_real64) ! sensor sensitivity + call readline_ascii (wch0%unit) ! sensor unit + call readline_real64 (wch0%period, -12345.0_real64) ! natural period + call readline_real64 (wch0%damp, -12345.0_real64) ! damping constant + call readline_integer(wch0%ampl, -12345) ! amplification factor + call readline_real64 (wch0%step, 1.0_real64) ! 1 digit step in voltage + call readline_real64 (wch0%lat, -12345.0_real64) ! station latitude + call readline_real64 (wch0%lon, -12345.0_real64) ! station longitude + call readline_real64 (wch0%elev, -12345.0_real64) ! station elevation + + ! conversion coefficient + wch0%conv = wch0%step / (wch0%sens * 10**(dble(wch0%ampl)/20.0_real64)) + + !! expand array + wch = [wch, wch0] + + end do + + contains + + subroutine readline_skip() + integer :: idx + !---- + idx=index(line,' ') + line=adjustl(line(idx:len_trim(line))) + + end subroutine readline_skip + + subroutine readline_ascii(var) + character(*), intent(out) :: var + integer :: idx + !---- + idx=index(line,' ') + read(line(1:idx),'(A)') var + line=adjustl(line(idx:len_trim(line))) + + end subroutine readline_ascii + + subroutine readline_real64(var, default_var) + real(real64), intent(out) :: var + real(real64), intent(in) :: default_var + integer :: idx + + idx=index(line,' ') + if( trim(line(1:idx)) == "*" ) then + var= default_var + else + read(line(1:idx),*) var + end if + line=adjustl(line(idx:len_trim(line))) + + end subroutine readline_real64 + + subroutine readline_integer(var, default_var) + integer, intent(out) :: var + integer, intent(in) :: default_var + integer :: idx + + idx=index(line,' ') + if( trim(line(1:idx)) == "*" ) then + var = default_var + else + read(line(1:idx),*) var + end if + line=adjustl(line(idx:len_trim(line))) + end subroutine readline_integer + + end subroutine winch__read_tbl + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Return all channel IDs in the list + !-- + subroutine winch__get_all_chid(wch, chid) + + type(winch__hdr), intent(in) :: wch(:) + character(4), intent(out), allocatable :: chid(:) + + allocate(chid(size(wch))) + chid(:) = wch(:)%achid + + end subroutine winch__get_all_chid + + !----------------------------------------------------------------------------------------------! + !> Return all station names in the channel list without duplications + !-- + subroutine winch__get_all_stnm(wch, stnm) + + type(winch__hdr), intent(in) :: wch(:) + character(*), intent(out), allocatable :: stnm(:) + !-- + integer :: i, j + character(len(stnm)) :: stnm0 + logical :: found + integer :: nch, nst + !---- + + nch = size(wch) + + if( nch <= 0 ) then + nst = 0 + return + end if + + nst = 0 + do i=1, nch + if( .not. allocated(stnm) ) then + nst = 1 + allocate(stnm(nst)) + stnm(nst) = trim(adjustl(wch(i)%stnm)) + end if + + stnm0 = trim(adjustl(wch(i)%stnm)) + found = .false. + + do j=1, nst + found = trim(adjustl(stnm0)) == trim(adjustl(stnm(j))) + if( found ) exit + end do + + if( .not. found ) then + nst = nst + 1 + stnm = [stnm, stnm0] !! expand the array + end if + + end do + + end subroutine winch__get_all_stnm + !----------------------------------------------------------------------------------------------! + + !----------------------------------------------------------------------------------------------! + !> Return all station names in the channel list without duplications + !-- + subroutine winch__get_all_cmpnm(wch, cmpnm) + + type(winch__hdr), intent(in) :: wch(:) + character(*), intent(out), allocatable :: cmpnm(:) + !-- + integer :: i, j + character(len(cmpnm)) :: cmpnm0 + logical :: found + integer :: nch, ncmp + !---- + + nch = size(wch) + + if( nch <= 0 ) then + ncmp = 0 + return + end if + + ncmp = 0 + do i=1, nch + if( .not. allocated(cmpnm) ) then + ncmp = 1 + allocate(cmpnm(ncmp)) + cmpnm(ncmp) = trim(adjustl(wch(i)%cmpnm)) + end if + + cmpnm0 = trim(adjustl(wch(i)%cmpnm)) + found = .false. + + do j=1, ncmp + found = trim(adjustl(cmpnm0)) == trim(adjustl(cmpnm(j))) + if( found ) exit + end do + + if( .not. found ) then + ncmp = ncmp + 1 + cmpnm = [cmpnm, cmpnm0] !! expand the array + end if + + end do + + end subroutine winch__get_all_cmpnm + !----------------------------------------------------------------------------------------------! + + + !----------------------------------------------------------------------------------------------! + !> Return channel ID corresponding to given station and component name + !-- + subroutine winch__st2chid(wch, stnm, cmpnm, chid, ikey) + + type(winch__hdr), intent(in) :: wch(:) + character(*), intent(in) :: stnm !< station name + character(*), intent(in) :: cmpnm !< component name + character(4), intent(out) :: chid !< channel id + integer, intent(out) :: ikey !< location in the wch(:), -1 for failure + !-- + integer :: i + !---- + + do i=1, size(wch) + if( trim(adjustl( wch(i)%stnm )) == trim(adjustl(stnm)) .and. & + trim(adjustl( wch(i)%cmpnm )) == trim(adjustl(cmpnm)) ) then + chid = wch(i)%achid + ikey = i + return + end if + + end do + + ! not found + ikey = -1 + chid = 'XXXX' + return + + end subroutine winch__st2chid + !----------------------------------------------------------------------------------------------! + + +end module m_winch +!-------------------------------------------------------------------------------------------------! \ No newline at end of file diff --git a/src/m_wsac.f90 b/src/m_wsac.f90 new file mode 100644 index 0000000..42067f2 --- /dev/null +++ b/src/m_wsac.f90 @@ -0,0 +1,485 @@ +!------------------------------------------------------------------------------------------------! +!> Write SAC-formatted seismograms +!! +!! @copyright +!! Copyright (c) 2019 Takuto Maeda. All rights reserved. +!! +!! @license +!! This software is released under the MIT license. See LICENSE for details. !-- +module m_wsac + + use iso_fortran_env, only: DP=>real64, SP=>real32, error_unit, input_unit + + implicit none + private + + !! -- Public Procedures + public :: sac__hdr ! sac data type + public :: sac__write ! write sac datafile + public :: sac__init ! initialize sac data type + public :: sac__whdr ! read header + + !----------------------------------------------------------------------------------------------! + !> Sac Header Type Definition + !-- + type sac__hdr + + !! var name description record# + real(DP) :: delta ! sampling interval (001) + real(DP) :: depmin ! minimum value (002) + real(DP) :: depmax ! maximum value (003) + real(DP) :: b ! begenning independent value (006) + real(DP) :: e ! ending independent value (007) + real(DP) :: o ! event origin time (008) + real(DP) :: a ! first arrival time (009) + real(DP) :: t0 ! time picks (011) + real(DP) :: t1 ! time picks (012) + real(DP) :: t2 ! time picks (013) + real(DP) :: t3 ! time picks (014) + real(DP) :: t4 ! time picks (015) + real(DP) :: t5 ! time picks (016) + real(DP) :: t6 ! time picks (017) + real(DP) :: t7 ! time picks (018) + real(DP) :: t8 ! time picks (019) + real(DP) :: t9 ! time picks (020) + real(DP) :: stla ! station latitude (032) + real(DP) :: stlo ! station longitude (033) + real(DP) :: stel ! station elevation (m) (034) + real(DP) :: stdp ! station depth (m) (035) + real(DP) :: evla ! event latitude (036) + real(DP) :: evlo ! event longitude (037) + real(DP) :: evel ! event elevation (m) (038) + real(DP) :: evdp ! event depth (m) (039) + real(DP) :: mag ! event magnitude (040) + real(DP) :: user0 ! user header (041) + real(DP) :: user1 ! user header (042) + real(DP) :: user2 ! user header (043) + real(DP) :: user3 ! user header (044) + real(DP) :: user4 ! user header (045) + real(DP) :: user5 ! user header (046) + real(DP) :: user6 ! user header (047) + real(DP) :: user7 ! user header (048) + real(DP) :: user8 ! user header (049) + real(DP) :: user9 ! user header (050) + real(DP) :: dist ! distance (km) (051) + real(DP) :: az ! azimuth (deg) (052) + real(DP) :: baz ! back azimuth (deg) (053) + real(DP) :: gcarc ! angular distance (deg) (054) + real(DP) :: depmen ! mean value (057) + real(DP) :: cmpaz ! component azimuth (058) + real(DP) :: cmpinc ! component incident angle (059) + integer :: nzyear ! reference time, year (071) + integer :: nzjday ! reference time, julian day (072) + integer :: nzhour ! reference time, hour (073) + integer :: nzmin ! reference time, minute (074) + integer :: nzsec ! reference time, second (075) + integer :: nzmsec ! reference time, millisecond (076) + integer :: nvhdr ! header version (077) + integer :: npts ! number of data points (080) + integer :: iftype ! type of file (086) + integer :: idep ! type of dependent var. (087) + integer :: ievtyp ! event type (093) + logical :: leven ! is evenly spaced file (106) + logical :: lpspol ! is positive polarity (107) + logical :: lovrok ! is overwrite ok? (108) + logical :: lcalda ! is calc distance azimuth (109) + character(8) :: kstnm ! station name (111) + character(16) :: kevnm ! event name (113) + character(8) :: khole ! hole name (117) + character(8) :: ko ! origin time identification (119) + character(8) :: ka ! time pick name (121) + character(8) :: kt0 ! time pick name (123) + character(8) :: kt1 ! time pick name (125) + character(8) :: kt2 ! time pick name (127) + character(8) :: kt3 ! time pick name (129) + character(8) :: kt4 ! time pick name (131) + character(8) :: kt5 ! time pick name (133) + character(8) :: kt6 ! time pick name (135) + character(8) :: kt7 ! time pick name (137) + character(8) :: kt8 ! time pick name (139) + character(8) :: kt9 ! time pick name (141) + character(8) :: kf ! fini identification (143) + character(8) :: kuser0 ! user area (145) + character(8) :: kuser1 ! user area (147) + character(8) :: kuser2 ! user area (149) + character(8) :: kcmpnm ! component name (151) + character(8) :: knetwk ! network name (153) + character(8) :: kdatrd ! date data onto comp. (155) + character(8) :: kinst ! instrument (157) + + !! Unofficial header at unused blocks + real(DP) :: user10 ! user header (064) + real(DP) :: user11 ! user header (065) + real(DP) :: user12 ! user header (066) + real(DP) :: user13 ! user header (067) + real(DP) :: user14 ! user header (068) + real(DP) :: user15 ! user header (069) + real(DP) :: user16 ! user header (070) + integer :: iuser0 ! user header (098) + integer :: iuser1 ! user header (099) + integer :: iuser2 ! user header (100) + integer :: iuser3 ! user header (101) + integer :: iuser4 ! user header (102) + integer :: iuser5 ! user header (103) + integer :: iuser6 ! user header (104) + integer :: iuser7 ! user header (105) + logical :: luser0 ! user header (110) + + !! associated information from sac header + integer :: nzmonth ! month of begin time from nzjday + integer :: nzday ! day of begin time from nzjday + + integer :: tim ! absolute begin time from 1970/1/1 0:0:0 in second + + end type sac__hdr + !----------------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------------- + !> Write SAC file + !! + !! @par Usage + !! call sac__write( char filename, type__header, real data(:), logical sw ) + !! data can be single or double precisions + !! if sw = true, the existing file is automatically replaced. + !-- + interface sac__write + + module procedure wsac_d, wsac_s + + end interface + !---------------------------------------------------------------------------------------------- + +contains + + !---------------------------------------------------------------------------------------------- + !> Write SAC file + !! -- + subroutine wsac_d( fn_sac, ss, dat, overwrite ) + + character(*), intent(in) :: fn_sac + type(sac__hdr), intent(in) :: ss + real(DP), intent(in) :: dat(:) + logical, intent(in), optional :: overwrite + !-- + real(SP), allocatable :: fdat(:) + !---- + + allocate(fdat(1:ss%npts)) + fdat(1:ss%npts) = real(dat(1:ss%npts)) + + if( present( overwrite) ) then + call wsac_s( fn_sac, ss, fdat, overwrite) + else + call wsac_s( fn_sac, ss, fdat ) + end if + + deallocate( fdat ) + end subroutine wsac_d + !---------------------------------------------------------------------------------------------- + + !---------------------------------------------------------------------------------------------- + !> Write SAC file + subroutine wsac_s( fn_sac, ss, dat, overwrite ) + + character(*), intent(in) :: fn_sac + type(sac__hdr), intent(in) :: ss + real(SP), intent(in) :: dat(:) + logical, intent(in), optional :: overwrite + !-- + logical :: isexist + integer :: io + character(1) :: yn + !---- + + !! overwrite check + inquire( file = fn_sac, exist=isexist ) + if( isexist ) then + if( present( overwrite) ) then + if( .not. overwrite ) then + write(error_unit,*) 'wsac: file '//trim(fn_sac)//' exists.' + write(error_unit,*) 'wsac: could not overwrite the file.' + write(error_unit,*) 'wsac: return without success' + write(error_unit,*) + return + end if + else + write(error_unit,*) 'wsac: file '//trim(fn_sac)//' exists.' + write(error_unit,*) 'wsac: Overwrite ? (y/n)' + read(input_unit,'(A)') yn + if( yn /= 'y' .and. yn /='Y' ) then + write(error_unit,*) 'wsac: could not overwrite the file.' + write(error_unit,*) 'wsac: return without success' + write(error_unit,*) + return + end if + end if + end if + + open( newunit=io, file=trim(fn_sac), action='write', access='stream', form='unformatted') + + call sac__whdr(io, ss) + + write( io ) dat(1:ss%npts) + close( io ) + + end subroutine wsac_s + !---------------------------------------------------------------------------------------------- + + !---------------------------------------------------------------------------------------------- + !> Write SAC data header from pre-opened file io + !-- + subroutine sac__whdr(io, ss) + + integer, intent(in) :: io + type(sac__hdr), intent(in) :: ss + !-- + real(SP) :: fheader(70) + integer :: iheader(71:105) + logical :: lheader(106:110) + character(4) :: aheader(111:158) + integer :: i + !---- + + !! header initialize + fheader(1:70) = -12345.0 + iheader(71:105) = -12345 + lheader(106:110) = .false. + do i=111, 157, 2 + aheader( i ) = '-123' + aheader( i+1 ) = '45' + end do + + + ! Copy header data to temprary arrays + fheader( 1) = real( int( ss % delta * 1d7 ) ) / 1e7 + fheader( 2) = real( ss % depmin ) + fheader( 3) = real( ss % depmax ) + fheader( 6) = real( ss % b ) + fheader( 7) = real( ss % e ) + fheader( 8) = real( ss % o ) + fheader( 9) = real( ss % a ) + fheader( 11) = real( ss % t0 ) + fheader( 12) = real( ss % t1 ) + fheader( 13) = real( ss % t2 ) + fheader( 14) = real( ss % t3 ) + fheader( 15) = real( ss % t4 ) + fheader( 16) = real( ss % t5 ) + fheader( 17) = real( ss % t6 ) + fheader( 18) = real( ss % t7 ) + fheader( 19) = real( ss % t8 ) + fheader( 20) = real( ss % t9 ) + fheader( 32) = real( ss % stla ) + fheader( 33) = real( ss % stlo ) + fheader( 34) = real( ss % stel ) + fheader( 35) = real( ss % stdp ) + fheader( 36) = real( ss % evla ) + fheader( 37) = real( ss % evlo ) + fheader( 38) = real( ss % evel ) + fheader( 39) = real( ss % evdp ) + fheader( 40) = real( ss % mag ) + fheader( 41) = real( ss % user0 ) + fheader( 42) = real( ss % user1 ) + fheader( 43) = real( ss % user2 ) + fheader( 44) = real( ss % user3 ) + fheader( 45) = real( ss % user4 ) + fheader( 46) = real( ss % user5 ) + fheader( 47) = real( ss % user6 ) + fheader( 48) = real( ss % user7 ) + fheader( 49) = real( ss % user8 ) + fheader( 50) = real( ss % user9 ) + fheader( 51) = real( ss % dist ) + fheader( 52) = real( ss % az ) + fheader( 53) = real( ss % baz ) + fheader( 54) = real( ss % gcarc ) + fheader( 57) = real( ss % depmen ) + fheader( 58) = real( ss % cmpaz ) + fheader( 59) = real( ss % cmpinc ) + fheader( 64) = real( ss % user10 ) + fheader( 65) = real( ss % user11 ) + fheader( 66) = real( ss % user12 ) + fheader( 67) = real( ss % user13 ) + fheader( 68) = real( ss % user14 ) + fheader( 69) = real( ss % user15 ) + fheader( 70) = real( ss % user16 ) + + iheader( 71) = ss % nzyear + iheader( 72) = ss % nzjday + iheader( 73) = ss % nzhour + iheader( 74) = ss % nzmin + iheader( 75) = ss % nzsec + iheader( 76) = ss % nzmsec + iheader( 77) = ss % nvhdr + iheader( 80) = ss % npts + iheader( 86) = ss % iftype + iheader( 87) = ss % idep + iheader( 93) = ss % ievtyp + + iheader( 98) = ss % iuser0 + iheader( 99) = ss % iuser1 + iheader(100) = ss % iuser2 + iheader(101) = ss % iuser3 + iheader(102) = ss % iuser4 + iheader(103) = ss % iuser5 + iheader(104) = ss % iuser6 + iheader(105) = ss % iuser7 + + lheader(106) = ss % leven + lheader(107) = ss % lpspol + lheader(108) = ss % lovrok + lheader(109) = ss % lcalda + + lheader(110) = ss % luser0 + + aheader(111) = ss%kstnm(1:4); aheader(112) = ss%kstnm(5:8) + aheader(113) = ss%kevnm(1:4); aheader(114) = ss%kevnm(5:8) + aheader(115) = ss%kevnm(9:12); aheader(116) = ss%kevnm(13:16) + aheader(117) = ss%khole(1:4); aheader(118) = ss%khole(5:8) + aheader(119) = ss%ko(1:4); aheader(120) = ss%ko(5:8) + aheader(121) = ss%ka(1:4); aheader(122) = ss%ka(5:8) + aheader(123) = ss%kt0(1:4); aheader(124) = ss%kt0(5:8) + aheader(125) = ss%kt1(1:4); aheader(126) = ss%kt1(5:8) + aheader(127) = ss%kt2(1:4); aheader(128) = ss%kt2(5:8) + aheader(129) = ss%kt3(1:4); aheader(130) = ss%kt3(5:8) + aheader(131) = ss%kt4(1:4); aheader(132) = ss%kt4(5:8) + aheader(133) = ss%kt5(1:4); aheader(134) = ss%kt5(5:8) + aheader(135) = ss%kt6(1:4); aheader(136) = ss%kt6(5:8) + aheader(137) = ss%kt7(1:4); aheader(138) = ss%kt7(5:8) + aheader(139) = ss%kt8(1:4); aheader(140) = ss%kt8(5:8) + aheader(141) = ss%kt9(1:4); aheader(142) = ss%kt9(5:8) + aheader(143) = ss%kf(1:4); aheader(143) = ss%kf(5:8) + aheader(145) = ss%kuser0(1:4); aheader(146) = ss%kuser0(5:8) + aheader(147) = ss%kuser1(1:4); aheader(148) = ss%kuser1(5:8) + aheader(149) = ss%kuser2(1:4); aheader(150) = ss%kuser2(5:8) + aheader(151) = ss%kcmpnm(1:4); aheader(152) = ss%kcmpnm(5:8) + aheader(153) = ss%knetwk(1:4); aheader(154) = ss%knetwk(5:8) + aheader(155) = ss%kdatrd(1:4); aheader(156) = ss%kdatrd(5:8) + aheader(157) = ss%kinst(1:4); aheader(158) = ss%kinst(5:8) + + !! write + write( io ) fheader(1:70) + write( io ) iheader(71:105) + write( io ) lheader(106:110) + write( io ) aheader(111:158) + + end subroutine sac__whdr + !---------------------------------------------------------------------------------------------- + + !---------------------------------------------------------------------------------------------- + !> Initialize SAC header + subroutine sac__init( ss ) + + type(sac__hdr), intent(inout) :: ss + !-- + real(SP) :: ferr = -12345.0_SP + integer :: ierr = -12345 + character(6) :: cerr = '-12345' + !---- + + ss%delta = ferr + ss%depmin = ferr + ss%depmax = ferr + ss%b = ferr + ss%e = ferr + ss%o = ferr + ss%a = ferr + ss%t0 = ferr + ss%t1 = ferr + ss%t2 = ferr + ss%t3 = ferr + ss%t4 = ferr + ss%t5 = ferr + ss%t6 = ferr + ss%t7 = ferr + ss%t8 = ferr + ss%t9 = ferr + ss%stla = ferr + ss%stlo = ferr + ss%stel = ferr + ss%stdp = ferr + ss%evla = ferr + ss%evlo = ferr + ss%evel = ferr + ss%evdp = ferr + ss%mag = ferr + ss%user0 = ferr + ss%user1 = ferr + ss%user2 = ferr + ss%user3 = ferr + ss%user4 = ferr + ss%user5 = ferr + ss%user6 = ferr + ss%user7 = ferr + ss%user8 = ferr + ss%user9 = ferr + ss%dist = ferr + ss%az = ferr + ss%baz = ferr + ss%gcarc = ferr + ss%depmen = ferr + ss%cmpaz = ferr + ss%cmpinc = ferr + ss%nzyear = ierr + ss%nzjday = ierr + ss%nzhour = ierr + ss%nzmin = ierr + ss%nzsec = ierr + ss%nzmsec = ierr + ss%nvhdr = 6 + ss%npts = ierr + ss%iftype = 01 + ss%ievtyp = ierr + ss%idep = ierr + ss%leven = .true. + ss%lpspol = .false. + ss%lovrok = .true. + ss%lcalda = .true. + ss%kstnm = cerr + ss%kcmpnm = cerr + ss%kevnm = cerr + ss%khole = cerr + ss%ko = cerr + ss%ka = cerr + ss%kt0 = cerr + ss%kt1 = cerr + ss%kt2 = cerr + ss%kt3 = cerr + ss%kt4 = cerr + ss%kt5 = cerr + ss%kt6 = cerr + ss%kt7 = cerr + ss%kt8 = cerr + ss%kt9 = cerr + ss%kf = cerr + ss%kuser0 = cerr + ss%kuser1 = cerr + ss%kuser2 = cerr + ss%knetwk = cerr + ss%kdatrd = cerr + ss%kinst = cerr + + ss%nzmonth = ierr + ss%nzday = ierr + ss%tim = ierr + + !! inoficial headers + ss%user10 = ferr + ss%user11 = ferr + ss%user12 = ferr + ss%user13 = ferr + ss%user14 = ferr + ss%user15 = ferr + ss%user16 = ferr + ss%iuser0 = ierr + ss%iuser1 = ierr + ss%iuser2 = ierr + ss%iuser3 = ierr + ss%iuser4 = ierr + ss%iuser5 = ierr + ss%iuser6 = ierr + ss%iuser7 = ierr + ss%luser0 = .false. + + end subroutine sac__init + !---------------------------------------------------------------------------------------------- + +end module m_wsac +!------------------------------------------------------------------------------------------------ diff --git a/src/makefile b/src/makefile new file mode 100644 index 0000000..85e4c20 --- /dev/null +++ b/src/makefile @@ -0,0 +1,46 @@ +## Makefile for win package +## +## @copyright +## Copyright (c) 2019 Takuto Maeda. All rights reserved. +## +## @license +## This software is released under the MIT license. See LICENSE for details. + +FC = gfortran +FFLAGS = -pthread -ffast-math -Ofast -march=native +#FFLAGS = -fbounds-check -fbacktrace -O0 +all: ../fchinf.x ../fwck.x ../fdewin_a.x ../fdewin_s.x ../fwin2sac.x + +.SUFFIXES: +.SUFFIXES: .f90 .o + +.f90.o: + $(FC) -c $(FFLAGS) $< + +fchinf.o: fchinf.f90 m_win.o m_winch.o m_getopt.o m_util.o +fdewin_a.o: fdewin_a.f90 m_win.o m_getopt.o +fdewin_s.o: fdewin_s.f90 m_win.o m_getopt.o +fwin2sac.o: fwin2sac.f90 m_win.o m_getopt.o m_system.o m_util.o m_wsac.o m_util.o m_winch.o +fwck.o: fwck.f90 m_win.o +m_getopt.o: m_getopt.f90 m_system.o +m_winch.o: m_winch.f90 m_win.o +m_wsac.o: m_wsac.f90 m_util.o +m_win.o: m_win.f90 m_util.o + +../fchinf.x: fchinf.o m_win.o m_getopt.o m_system.o m_util.o m_winch.o m_util.o + $(FC) $(FFLAGS) $^ -o $@ + +../fwck.x: fwck.o m_win.o m_util.o + $(FC) $(FFLAGS) $^ -o $@ + +../fdewin_a.x: fdewin_a.o m_win.o m_getopt.o m_system.o m_util.o + $(FC) $(FFLAGS) $^ -o $@ + +../fdewin_s.x: fdewin_s.o m_win.o m_getopt.o m_system.o m_util.o + $(FC) $(FFLAGS) $^ -o $@ + +../fwin2sac.x: fwin2sac.o m_win.o m_getopt.o m_system.o m_util.o m_wsac.o m_util.o m_winch.o + $(FC) $(FFLAGS) $^ -o $@ + +clean: + /bin/rm -f *.mod *.o ../*.x