program wd2sac
!-----------------------------------------------------------------------------------------------
! waveform data with WD (extended WIN ; HAKUSAN ls8000wd) format is transformed to a SAC format
!   orignal version was written by R.Hino?
!   modified for gfortran, Y. Ishihara
!-----------------------------------------------------------------------------------------------

include 'wvprm.inc'
include 'wvinf.inc'
include 'wvdtc.inc'
include 'wvfln.inc'

character*8 stnd(nchdef)	! (i) reading station code
character*4 cmpd(nchdef)	! (i)          & component 
dimension xwv(nmemdef),ywv(nmemdef)
character*100 ch_file,ofile

open(51,file='wd2sac.prm',status='old')
read(51,*)ifile
read(51,*)ch_file
read(51,*)ndt
read(51,*)(dfile(i),i=1,ndt)
read(51,*)nchd		! number of channel to transfer

do ich=1,nchd
  read(51,*)stnd(ich),cmpd(ich)
enddo

close(51)

do ii=1,ndt
  call rdwv_wd( ifile,dfile,ii,ch_file,nch,wave,hp,st,wv,nchd,stnd,cmpd)
  
  do ich=1,nch
    write(ofile,'(a12,a1,a3,a1,a2)')dfile(ii),'.',st(ich)%code,'.',st(ich)%comp
	write(6,*)ofile(1:80)
	sampt=1./wv(ich)%sampf
	iyy=wv(ich)%yy+2000
	imm=wv(ich)%mm
	idd=wv(ich)%dd
	ihr=wv(ich)%hr
	imn=wv(ich)%mn
	isec=int(wv(ich)%sec)
	bs=wv(ich)%bsec
	write(*,*) 'bsec',bs
	msec=1000*(wv(ich)%sec-isec)
	call calijday(iyy,imm,idd,ijday)
	write(6,*)'ijday=',ijday
	write(6,*)'location : ',st(ich)%locx,st(ich)%locy,st(ich)%locz
	
	do i=1,wv(ich)%ndata
	  xwv(i)=wave(i,ich)*1.0e+3 ! e.g., *1.0e+3 for m/s ==> mm/s kPa ==> Pa
	enddo
	
	call wrsac (ofile,xwv,ywv,wv(ich)%ndata,sampt,iyy,ijday,ihr,imn,isec,msec, &
	            st(ich)%code,st(ich)%comp,st(ich)%locx,st(ich)%locy,st(ich)%locz)
  enddo  
enddo

stop
end



!---------------------------------------------------------------------------------------------------
subroutine wrsac(ofile,xwv,ywv,ndata,sampt,iyy,ijday,ihh,imi,isec,imsec,stnm,cmpn,slat,slon,sheight)
!---------------------------------------------------------------------------------------------------
! wrsac() make a SAC data file
!
real xwv(*),ywv(*)
character*(*) ofile
character*(*) stnm
character*(*) cmpn

bsec=0.
call newhdr                             ! create a new header
call setihv('IEVTYP','IUNKN',NERR)      ! set type of event
call setihv('IFTYPE','ITIME',NERR)      !          of file (ITIME:Time Series)
call setihv('IZTYPE','IB',NERR)         ! reference time equivalence (IB:Begin Time)
call setihv('IDEP','IUNKN',NERR)        ! type of dependent variable (IVEL:nm/sec, IUNKN:Unknown)
call setnhv('NZYEAR',iyy,  NERR)        !       year
call setnhv('NZJDAY',ijday,NERR)        !       jiulain day
call setnhv('NZHOUR',ihh,  NERR)        !       hour
call setnhv('NZMIN', imi,  NERR)        !       minute
call setnhv('NZSEC', isec, NERR)        !       second
call setnhv('NZMSEC',imsec,NERR)        !       mili-second
call setfhv('B',     bsec, NERR)        !
call setlhv('LEVEN',.TRUE.,NERR)        ! data is evenly spaced
call setfhv('DELTA',sampt,NERR)         ! increment between spaced
call setnhv('NPTS', ndata,NERR)         ! total data number
call setkhv('KSTNM', stnm,NERR)         ! station code
call setkhv('KCMPNM',cmpn,NERR)         ! station code
call setkhv('KINST','Model-2 ',NERR)     ! generic code of recordings
call setfhv('STLA',slat,NERR)           ! station location LATITUDE
call setfhv('STLO',slon,NERR)           !                  LONGITUDE
call setfhv('STEL',sheight,NERR)        !                  ELEVATION

call wsac0(ofile,ywv,xwv,NERR)

return
end



!-----------------------------------------------------------------------------
subroutine rdwv_wd( ifile,dfile,idt,ch_file,nch,wave,hp,st,wv,nchd,stnd,cmpd)
!-----------------------------------------------------------------------------
include 'wvprm.inc'
include 'wvinf.inc'
include 'wvdtc.inc'
include 'wvfln.inc'

character*8 stnd(nchdef)	! (i) reading station code
character*4 cmpd(nchdef)	! (i)              & component 

