c     file sac_io.f
c          ========
c
c     Revision by Karl Koch, 24Nov99
c
c     version 1, 20-Aug-98
c
c     SAC I/O routines for codeco.
c     Original routines by Hugues Dufumier
c     Major revision by K. Koch (19-Aug-98) and K. Stammler (20-Aug-98)



c==============================================================================



      subroutine sacin( infile, iy, ierror )
c
c     Program for conversions between GSE1.0, GSE2.0 and SAC data files.
c
c     Program modified by Hugues Dufumier to add Sac input/outputs.
c     Some bugs fixed and extension to more than 80 characters. January 98 /hd
c
c     ! Portability: Modify write(6,*) statements if not supported on HP or VMS
c                    and adapt subroutine INTAND to your computer !
c
c     Warning: Some bugs may persist with the compression schemes
c              GSE1.0 CMP7 & CMP8 and GSE2.0 CM8.
c              It is recommended to rather use CMP6 and CM6 compressions.
c              If not, for CMP7 & CMP8, prefer 0 or 1 differences to 2.

      implicit none

c     -- get header variables (hdr_...) and constants (c_...)
      include 'codeco_common.f'

c     -- parameters
      character*80   infile          ! name of input file (input)
      integer        iy(c_sigsize)   ! sample array (output)
      integer        ierror          ! return status (output)

c     -- local variables
      logical        sac_found       ! identified as SAC
      real           dt              ! sample distance in sec
      real           begin           ! start time of waveform segment
      real           sec             ! variable to hold seconds (real)
      integer*4 i, IRU, ierr
      real y(c_sigsize)

c     -- executable code

      if  (hdr_debug .gt. 2)  write(*,*)  '-- entering sacin'
      sac_found = .false.

c     -- test if this is a SAC binary file
      IRU=3
      call BRSAC(IRU,c_sigsize,infile,y,ierror)
      if  (ierror .eq. 0)  then
         if  (hdr_debug .gt. 0)  write(*,*)
     &      '-- detected SAC binary format.'
         hdr_ifmtname='SAC'
         hdr_ifmtnum=-1
         sac_found = .true.
      endif

c     -- test if this is a SAC ascii file
      if  (.not. sac_found)  then
         call ARSAC(IRU,c_sigsize,infile,y,ierror)
         if  (ierror .eq. 0)  then
            if  (hdr_debug .gt. 0)  write (*,*)
     &         '-- detected SAC ASCII format.'
            hdr_ifmtname='SACA'
            hdr_ifmtnum=-2
            sac_found = .true.
         endif
      endif

c     -- no SAC file, return
      if  (.not. sac_found)  then
         if  (hdr_debug .gt. 0)  write(*,*)
     &      '-- sacin: no sac file', infile
         hdr_ifmtname = 'UNK'
         hdr_ifmtnum = 0
         ierror = 999
         if  (hdr_debug .gt. 3)  write (*,*) '-- leaving sacin'
         return
      endif

c     -- extract number of samples and convert samples to integer
      CALL GETNHV('NPTS    ',hdr_nsamp,IERR)
      if  (hdr_debug .gt. 0)  write (*,*)
     &   '-- no.samples, err', hdr_nsamp,ierr
      do i=1,hdr_nsamp
         iy(i)=int(y(i))
      enddo

c     -- read Sac header values

      CALL GETFHV('DELTA   ',dt,ierr)
      CALL GETFHV('B       ',begin,ierr)
      hdr_smprate = 1.0 / dt

c     -- read the record time
      CALL GETNHV('NZYEAR  ',hdr_year,ierr)
      CALL GETNHV('NZJDAY  ',hdr_jday,ierr)
      CALL GETNHV('NZHOUR  ',hdr_hour,ierr)
      CALL GETNHV('NZMIN   ',hdr_min,ierr)
      CALL GETNHV('NZSEC   ',hdr_sec,ierr)
      CALL GETNHV('NZMSEC  ',hdr_msec,ierr)
      if  (hdr_debug .gt. 0)  write (*,*) '-- got time ',
     &   hdr_year, hdr_jday, hdr_hour, hdr_min, hdr_sec,
     &   hdr_msec, ierr

c     -- retrieve record time including begin shift
      sec=begin+float(hdr_sec)+float(hdr_msec)/1000.
      hdr_sec=int(sec)
      hdr_msec=1000.*(sec-float(hdr_sec))
      if (hdr_msec .lt.0) then
          hdr_sec=hdr_sec -1
