!**********************************************************
!   XLDHG : Extract Line data                             *
!            from HGAM bird-mag and GPS position data     *
!               with diurnal correction by GSmag data,    *
!                    and calculation of IGRF-Residuals    *
!   This program assumes no survey flight over midnight.  *
!**********************************************************
  implicit none
  integer :: lwkdir, lrtrim
  character(80) :: wdr, flog, finf, fgst, fnam, fpos, fout
  character(150) :: line;  character(80) :: str;  character(60) :: ss
  character(20) :: sdtm;  character(9) :: ps;  character(8) :: name
  character(5) :: tzone
  integer :: ldr, ips, nps, ngn, intv, kam, is, it, kps, kst, IER
  integer :: k, ih, im, nrec, iymds, iymd, kt, mintz, kintv
  real :: year, dalt, sec, falt, falta, faltb
  real :: dvc, f, tmag1, tmag2, tres1, tres2, tigrf
  real(8) :: wtm, tms, tme, tmsr, tmer
  real(8) :: sect, seca, secb, secst, seced
  real(8) :: dido, dkdo, didoa, dkdoa, didob, dkdob

  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(' Position data specification ?')
  call gparmi('   PNAV-DGPS (0) or Real-time (1) ==> ', ips)
  if (ips == 0) then
    nps = 80;  ps = 'PNAV-DGPS'
  else if (ips == 1) then
    nps = 84;  ps = 'Real-time'
  else
    call abendm('invalid value')
  endif
  call premsg(' Time Spec. of Position data ?')
  call gparmi('  Localtime sama as Obs.Mag.data [0] or UTC[1] ==> ', kt)
  if (kt == 0) then
    mintz = 0
  else if (kt == 1) then
    call premsg('    UTC time will be converted into Localtime.')
    call gparma('    Localtime-Zone ?  "[+/-]HHMM" (default: +0900) ==> ',  &
                5, tzone)
    if (tzone == '     ') tzone = '+0900'
    if ((tzone(1:1) /= '+') .and. (tzone(1:1) /= '-'))  &
      call abendm('invalid timezone specification')
    read(tzone(2:5),'(2i2)') ih, im
    if ((ih<0) .or. (ih>23) .or. (im<0) .or. (im>59))  &
      call abendm('invalid timezone specification')
    mintz = ih*60 + im
    if (tzone(1:1) == '-') mintz = -mintz
  else
    call abendm('invalid specification')
  endif
  call gparma('Enter GSmag data filename ==> ', 81-ldr, fgst(ldr:80))
  open(20,file=fgst,status='old')
  call gparma('Enter Output filename ==> ', 81-ldr, fout(ldr:80))
  open(30,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 <= 0) intv = 1
  call clspin()
  call dvcmag(0, 0.d0, dvc)
  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
  if (kt == 0) then
    write(6,'(a)') '  Input Pos. data Time : Localtime (same as Obs.Mag.)'
  else
    write(6,'(a)')   '  Input Pos. data Time : UTC'
    write(6,'(a,a)') '    will be converted into Localtime of Zone=', tzone
  endif
  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')
  kam = 0; kps = 0
  do
    read(10,'(a)',iostat=IER) str
    if (IER < 0) exit
    if (str(1:1) == '=') then
      is = 2
      do while ((is <= 80) .and. (str(is:is) == ' ')); is = is + 1; enddo
      if (is > 80) call abendm('Lines.Inf. data illegal')
      it = is + 1
      do while ((it <= 80) .and. (str(it:it) /= ' ')); it = it + 1; enddo
      if (it > 80) call abendm('Lines.Inf. data illegal')
      fnam(ldr:80) = str(is:it-1)
      is = it + 1
      do while ((is <= 80) .and. (str(is:is) == ' ')); is = is + 1; enddo
      if (is > 80) call abendm('Lines.Inf. data illegal')
      it = is + 1
      do while ((it <= 80) .and. (str(it:it) /= ' ')); it = it + 1; enddo
      if (it > 80) call abendm('Lines.Inf. data illegal')
      fpos(ldr:80) = str(is:it-1)
      if (kam == 1) close(1);  kam = 1
      write(6,'(2a)') ' Input AM filename  : ', fnam(1:lrtrim(fnam))
      write(9,'(2a)') ' Input AM filename  : ', fnam(1:lrtrim(fnam))
      open(1,file=fnam,status='old')
      read(1,'(1x,i8,1x,2(i2,1x),f5.2,2f10.2)',iostat=IER)   &
        iymds, ih, im, sec, tmag1, tmag2            ! get first HGAM data
      nrec = 1
      if (IER < 0) call abendm('EOF reached on HGAM input')
      wtm = dble(ih*10000 + im*100) + sec
      sect = dble(ih*3600 + im*60.) + sec
      if (kps == 1) close(2); kps = 1
      write(6,'(2a)') ' HGAM data filename : ', fpos(1:lrtrim(fpos))
      write(9,'(2a)') ' HGAM data filename : ', fpos(1:lrtrim(fpos))
      open(2,file=fpos,status='old')
      read(2,'(a)',iostat=IER) ss                   ! get first HPOS data
      if (IER < 0) call abendm('EOF reached on HPOS input')
      k = index(ss, ':')
      if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.'))   &
        call abendm('HPOS data illegal format')
      ss(k:k) = ' ';  ss(k+3:k+3) = ' '
      read(ss,*) ih, im, sec, didoa, dkdoa, falta
      if ((ih < 0) .or. (im < 0) .or. (sec < 0.) .or.   &
          (ih >= 24) .or. (im >= 60) .or. (sec >= 60.))   &
        call abendm('HPOS data illegal format')
      seca = dble((ih*60 + im + mintz) * 60) + sec
      if (seca >= 86400.) seca = seca - 86400.
      if (seca < 0.) seca = seca + 86400.
      read(2,'(a)',iostat=IER) ss                   ! get next HPOS data
      if (IER < 0) call abendm('EOF reached on HPOS input')
      k = index(ss, ':')
      if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.'))   &
        call abendm('HPOS data illegal format')
      ss(k:k) = ' ';  ss(k+3:k+3) = ' '
      read(ss,*) ih, im, sec, didob, dkdob, faltb
      secb = dble((ih*60 + im + mintz) * 60) + sec
      if (secb >= 86400.) secb = secb - 86400.
      if (secb < 0.) secb = secb + 86400.
      if (secb <= seca) call abendm('HPOS data [sec] illegal sequence')
      cycle
    endif
    if (kam == 0) call abendm('input file not selected')
    is = 1
    do while ((is <= 80) .and. (str(is:is) == ' ')); is = is + 1; enddo
    if (is > 80) call abendm('Lines.Inf. data illegal')
    it = is + 1
    do while ((it <= 80) .and. (str(it:it) /= ' ')); it = it + 1; enddo
    if (it > 80) call abendm('Lines.Inf. data illegal')
    name = str(is:it-1)
    read(str(it:80),*) tms, tme
    write(str,'(a8,1x,i9,1x,f9.2,1x,f9.2)') name, iymd, tms, tme
    call premsg('   Line: ' // str(1:38))
    write(30,'(a1,a)') '&', str(1:38)
    if ((tms < 0.) .or. (tme < 0.) .or. (tme <= tms))   &
      call abendm('illegal TIME spec')
    im = int(tms/100.); sec = sngl(tms - im*100.)
    ih = im/100; im = im - ih*100
    if ((ih > 23) .or. (im > 59) .or. (sec >= 60.))   &
      call abendm('illegal TIME spec')
    secst = dble(ih*3600 + im*60) + sec
    im = int(tme/100.); sec = sngl(tme - im*100.)
    ih = im/100; im = im - ih*100
    if ((ih > 23) .or. (im > 59) .or. (sec >= 60.))   &
      call abendm('illegal TIME spec')
    seced = dble(ih*3600 + im*60) + sec
    kst = 0; kintv = 1
    do while (sect < secst)
      read(1,'(1x,i8,1x,2(i2,1x),f5.2,2f10.2)',iostat=IER)   &
        iymd, ih, im, sec, tmag1, tmag2             ! get next HGAM data
      if (IER < 0) call abendm('EOF reached on AM input')
      if (iymd /= iymds) call abendm('date change during a file')
      wtm = dble(ih*10000 + im*100) + sec
      sect = dble(ih*3600 + im*60) + sec
      nrec = nrec + 1
    enddo
    if (sect > seced) call abendm('Nodata Extracted')
    do while (wtm <= tme)
      if (seca > sect) then
        write(ss,'(a10,f9.1,a26,f9.1)')  &
          'HPOS (sec=', seca, ') passed HGAM time [sec=]', sect
        call abendm(ss)
      endif
      do while (secb < sect)
        seca = secb;  didoa = didob;  dkdoa = dkdob;  falta = faltb
        read(2,'(a)',iostat=IER) ss
        if (IER < 0) call abendm('EOF reached on HPOS data input')
        k = index(ss, ':')
        if ((k == 0) .or. (ss(k+3:k+3) /= ':') .or. (ss(k+6:k+6) /= '.'))   &
          call abendm('HPOS data illegal format')
        ss(k:k) = ' ';  ss(k+3:k+3) = ' '
        read(ss,*) ih, im, sec, didob, dkdob, faltb
        secb = dble((ih*60 + im + mintz) * 60) + sec
        if (secb >= 86400.) secb = secb - 86400.
        if (secb < 0.) secb = secb + 86400.
        if (secb <= seca) then
          write(ss,'(a10,f9.1,a21,f9.1)')  &
            'HPOS (sec=', secb, ') illegal seq. after ', seca
          call abendm(ss)
        endif
      enddo
      kintv = kintv - 1
      if (kintv == 0) then
        f = real((sect-seca) / (secb-seca))
        dido = didoa + (didob-didoa)*f;  dkdo = dkdoa + (dkdob-dkdoa)*f
        falt = falta + (faltb-falta)*f
        call dvcmag(iymd, sect, dvc)
        if (tmag1 >= 90000.) then
          tres1 = 9999.9
        else
          call igrfc(sngl(dido), sngl(dkdo), falt, tigrf)
          tmag1 = tmag1 - dvc;  tres1 = tmag1 - tigrf
        endif
        if (tmag2 >= 90000.) then
          tres2 = 9999.9
        else
          call igrfc(sngl(dido), sngl(dkdo), falt+dalt, tigrf)
          tmag2 = tmag2 - dvc;  tres2 = tmag2 - tigrf
        endif
        write(30,'(2(i8,1x),f9.2,1x,i2,1x,f11.7,1x,f12.7,1x,f7.2,$)')   &
          nrec, iymd, wtm, nps, dido, dkdo, falt
        write(30,'(5(1x,f8.2))') tmag1, tres1, tmag2, tres2, (tmag1-tmag2)
        tmer = wtm;  kintv = intv
        if (kst == 0) then
          tmsr = wtm;  kst = 1
        endif
      endif
      read(1,'(1x,i8,1x,2(i2,1x),f5.2,2f10.2)',iostat=IER)   &
        iymd, ih, im, sec, tmag1, tmag2             ! get next HGAM data
      if (IER < 0) call abendm('EOF reached on HGAM input')
      if (iymd /= iymds) call abendm('date change during a file')
      wtm = dble(ih*10000 + im*100) + sec
      sect = dble(ih*3600 + im*60) + sec
      nrec = nrec + 1
    enddo
    write(str,'(a17,f9.2,1x,f9.2,a10,i4)')   &
      '    Extracted :  ', tmsr, tmer, '    intv =', intv
    write(6,'(a)') str(1:50)
    write(9,'(a)') str(1:50)
  enddo
  close(1);  close(2);  close(10);  close(20);  close(30)
  endfile(9)
  rewind(9)

  if (flog(ldr:ldr) /= ' ') 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
    if (kt == 0) then
      write(99,'(a)') '  Input Pos. data Time : Localtime (same as Obs.Mag.)'
    else
      write(99,'(a)')   '  Input Pos. data Time : UTC'
      write(99,'(a,a)') '    will be converted into Localtime of Zone=', tzone
    endif
    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
    do
      read(9,'(a)',iostat=IER) line
      if (IER < 0) exit
      write(99,'(a)') line(1:lrtrim(line))
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end
!---
  subroutine dvcmag(iymd, sect, dvc)
  implicit none
  integer,intent(in) :: iymd
  real(8),intent(in) :: sect
  real,intent(out) :: dvc
  save
  character(40) :: buf; character(17) :: dtm
  integer :: jymd, jymdn, jhms, jh, jm, js, kgs=0, IER
  real :: base, vm, dvc1, dvc2
  real(8) :: secn, sec1, sec2

  if (iymd == 0) then
    jymd = 0;  base = 0.
    do
      read(20,'(a)',iostat=IER) buf
      if (IER < 0) call abendm('EOF on GSmag data file')
      if (buf(1:1) == '/') then
        if ((buf(1:5) == '/date') .or. (buf(1:5) == '/Date')) then
          read(buf(7:40),*) jymdn
          if (jymdn < jymd) call abendm('GSmag data out of order')
          jymd = jymdn;  kgs = 0;  cycle
        else if ((buf(1:5) == '/base') .or. (buf(1:5) == '/Base')) then
          read(buf(7:40),*) base
          do while (base > 200000.)
            base = base / 10.
          enddo
          cycle
        endif
      else
        read(buf,*) jhms, vm
        jm = jhms / 100;  js = jhms - jm*100
        jh = jm / 100;  jm = jm - jh*100
        secn = dble(jh*3600 + jm*60 + js)
        do while (vm > 200000.)
          vm = vm / 10.
        enddo
        if (kgs == 0) then
          sec2 = secn;  dvc2 = vm - base;  kgs = 1;  cycle
        else
          sec1 = sec2;  dvc1 = dvc2
          sec2 = secn;  dvc2 = vm - base;  kgs = 2;  exit
        endif
      endif
    enddo
    dvc = 0.; return
  endif
  do while (jymd < iymd)
    do
      read(20,'(a)',iostat=IER) buf
      if (IER < 0) call abendm('EOF on GSmag data file')
      if (buf(1:1) == '/') then
        if ((buf(1:5) == '/date') .or. (buf(1:5) == '/Date')) then
          read(buf(7:40),*) jymdn
          if (jymdn < jymd) call abendm('GSmag data out of order')
          jymd = jymdn;  kgs = 0;  cycle
        else if ((buf(1:5) == '/base') .or. (buf(1:5) == '/Base')) then
          read(buf(7:40),*) base
          do while (base > 200000.)
            base = base / 10.
          enddo
          cycle
        endif
      else
        read(buf,*) jhms, vm
        jm = jhms / 100;  js = jhms - jm*100
        jh = jm / 100;  jm = jm - jh*100
        secn = dble(jh*3600 + jm*60 + js)
        do while (vm > 200000.)
          vm = vm / 10.
        enddo
        if (kgs == 0) then
          sec2 = secn;  dvc2 = vm - base;  kgs = 1;  cycle
        else
          sec1 = sec2;  dvc1 = dvc2
          sec2 = secn;  dvc2 = vm - base;  kgs = 2;  exit
        endif
      endif
    enddo
  enddo
  do while (sect > sec2)
    do
      if ((jymd > iymd) .or. (sec1 > sect)) then
        jhms = idnint(sect);  jm = jhms / 60;
        js = jhms - jm*60;  jh = jm / 60;  jm = jm - jh*60
        write(dtm,'(i8,a1,i2.2,a1,i2.2,a1,i2.2)') iymd, '-', jh,':',jm,':',js
        call abendm('GSmag data not available for ' // dtm)
      endif
      read(20,'(a)',iostat=IER) buf
      if (IER < 0) call abendm('EOF on GSmag data file')
      if (buf(1:1) == '/') then
        if ((buf(1:5) == '/date') .or. (buf(1:5) == '/Date')) then
          read(buf(7:40),*) jymdn
          if (jymdn < jymd) call abendm('GSmag data out of order')
          jymd = jymdn;  kgs = 0;  cycle
        else if ((buf(1:5) == '/base') .or. (buf(1:5) == '/Base')) then
          read(buf(7:40),*) base
          do while (base > 200000.)
            base = base / 10.
          enddo
          cycle
        endif
      else
        read(buf,*) jhms, vm
        jm = jhms / 100;  js = jhms - jm*100
        jh = jm / 100;  jm = jm - jh*100
        secn = dble(jh*3600 + jm*60 + js)
        do while (vm > 200000.)
          vm = vm / 10.
        enddo
        if (kgs == 0) then
          sec2 = secn;  dvc2 = vm - base;  kgs = 1;  cycle
        else
          sec1 = sec2;  dvc1 = dvc2
          sec2 = secn;  dvc2 = vm - base;  kgs = 2;  exit
        endif
      endif
    enddo
  enddo
  dvc = real((dvc1*(sec2-sect) + dvc2*(sect-sec1)) / (sec2-sec1))
  return
  end subroutine