!**************************************************
!   FCOMP : Compensate Survey Flight Data         *
!             with removing linear Trend surface  *
!**************************************************
!           must be linked with sml5.o
!--------------------------------------------------
  implicit none
  integer :: lwkdir, lrtrim
  real :: coef(9), ctrd(3)
  character(80) :: wdr, flog, fcmp, fnam, fout
  character(150) :: buf; character(20) :: sdtm
  integer :: ldr, iymd, km, IER
  real :: sn, fid, falt, tres, fx, fy, fz, utc, trend, corr, tcor, rand
  real(8) :: wtm, tmag, dido, dkdo, dsi, dsk

  call premsg('FCOMP : Compensate Survey Flight Data')
  call premsg('         with removing linear Trend surface')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fcmp = wdr;  fnam = wdr;  fout = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter CompBox data filename ==> ', 81-ldr, fcmp(ldr:80))
  open(1,file=fcmp,status='old')
  call gparma('Enter Input data filename   ==> ', 81-ldr, fnam(ldr:80))
  open(2,file=fnam,status='old')
  call gparma('Enter Output data filename  ==> ', 81-ldr, fout(ldr:80))
  call clspin()
  open(3,file=fout,status='new')
  dsi = 0.;  dsk = 0.;  sn = 0.
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) cycle
    read(buf,*) fid,iymd,wtm,km,dido,dkdo,falt,tmag,tres,fx,fy,fz,utc
    if (km == (km/2*2)) call abendm('Input data already compensated')
    dsi = dsi + dido;  dsk = dsk + dkdo;  sn = sn + 1.
  enddo
  rewind(1)
  dsi = dsi / sn;  dsk = dsk / sn
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'FCOMP : Compensate Survey Flight Data'
  write(6,'(a)') '         with removing linear Trend surface'
  write(6,'(a,a)') ' CompBox filename : ', fcmp(1:lrtrim(fcmp))
  write(6,'(a,a)') ' Input filename   : ', fnam(1:lrtrim(fnam))
  write(6,'(a,a)') ' Output filename  : ', fout(1:lrtrim(fout))
  write(6,'(a)')   ' ---'
  write(6,*) '  Average Lat./Long. : ', sngl(dsi), sngl(dsk)

  call sm5opn()
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) cycle
    read(buf,*) fid,iymd,wtm,km,dido,dkdo,falt,tmag,tres,fx,fy,fz,utc
    call sm5ex(sngl((dido-dsi)*60.), sngl((dkdo-dsk)*60.), fx, fy, fz, tres)
  enddo
  close(1)
  call sm5cls(ctrd, 3, coef, 9)
  write(6,*) ' Coef. of Trend Surface : ', ctrd
  write(6,*) ' Compensation coefficients : '
  write(6,*) '  (x)   ', coef(1), coef(2), coef(3)
  write(6,*) '  (xx)  ', coef(4), coef(7), coef(9)
  write(6,*) '  (xy)  ', coef(5), coef(8), coef(6)

  call premsg('Processing for Survey data ...')
  do
    read(2,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) then
      write(3,'(a)') buf(1:lrtrim(buf))
    else
      read(buf,*) fid,iymd,wtm,km,dido,dkdo,falt,tmag,tres,fx,fy,fz,utc
      call sm5rv(0., 0., fx, fy, fz, corr)
      call sm5rv(sngl((dido-dsi)*60.), sngl((dkdo-dsk)*60.), 0.,0.,0., trend)
      corr = corr - ctrd(1);  tmag = tmag - corr
      tcor = tres - corr;  rand = tcor - trend
      write(3,'(i8,1x,i8,1x,f9.2,1x,i2,1x,f11.7,1x,f12.7,1x,f7.2,$)')   &
        nint(fid),iymd,wtm,(km-1),dido,dkdo,falt
      write(3,'(2(1x,f8.2), 3(1x,f7.3), 1x,f9.2, 4(1x,f8.2))')   &
        tmag,tcor,fx,fy,fz,utc,tres,corr,rand,trend
    endif
  enddo
  close(2)
  close(3)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'FCOMP : Compensate Survey Flight Data'
    write(99,'(a)') '         with removing linear Trend surface'
    write(99,'(a,a)') ' CompBox filename : ', fcmp(1:lrtrim(fcmp))
    write(99,'(a,a)') ' Input filename   : ', fnam(1:lrtrim(fnam))
    write(99,'(a,a)') ' Output filename  : ', fout(1:lrtrim(fout))
    write(99,'(a)')   ' ---'
    write(99,*) '  Average Lat./Long. : ', sngl(dsi), sngl(dsk)
    write(99,*) '  Coef. of Trend Surface : ', ctrd
    write(99,*) '  Compensation coefficients : '
    write(99,*) '   (x)   ', coef(1), coef(2), coef(3)
    write(99,*) '   (xx)  ', coef(4), coef(7), coef(9)
    write(99,*) '   (xy)  ', coef(5), coef(8), coef(6)
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end