c --- just in case you worry about the 999: the float to integer
c     operation from above cuts for negative values towards 0
          hdr_msec=999+hdr_msec
      endif
      hdr_min=hdr_min+int(float(hdr_sec)/60.)
      hdr_sec=mod(hdr_sec,60)
      if (hdr_sec .lt.0) then
          hdr_min=hdr_min-1
          hdr_sec=60+hdr_sec
      endif
      hdr_hour=hdr_hour+int(float(hdr_min)/60.)
      hdr_min=mod(hdr_min,60)
      if (hdr_min .lt.0) then
          hdr_hour=hdr_hour-1
          hdr_min=60+hdr_min
      endif
      hdr_jday=hdr_jday+int(float(hdr_hour)/24.)
      hdr_hour=mod(hdr_hour,24)
      if (hdr_hour .lt.0) then
          hdr_jday=hdr_jday-1
          hdr_hour=24+hdr_hour
      endif
      if (hdr_jday.lt.1) then
          hdr_year=hdr_year-1
          call julien( 31, 12, hdr_year, hdr_jday, 1 )
      else
          call julien( hdr_day, hdr_month, hdr_year, hdr_jday, 0 )
      endif
      if  (hdr_debug .gt. 0)  write (*,*) '-- got start time ',
     &   hdr_year, hdr_jday, hdr_hour, hdr_min, hdr_sec,
     &   hdr_msec, begin,sec

c     -- read other header values
      CALL GETKHV('KSTNM   ',hdr_station,ierr)
      CALL GETKHV('KCMPNM  ',hdr_chan,ierr)
      CALL GETKHV('KINST   ',hdr_instr,ierr)
      if  (hdr_debug .gt. 0)  write (*,*) '-- and ',
     &   hdr_station, hdr_chan, hdr_instr, ierr
      CALL GETFHV('STLA    ',hdr_stalat,ierr)
      CALL GETFHV('STLO    ',hdr_stalon,ierr)
      CALL GETFHV('STEL    ',hdr_staelev,ierr)
      CALL GETFHV('STDP    ',hdr_stadepth,ierr)
      if  (hdr_debug .gt. 0)  write (*,*) '-- and ',
     &   hdr_smprate, hdr_stalat, hdr_stalon, hdr_staelev,
     &   hdr_stadepth, ierr
      CALL GETFHV('BAZ     ',hdr_beamaz,ierr)
      CALL GETFHV('CMPAZ   ',hdr_hang,ierr)
      CALL GETFHV('CMPINC  ',hdr_vang,ierr)
c     calibration factor and period are usually not available from Sac
      hdr_calfac = 0.0
      hdr_calper = -1.0
      CALL GETFHV('USER0   ',hdr_calfac,IERR)
      CALL GETFHV('USER1   ',hdr_calper,IERR)
      if  (hdr_debug .gt. 0)  write (*,*) '-- and ',
     &   hdr_beamaz, hdr_hang, hdr_vang, hdr_calfac,
     &   hdr_calper, ierr
      if  (hdr_stalat .le. -12345.0)  hdr_stalat = -999.0
      if  (hdr_stalon .le. -12345.0)  hdr_stalon = -999.0
      if  (hdr_staelev .le. -12345.0)  hdr_staelev = -999.0
      if  (hdr_stadepth .le. -12345.0)  hdr_stadepth = -999.0
      if  (hdr_beamaz .le. -12345.0)  hdr_beamaz = -999.

cks      chid='        '
cks      hdr_idiff = 0

      if  (hdr_debug .gt. 3)  write (*,*) '-- leaving sacin'
      close(3)
      return
      end



c ==============================================================================



      subroutine sacout( outfile, iy, wavecnt )

c     Program to convert GSE2.0-datafiles into other GSE or Sac format.
c                        ------
c
c     Program adapted by Hugues Dufumier to add output in Sac format.
c     1997, EOST Ecole et Observatoire Sciences de la Terre, Strasbourg
c
c     Converts INT, CM2 and CM6 formats in INT, CM6 or SAC[A] ones.
c
c     Usage: gse2sac InputFile OutputFile OutputFormat (INT, CM6 or SAC[A])
c
c     For Sac outputs, all channels will be converted in single files.
c     Station and Component names will be appended: OutputFile.STA.CMP
c

      implicit none

c     -- get header variables
      include 'codeco_common.f'

c     -- parameters
      character      outfile*80      ! name of output file (input)
      integer        iy(c_sigsize)   ! sample array (input)
      integer        wavecnt         ! waveform counter (input)

c     -- local variables
      integer     ierr               ! return code
      integer     i                  ! counter