dimension jdate(6)
dimension iwave(nmemdef)		! waveform buffer 
character*100 wv_file,ch_file
character*50 jfile,cfile
character*128 chead

external win_open        !$pragma C(win_open)
external wd_read_head    !$pragma C(wd_read_head)
external win_read_1trace !$pragma C(win_read_1trace)
external win_close       !$pragma C(win_close)
	
iunit=51
ichr=0
ibufsiz=nmemdef

do i=1,nchd
  call lnfilenm(wv_file,ifile,dfile(idt))
  write(6,*)'ifile,dfile=',ifile,dfile(idt)
  
  call get_chinf1(ch_file,stnd(i),cmpd(i),ichno,pkand,db,vdig,delay,pickf,h,xlat,xlon,xhei,ier)
  write(6,*)idt,wv_file(1:60),stnd(i),cmpd(i),ier,xlat,xlon
  
  if (ier.eq.0) then
    lunit=50
	call win_open(lunit,wv_file,0,iflag)
	call wd_read_head(lunit,chead)	
	call win_read_1trace(lunit,ichno,ibufsiz,iwave,jdate,nsample,length,iflag)
	call win_close(lunit)
	
	if (iflag.ne.-1) then 
	  if (iflag.eq.-2) then
	    write(6,*)'--> Data array is not enough'
	  endif	  
	  	  
	  ichr=ichr+1
	  st(ichr)%code=stnd(i)
	  st(ichr)%comp=cmpd(i)
	  wv(ichr)%yy=jdate(1)
	  wv(ichr)%mm=jdate(2)
	  wv(ichr)%dd=jdate(3)
	  wv(ichr)%hr=jdate(4)
	  wv(ichr)%mn=jdate(5)
	  wv(ichr)%sec=float(jdate(6))+delay/1000.
	  call dslvchead(chead,wv,ichr)
	  wv(ichr)%bsec=wv(ichr)%sec
	  wv(ichr)%ndata=nsample*length
	  wv(ichr)%sampf=float(nsample)
	  wv(ichr)%gain=vdig/pkand/10**(db/20.)  
!	   wv(ichr)%unit='m/s'                         ! for velocity seismogram
      wv(ichr)%unit='Pa'                           ! for atmospheric pressure
	  st(ichr)%locx=xlat
	  st(ichr)%locy=xlon
	  st(ichr)%locz=xhei
	  write(6,*)jdate
	  
	  do j=1,wv(ichr)%ndata
	    wave(j,ichr)=float(iwave(j))*wv(ichr)%gain ! digital count ==> physical quantity
	  enddo
	  
	endif	
  endif	
enddo

nch=ichr

return	
end



!----------------------------------------------------------------------------------------------------
subroutine get_chinf1(ch_file, stn,comp, ichno, pkand, db, vdig,delay, pickf, h, xlat,xlon,xhei,ier)
!----------------------------------------------------------------------------------------------------

!ier     = -999  :       requested ch not exist

character ch_file*(*),stn*(*),comp*(*)
character cbuf*120,stnb*4,compb*2

iu=99
open(iu,file=ch_file,iostat=ier)

if(ier.ne.0) return
  icon=0
  
  do while(icon.eq.0)
    read(iu,'(a)') cbuf
	if(cbuf(1:1).ne.'#') icon=1
  end do
  
  backspace(iu)
  icon=0
  
  do while(icon.eq.0)
    read(iu,'(a)',iostat=ier) cbuf
	if(ier.eq.0) then
	  call pickinfo1(cbuf,ichno,stnb,compb,pkand,db,vdig,delay,pickf,h,xlat,xlon,xhei)
	  if(stnb.eq.stn(1:4).and.compb.eq.comp(1:2)) icon=1
	    else
		  icon=1
    end if
  end do
close(iu)

return
end



!--------------------------------------------------------------------------------------
subroutine pickinfo1(cbuf,ichno,stnb,compb,pkand,db,vdig,delay,pickf,h,xlat,xlon,xhei)
!--------------------------------------------------------------------------------------
character cbuf*(*),stnb*(*),compb*(*),cb*1,tab*1
       
read(cbuf(1:18),'(z4,i2,i4,1x,a4,1x,a2)')ichno,ii,idelay,stnb,compb

delay=idelay        
ier=0
cb=' '
tab='	'
i=19

do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) iamp
        
do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do
        
do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) iadbit

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do
        
do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) ipkand

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do

pkand=ipkand

do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),'(a3)') unit

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do

do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) pickf

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do

do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) h

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do

do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) db

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do
 
do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
  i=i+1
end do

read(cbuf(i:),*) vdig

do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
  i=i+1
end do

