!--------( XSLINA : extract StdLIN data from DPAM data      --------
!--------                        with averaged re-sampling )--------
  implicit none
  integer :: lwkdir, lrtrim
  real :: hgeoid
  character(80) :: wdr, flog, fnam1, fnam2
  character(150) :: buf;  character(20) :: sdtm, sgmdl, sgmdln
  character(24) :: ca, can;  character(3) :: cg, cgn
  integer :: ldr, kg, ka, kgn, kan, kr, nd, iymd, km, l, IER
  integer :: ihm, ihsec, jhsec
  real :: tint, sec, fid
  real(8) :: wtm, slat, slon, salt, smag, rlat, rlon, ralt, rmag
  real(8) :: vlat, vlon, valt, tlat, tlon, talt, wlat, wlon, walt, tmag, vmag

  call premsg('XSLINA : extract StdLIN data from DPAM data '   &
              // 'with averaged re-sampling')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input AM line data filename ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  call gparmi(' Geodetic ?  in WGS(1) or TokyoDatum(2) ==> ', kg)
  if (kg == 1) then
    cg = 'WGS'
    call gparmi(' Alt. ?  above Ellipsoid(1) or Geoid(2) ==> ', ka)
    if (ka == 1) then
      ca = 'above ITRF-GRS'
    else if (ka == 2) then
      sgmdl = '                    '
      call gparma('  Geoid model ?  0:NGA, or others [Enter NAME]  ==> ',   &
                  20, sgmdl)
      if (sgmdl == '0') sgmdl = 'NGA'
      call sgeoid(sgmdl)
      l = lrtrim(sgmdl)
      if (l <= 12) then
        ca = 'above Geoid ' // sgmdl(1:l)
      else
        ca = 'above Geoid ' // sgmdl(1:9) // '...'
      endif
    else
      call abendm('Alt. reference unknown')
    endif
  else if (kg == 2) then
    cg = 'TYO'
    call premsg(' Input Alt. data is treated as "above Geoid NGA"')
    ca = 'above Geoid NGA';  sgmdl = 'NGA';  call sgeoid(sgmdl)
  else
    call abendm('Geodetic reference unknown')
  endif

  call gparmf('Averaging interval (in sec.) ==> ', tint)
  if ((tint < 0.1) .or. (tint > 5.)) call abendm('invalid value')
  call gparma('Enter Output StdLIN data filename ==> ', 81-ldr, fnam2(ldr:80))
  open(2,file=fnam2,status='new')
  call gparmi(' Geodetic ?  in WGS(1) or TokyoDatum(2) ==> ', kgn)
  if (kgn == 1) then
    cgn = 'WGS'
    call prompt(' Output Alt. ?  above Ellipsoid(1) or Geoid(2)')
    call gparmi(' ==> ', kan)
    if (kan.eq.1) then
      can = 'above ITRF-GRS  '
    else if (kan.eq.2) then
      sgmdln = '                    '
      call gparma('  Geoid model ?  0:NGA, or others [Enter NAME]  ==> ',   &
                  20, sgmdln)
      if (sgmdln.eq.'0') sgmdln = 'NGA'
      call sgeoid(sgmdln)
      l = lrtrim(sgmdln)
      if (l.le.12) then
        can = 'above Geoid ' // sgmdln(1:l)
      else
        can = 'above Geoid ' // sgmdln(1:9) // '...'
      endif
    else
      call abendm('Alt. reference unknown')
    endif
  else if (kgn == 2) then
    call premsg(' Output Alt. data is "above Geoid NGA"')
    can = 'above Geoid NGA';  sgmdln = 'NGA';  call sgeoid(sgmdln)
  else
    call abendm('Geodetic reference unknown')
  endif
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'XSLINA : extract StdLIN data from DPAM data '   &
                 // 'with averaged re-sampling'
  write(6,'(a,a)') ' DPAM line data filename : ', fnam1(1:lrtrim(fnam1))
  write(6,'(8x,4a)')   &
     ' 5: ', cg, '-Lat.(Deg.)    7: Alt.(m) ', ca,   &
     ' 6: ', cg, '-Lon.(Deg.)    9: Mag.(nT)', ''
  write(6,'(a,f6.2,a)') '  Averaging interval :', tint, ' sec.'
  write(6,'(a,a)') ' Output StdLIN data filename : ', fnam2(1:lrtrim(fnam2))
  write(6,'(8x,4a)') 'Lat./Lon. in ', cgn, ';  Alt. ', can
  write(2,'(a)')  '# output from XSLIN'
  write(2,'(4a)') '#   Lat./Lon. in ', cgn, ';  Alt. ', can
  write(2,'(a,a)') '# org. line data filename : ', fnam1(1:lrtrim(fnam1))

  kr = 0            ! kr==1 if averaged re-sampling activated
                    !  re-sampled data is determined as the mid point of
                    !   2 consecutive areraging time intervals
  nd = 0            ! nd : number of data in averaging
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1)=='&') .or. (buf(1:1)=='%') .or. (buf(1:1)=='#')) then
      if ((nd /= 0) .and. (kr == 1)) then   ! output last data of prev.line
        slat = slat / float(nd);  slon = slon / float(nd)
        salt = salt / float(nd);  smag = smag / float(nd)
        write(2,'(1x,f11.5,a2,f11.5,a2,f8.2,a2,f8.2,a2)')   &
          (rlat+slat)/2., 'N ', (rlon+slon)/2., 'E ',   &
          (ralt+salt)/2., 'm ', (rmag+smag)/2., 'nT'
      endif
      kr = 0;  nd = 0
      call premsg(buf(1:lrtrim(buf)))
      write(2,'(a)') buf(1:lrtrim(buf))
      cycle
    endif
    read(buf,*) fid, iymd, wtm, km, vlat, vlon, valt, tmag, vmag
    ihm = int(wtm) / 100;  sec = real(wtm - dble(ihm*100))
    ihm = ihm - ihm/100*40;  sec = sec + float(ihm*60)
    ihsec = int(sec/tint);  vlat = vlat * 60.;  vlon = vlon * 60.