c     integer     y2                 ! 2-digit year
      real        dt                 ! sample distance in sec
      real        y(c_sigsize)       ! real array
      character FSAC*90
      integer*4 IRU, LN, NERR

c     -- functions
      integer   trimlen

c     -- executable code

      if  (hdr_debug .gt. 1)  write (*,*) '-- entering sacout'

      dt = 1.0 / hdr_smprate

c     -- build sac filename
      fsac = ' '
c     -- use station and time info in filename
c     y2 = hdr_year - 1900
c     if  (y2 .ge. 100)  y2 = y2 - 100
c     write( fsac, '(4a,3i2.2,a,3i2.2,2a)' )
c    &   outfile(1:trimlen(outfile)), '.',
c    &   hdr_station(1:trimlen(hdr_station)), '_', y2, hdr_month,
c    &   hdr_day, '_', hdr_hour, hdr_min, hdr_sec, '.',
c    &   hdr_chan(1:trimlen(hdr_chan))
c     -- just number files
      write( fsac, '(2a,i3.3)' )
     &   outfile(1:trimlen(outfile)), '.', wavecnt
      if  (hdr_debug .gt. 2)  write (*,*) '-- create sac file ',
     &   outfile(1:trimlen(outfile)), wavecnt

c     -- writes in header the useful values
      CALL NEWHDR()
      CALL SETNHV('NPTS    ',hdr_nsamp,IERR)
      CALL SETLHV('LEVEN   ',.TRUE.,IERR)
      CALL SETLHV('LOVROK  ',.TRUE.,IERR)
      CALL SETLHV('LCALDA  ',.TRUE.,IERR)
      CALL SETFHV('B       ',0.,IERR)
      CALL SETFHV('E       ',FLOAT(hdr_nsamp-1)*dt,IERR)
      CALL SETFHV('DELTA   ',dt,IERR)
      CALL SETIHV('IFTYPE  ','ITIME',IERR)
      CALL SETIHV('IZTYPE  ','IB',IERR)
c     -- writes the record time
      CALL SETNHV('NZYEAR  ',hdr_year,IERR)
      CALL SETNHV('NZJDAY  ',hdr_jday,IERR)
      CALL SETNHV('NZHOUR  ',hdr_hour,IERR)
      CALL SETNHV('NZMIN   ',hdr_min,IERR)
      CALL SETNHV('NZSEC   ',hdr_sec,IERR)
      CALL SETNHV('NZMSEC  ',hdr_msec,IERR)
c     -- writes other header values
      CALL SETKHV('KSTNM   ',hdr_station,IERR)
      CALL SETKHV('KCMPNM  ',hdr_chan,IERR)
      CALL SETKHV('KINST   ',hdr_instr,IERR)
      CALL SETKHV('KUSER0  ','codeco',IERR)
      CALL SETIHV('ISYNTH  ','IRLDTA',IERR)
      CALL SETFHV('CMPAZ   ',hdr_beamaz,IERR)
      CALL SETFHV('CMPINC  ',hdr_beamslo,IERR)
c     -- calibration factor and period reported in user0 and user1
c     -- change to displacement units
      hdr_calfac = hdr_calfac / 
     &   (atan(1.0)*8.0/hdr_calper)**hdr_calunit
      hdr_calunit = 0
      CALL SETFHV('USER0   ',hdr_calfac,IERR)
      CALL SETFHV('USER1   ',hdr_calper,IERR)
      call setihv('IFTYPE   ','ITIME   ',nerr)
      call setihv('IZTYPE   ','IB      ',nerr)

c     -- set logical header value
      call setlhv('LEVEN   ',.true.,nerr)
      call setlhv('LSPOL   ',.true.,nerr)
      call setlhv('LOVERK  ',.false.,nerr)
      call setlhv('LCALDA  ',.false.,nerr)
c     -- writes the data
      do i=1,hdr_nsamp
         y(i)=float(iy(i))
      enddo
C     CALL WSAC0(FSAC,Y,Y,IERR)
      IRU = 1
      LN = c_sigsize
      if (hdr_ofmtname.eq.'SACA') then
          call awsac(IRU,LN,FSAC(1:trimlen(FSAC)),Y)
      else
          call bwsac(IRU,LN,FSAC(1:trimlen(FSAC)),Y)
      endif
      IF(IERR.EQ.0) WRITE(6,'(3A,I6,A,F9.6,A)')'Written ',
     &   FSAC(1:trimlen(FSAC)),' :',hdr_nsamp,' points, dt=',dt,' s.'

      if  (hdr_debug .gt. 1)  write (*,*) '-- leaving sacout'
      return 
      end