if (compb.eq.'V'.or.compb.eq.'UD'.or.compb.eq.'Pa')then       
  do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
    i=i+1
  end do
  
  read(cbuf(i:),*) xlat

  do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
    i=i+1
  end do
  
  do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
    i=i+1
  end do
  
  read(cbuf(i:),*) xlon
  
  do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
    i=i+1
  end do
  
  do while(cbuf(i:i).eq.cb.or.cbuf(i:i).eq.tab)
    i=i+1
  end do
  
  read(cbuf(i:),*) xhei 
  
  do while(cbuf(i:i).ne.cb.and.cbuf(i:i).ne.tab)
    i=i+1
  end do
endif

return
end



!---------------------------------------
subroutine dslvchead(chead,wv,ichr)
!---------------------------------------
include 'wvprm.inc'
include 'wvinf.inc'

character*128 chead

integer*2 it1,idt(8)
integer*4 it
character*2 ch
equivalence (it1,ch)
dimension its(8,4)

do ii=1,4
  do jj=1,4
    is=18+(jj-1)*2+(ii-1)*8
	ch=chead(is:is+1)
!	call chtrn(ch)
    
	if (it1.lt.0) then 
	  it=it1+256*256
	else
	  it=it1
	endif
	
	i1=int(it/16/256)
	i2=int((it-i1*16*256)/256)
	i3=int((it-i1*16*256-i2*256)/16)
	i4= it-i1*16*256-i2*256-i3*16
	k1=(jj-1)*2+1
	k2=k1+1
	its(k1,ii)=i1*10+i2
	its(k2,ii)=i3*10+i4
  enddo
enddo

do ii=1,4
  write(6,'(8i2.2)')(its(j,ii),j=1,8)
enddo

do i=1,8
  idt(i)=its(i,4)-its(i,2)
enddo

call diftime(idt,its(1,2),its(2,2),dtrn)

do i=1,8
  idt(i)=its(i,2)-its(i,1)
enddo

call diftime(idt,its(1,2),its(2,2),dt1)
	
do i=1,8
  idt(i)=its(i,4)-its(i,3)
enddo
	
call diftime(idt,its(1,2),its(2,2),dt2)

idt(1)=wv(ichr)%yy-its(1,2)
idt(2)=wv(ichr)%mm-its(2,2)
idt(3)=wv(ichr)%dd-its(3,2)
idt(4)=wv(ichr)%hr-its(4,2)
idt(5)=wv(ichr)%mn-its(5,2)
isec=int(wv(ichr)%sec)
isec1=wv(ichr)%sec*100-isec*100
isec2=wv(ichr)%sec*10000-isec1*100
idt(6)=isec-its(6,2)
idt(7)=isec1-its(7,2)
idt(8)=isec2-its(8,2)

call diftime(idt,its(1,2),its(2,2),dt3)

dt=(dt2-dt1)*dt3/dtrn+dt1
wv(ichr)%sec=wv(ichr)%sec+dt
write(6,*)'dtrn,dt1,dt2,dt3=',dtrn,dt1,dt2,dt3
write(6,*)'Delt t = ',dt

return
end	



!--------------------------
subroutine chtrn(ch)
!--------------------------
character*2 ch
character*1 ch1,ch2

ch1=ch(1:1)
ch2=ch(2:2)
ch=ch2//ch1

return
end



!------------------------------------------
subroutine diftime(idt,iyear,imonth,dtrn)
!------------------------------------------
integer*2 idt(8)

idy=365

if (mod(iyear,4).eq.0) then
  idy=366
endif

imn=imonth

if (imn.eq.2.and.idy.eq.365) then
  idd=28
 else
  idd=29
endif

if (imn.eq.4.or.imn.eq.6.or.imn.eq.9.or.imn.eq.11) then
  idd=30
 else
  idd=31
endif
	
idate=idt(1)*idy+idt(2)*idd+idt(1)
dtrn=idt(8)/10000.+idt(7)/100.+idt(6)+idt(5)*60.+idt(4)*3600+idate*24*3600.

return
end



!-----------------------------------------
subroutine calijday(iyy,imm,idd,ijday)
!-----------------------------------------
dimension imnth(12)
data (imnth(i),i=1,12) /31,28,31,30,31,30,31,31,30,31,30,31/

if (mod(iyy,4).eq.0) then
  imnth(2)=29
endif

i=1
ijday=0

do while (imm.gt.i)
  ijday=ijday+imnth(i)
  i=i+1
enddo
	
ijday=ijday+idd

return
end



!------------------------------------------
subroutine lnfilenm (file0,file1,file2)
!------------------------------------------
! file0=file1//file2 (excluded space)

character file0*(*),file1*(*),file2*(*)

lc1=1

do while(file1(lc1:lc1).eq.' ')
  lc1=lc1+1
enddo	
        
lc2=lc1+1

do while(file1(lc2:lc2).ne.' ')
  lc2=lc2+1
enddo

ld1=1

do while(file2(ld1:ld1).eq.' ')
  ld1=ld1+1
enddo

ld2=ld1+1

do while (file2(ld2:ld2).ne.' ')
  ld2=ld2+1
enddo 

file0=file1(lc1:lc2-1)//file2(ld1:ld2-1)

return
end