c********************************************************
c   XLDHG : Extract Line data                           *
c            from HGAM bird-mag and GPS position data   *
c               with diurnal correction by GSmag data,  *
c                    and calculation of IGRF-Residuals  *
c********************************************************
      real*8 tms, tme, tmsr, tmer, tmt
      real*8 uts, utsa, utsb, utss, utse, uts1, uts2
      character*80 wdr, flog, finf, fgst, fnam, fpos, fout
      character ss*60, line*150, str*80, name*8, ps*9, sdtm*20
      call premsg('XLDHG : Extract Line data')
      call premsg('         from HGAM bird-mag and GPS position data')
      call premsg('            with diurnal correction by GSmag data,')
      call premsg('                 and calculation of IGRF-Residuals')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      finf = wdr
      fgst = wdr
      fnam = wdr
      fpos = wdr
      fout = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Lines.Inf. filename ==> ',
     *            81-ldr, finf(ldr:80))
      open(10,file=finf,status='old')
      call premsg(' Pos. data specification:')
      call gparmi('   PNAV-DGPS (0) or Real-time (1) ==> ', ips)
      if (ips.eq.0) then
        nps = 80
        ps = 'PNAV-DGPS'
      else if (ips.eq.1) then
        nps = 84
        ps = 'Real-time'
      else
        call abendm('invalid value')
      endif
      call gparma('Enter GSmag data filename ==> ',
     *            81-ldr, fgst(ldr:80))
      open(30,file=fgst,status='old')
      call gparma('Enter Output filename ==> ', 81-ldr, fout(ldr:80))
      open(2,file=fout,status='new')
      call gparmf(' Year of survey  : ', year)
      call gparmi(' Generation number of IGRF : ', ngn)
      call gparmf(' Sensor separation (2nd sensor relative alt. in m)'
     *            // ' ==> ', dalt)
      call gparmi(' Thin-out Ratio N (crop 1 out of N) ==> ', intv)
      if (intv.le.0) intv = 1
      call clspin()
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'XLDHG : Extract Line data'
      write(6,'(8x,a)') ' from HGAM bird-mag and GPS position data'
      write(6,'(8x,a)') '    with diurnal correction by GSmag data,'
      write(6,'(8x,a)') '         and calculation of IGRF-Residuals'
      write(6,'(2a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf))
      write(6,'(2a)') '  Pos. data spec.: ', ps
      write(6,'(2a)') ' GSmag data file : ', fgst(1:lrtrim(fgst))
      write(6,'(2a)') ' Output filename : ', fout(1:lrtrim(fout))
      write(6,'(a,f7.2)') '   Year of survey  : ', year
      write(6,'(a,i2)')   '   IGRF generation : ', ngn
      write(6,'(a,f7.1)') '  2nd sensor relative Alt. (m) :', dalt
      write(6,'(a,i3)')   '  Data thin-out ratio :', intv
      call gigrf(ngn, year)
      open(9,status='scratch')
c
      base = 0.
      jymd = 0
      read(10,'(a)') str
      if (str(1:1).ne.'=') call abendm('input file not selected')
   40 is = 1
   41 is = is + 1
      if (is.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(is:is).eq.' ') goto 41
      it = is
   42 it = it + 1
      if (it.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(it:it).ne.' ') goto 42
      fnam(ldr:80) = str(is:it-1)
      is = it
   43 is = is + 1
      if (is.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(is:is).eq.' ') goto 43
      it = is
   44 it = it + 1
      if (it.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(it:it).ne.' ') goto 44
      fpos(ldr:80) = str(is:it-1)
      write(6,'(2a)') ' Input AM filename  : ', fnam(1:lrtrim(fnam))
      open(1,file=fnam,status='old')
      write(9,'(2a)') ' Input AM filename  : ', fnam(1:lrtrim(fnam))
      write(6,'(2a)') ' Pos. data filename : ', fpos(1:lrtrim(fpos))
      open(20,file=fpos,status='old')
      write(9,'(2a)') ' Pos. data filename : ', fpos(1:lrtrim(fpos))
      read(20,'(a)') ss
      k = index(ss, ':')
      if (k.eq.0 .or.
     *    ss(k+3:k+3).ne.':' .or.
     *    ss(k+6:k+6).ne.'.') call abendm('Pos. data illegal format')
      ss(k:k) = ' '
      ss(k+3:k+3) = ' '
      read(ss,*) ih, im, sec, fidoa, fkdoa, falta
      if (ih.lt.0.or.ih.ge.24 .or.
     *    im.lt.0.or.im.ge.60 .or.
     *    sec.lt.0..or.sec.ge.60.)
     *     call abendm('Pos. data illegal format')
      if (ih.ge.15) ih = ih - 24
      utsa = dble(ih*3600 + im*60) + sec
      read(20,'(a)') ss
      k = index(ss, ':')
      if (k.eq.0 .or.
     *    ss(k+3:k+3).ne.':' .or.
     *    ss(k+6:k+6).ne.'.') call abendm('Pos. data illegal format')
      ss(k:k) = ' '
      ss(k+3:k+3) = ' '
      read(ss,*) ih, im, sec, fidob, fkdob, faltb
      if (ih.ge.15) ih = ih - 24
      utsb = dble(ih*3600 + im*60) + sec
      if (utsb.le.utsa) call abendm('Pos. data UTC illegal sequence')
      read(1,'(1x,i8,1x,2(i2,1x),f5.2,2f10.2)')
     *     iymd, ih, im, sec, tmag1, tmag2
      uts = dble((ih-9)*3600 + im*60) + sec
      tmt = dble((ih*100+im)*100) + sec
      nrec = 1
c
   30 read(10,'(a)',end=3) str
      if (str(1:1).eq.'=') then
        close(1)
        close(20)
        goto 40
      endif
      is = 0
   31 is = is + 1
      if (is.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(is:is).eq.' ') goto 31
      it = is
   32 it = it + 1
      if (it.gt.80) call abendm('Lines.Inf. data illegal')
      if (str(it:it).ne.' ') goto 32
      name = str(is:it-1)
      read(str(it:80),*) tms, tme
      write(str,'(a8,2x,f9.2,1x,f9.2)') name, tms, tme
      call premsg('   Line: ' // str(1:29))
      write(2,'(a1,a)') '&', str(1:29)
      im = int(tms) / 100
      ih = im / 100
      utss = tms - dble(im*40 + ih*2400 + 9*3600)
      im = int(tme) / 100
      ih = im / 100
      utse = tme - dble(im*40 + ih*2400 + 9*3600)
      kst = 0
c
   10 if (uts.gt.utse) goto 14
      if (uts.lt.utsa) goto 13
      if (uts.lt.utss) goto 13
   11 if (uts.le.utsb) goto 12
      utsa = utsb
      fidoa = fidob
      fkdoa = fkdob
      falta = faltb
      read(20,'(a)',end=18) ss
      k = index(ss, ':')
      if (k.eq.0 .or.
     *    ss(k+3:k+3).ne.':' .or.
     *    ss(k+6:k+6).ne.'.') call abendm('Pos. data illegal format')
      ss(k:k) = ' '
      ss(k+3:k+3) = ' '
      read(ss,*) ih, im, sec, fidob, fkdob, faltb
      if (ih.ge.15) ih = ih - 24
      utsb = dble(ih*3600 + im*60) + sec
      if (utsb.le.utsa) call abendm('Pos. data UTC illegal sequence')
      goto 11
   12 if (kst.eq.0) then
        tmsr = tmt
        kst = 1
      endif
      tmer = tmt
      f = (uts-utsa) / (utsb-utsa)
      fido = fidoa + (fidob-fidoa)*f
      fkdo = fkdoa + (fkdob-fkdoa)*f
      falt = falta + (faltb-falta)*f
   20 if (iymd.ne.jymd) goto 21
      if (uts.lt.uts1) goto 29
      if (uts.gt.uts2) goto 25
      dv = (dm1*(uts2-uts) + dm2*(uts-uts1)) / (uts2-uts1)
      if (tmag1.ge.90000.) then
        tres1 = 9999.99
      else
        call igrfc(fido, fkdo, falt, tigrf)
        tmag1 = tmag1 - dv
        tres1 = tmag1 - tigrf
      endif
      if (tmag2.ge.90000.) then
        tres2 = 9999.99
      else
        call igrfc(fido, fkdo, falt+dalt, tigrf)
        tmag2 = tmag2 - dv
        tres2 = tmag2 - tigrf
      endif
      write(2,1000) nrec, iymd, tmt, nps, fido, fkdo, falt,
     *              tmag1, tres1, tmag2, tres2, (tmag1-tmag2)
   13 do 100 n=2,intv
        read(1,'(1x)',end=19)
  100 continue
      read(1,'(1x,i8,1x,2(i2,1x),f5.2,2f10.2)',end=19)
     *     iymd, ih, im, sec, tmag1, tmag2
      uts = dble((ih-9)*3600 + im*60) + sec
      tmt = dble((ih*100+im)*100) + sec
      nrec = nrec + 1
      goto 10
   14 if (kst.eq.0) goto 15
      write(str,'(a17,f9.2,1x,f9.2,a10,i4)')
     *     '    Extracted :  ', tmsr, wtmr, '    intv =', intv
      write(6,'(a)') str(1:50)
      write(9,'(a)') str(1:50)
      goto 30
   15 call premsg('    No data Extracted')
      goto 30
c
   18 call abendm('EOF reached on PNAV input')
   19 call abendm('EOF reached on HGAM input')
c
   21 read(30,'(a)',end=28) str
   22 if (str(1:5).ne.'/Base') goto 23
      read(str(7:40),*) base
      if (base.ge.200000.) base = base / 10.
      if (jymd.eq.iymd) goto 24
      goto 21
   23 if (str(1:5).ne.'/Date') goto 21
      read(str(7:40),*) jymd
      if (jymd.ne.iymd .or. base.eq.0.) goto 21
   24 read(30,'(a)',end=28) str
      if (str(1:1).eq.'/') goto 22
      read(str,'(3i2,1x,i6)') jh,jm,js, mg
      uts2 = float((jh-9)*3600 + jm*60 + js)
      jt2 = jh*10000 + jm*100 + js
      dm2 = float(mg)/10. - base
   25 uts1 = uts2
      jt1 = jh2
      dm1 = dm2
      read(30,'(a)',end=28) str
      if (str(1:1).eq.'/') goto 22
      read(str,'(3i2,1x,i6)') jh,jm,js, mg
      uts2 = float((jh-9)*3600 + jm*60 + js)
      jt2 = jh*10000 + jm*100 + js
      dm2 = float(mg)/10. - base
      goto 20
c
   28 call abendm('EOF on GSmag data file')
   29 write(str,'(i9,1x,f10.2,2x,i6,1x,i6)') iymd, tmt, jt1, jt2
      call abendm('out of GSmag data range : ' // str(1:35))
c
    3 close(2)
      close(10)
      close(1)
      close(20)
      endfile(9)
      rewind(9)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)') 'XLDHG : Extract Line data'
        write(99,'(8x,a)') ' from HGAM bird-mag and GPS position data'
        write(99,'(8x,a)') '    with diurnal correction by GSmag data,'
        write(99,'(8x,a)') '         and calculation of IGRF-Residuals'
        write(99,'(2a)') ' Lines.Inf. file : ', finf(1:lrtrim(finf))
        write(99,'(2a)') '  Pos. data spec.: ', ps
        write(99,'(2a)') ' GSmag data file : ', fgst(1:lrtrim(fgst))
        write(99,'(2a)') ' Output filename : ', fout(1:lrtrim(fout))
        write(99,'(a,f7.2)') '   Year of survey  : ', year
        write(99,'(a,i2)')   '   IGRF generation : ', ngn
        write(99,'(a,f7.1)') '  2nd sensor relative Alt. (m) :', dalt
        write(99,'(a,i3)')   '  Data thin-out ratio :', intv
   92   read(9,'(a)',end=93) line
        write(99,'(a)') line(1:lrtrim(line))
        goto 92
   93   call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      close(9)
      stop
c
 1000 format(i7, 1x,i8, 1x,f9.2, 1x,i2, 1x,f9.5, 1x,f10.5, 1x,f6.1,
     *       2(1x,f8.2), 3(1x,f8.2), 1x,f9.2)
      end