!           [Let vlat/vlon/valt be in Output side spec]
    if (kg == 1) then
      if (kgn == 1) then
!               (vlat/vlon are in WGS, common IN/OUT)
        if (ka == 1) then
          if (kan == 1) then
!               [valt is above WGS Ellips, common IN/OUT]
          else
            valt = valt - hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above WGS Geoid (conv. from WGS Ellips)]
          endif
        else
          if (kan == 1) then
            valt = valt + hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above WGS Ellips (conv. from WGS geoid)]
          else
            if (sgmdl.ne.sgmdln) then
              call sgeoid(sgmdl)
              valt = valt + hgeoid(sngl(vlat), sngl(vlon))
              call sgeoid(sgmdln)
              valt = valt - hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above WGS Geoid (conv. from another model)]
            endif
          endif
        endif
      else
!               (Convert vlat/vlon from WGS to TYO)
        if (ka == 1) then
          valt = valt - hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above Geoid [TYO] (conv. from WGS Ellips)]
        else
          if (sgmdl.ne.sgmdln) then
            call sgeoid(sgmdl)
            valt = valt + hgeoid(sngl(vlat), sngl(vlon))
            call sgeoid(sgmdln)
            valt = valt - hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above Geoid [TYO] (conv. from another model)]
          endif
        endif
        call xw84td(vlat, vlon, valt, tlat, tlon, talt)
        vlat = tlat;  vlon = tlon
!               [now vlat/vlon are in TYO (conv. from WGS)]
      endif
    else
      if (kgn == 1) then
!               (Convert vlat/vlon from TYO to WGS)
        call xtw84d(vlat, vlon, valt, wlat, wlon, walt)
        vlat = wlat;  vlon = wlon
!               [now vlat/vlon are in WGS (conv. from TYO)]
        if (kan == 1) then
          valt = valt + hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above WGS Ellips (conv. from Geoid [TYO])
        else
          if (sgmdl.ne.sgmdln) then
            call sgeoid(sgmdl)
            valt = valt + hgeoid(sngl(vlat), sngl(vlon))
            call sgeoid(sgmdln)
            valt = valt - hgeoid(sngl(vlat), sngl(vlon))
!               [now valt is above WGS Geoid (conv. from Geoid [TYO])]
          endif
        endif
      else
!               (vlat/vlon are in TYO, common IN/OUT)
!               (valt is above Geoid [TYO}, common IN/OUT)
      endif
    endif
    if (nd == 0) then                   ! start averaging
      jhsec = ihsec;
      slat = vlat;  slon = vlon;  salt = valt;  smag = vmag
      nd = 1
    else
      if (ihsec == jhsec) then          ! cont. averaging
        slat = slat + vlat;  slon = slon + vlon
        salt = salt + valt;  smag = smag + vmag
        nd = nd + 1
      else
        slat = slat / float(nd);  slon = slon / float(nd)
        salt = salt / float(nd);  smag = smag / float(nd)
        if (kr == 1) then                 ! get 1 averaged sample 
          write(2,'(1x,f11.5,a2,f11.5,a2,f8.2,a2,f8.2,a2)')   &
            (rlat+slat)/2., 'N ', (rlon+slon)/2., 'E ',   &
            (ralt+salt)/2., 'm ', (rmag+smag)/2., 'nT'
        endif
        if (ihsec == (jhsec+1)) then    ! data continuous
                                        !  averaged re-sampling active
          kr = 1
          rlat = slat;  rlon = slon;  ralt = salt;  rmag = smag
          jhsec = ihsec
          slat = vlat;  slon = vlon;  salt = valt;  smag = vmag
          nd = 1
        else                            ! data non-continuous
                                        !  averaged re-sampling once inactive
          call premsg('*** non-continuous data')
          kr = 0;  jhsec = ihsec
          slat = vlat;  slon = vlon;  salt = valt;  smag = vmag
          nd = 1
        endif
      endif
    endif
  enddo

  if (nd /= 0) then
    slat = slat / float(nd);  slon = slon / float(nd)
    salt = salt / float(nd);  smag = smag / float(nd)
    if (kr == 1) then
      write(2,'(1x,f11.5,a2,f11.5,a2,f8.2,a2,f8.2,a2)')   &
        (rlat+slat)/2., 'N ', (rlon+slon)/2., 'E ',   &
        (ralt+salt)/2., 'm ', (rmag+smag)/2., 'nT'
    endif
  endif
  close(2)
  close(1)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'XSLINA : extract StdLIN data from DPAM data'   &
                      // ' with averaged re-sampling'
    write(99,'(a,a)') ' Input AM line data filename : ', fnam1(1:lrtrim(fnam1))
    write(99,'(8x,4a)')   &
       ' 5: ', cg, '-Lat.(Deg.)    7: Alt.(m) ', ca,   &
       ' 6: ', cg, '-Lon.(Deg.)    9: Mag.(nT)', ''
    write(99,'(a,f6.2,a)') '  Averaging interval :', tint, ' sec.'
    write(99,'(a,a)') ' Output StdLIN data filename : ', fnam2(1:lrtrim(fnam2))
    write(99,'(8x,4a)') 'Geodetic in ', cgn, ';  Alt. ', can
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end