!--------( GTRF : translate RF-residual between models )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), g(:,:)
  character(80) :: wdr, flog, fnam1, fnam3, buf
  character(72) :: out;  character(20) :: sdtm;  character(10) :: anul
  character(12) :: srf0, srf1
  character(8) :: area, area1, narea, salt
  integer :: ldr, i, j, kn, mrf0, mdl0, mrf1, mdl1
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  real :: vnul, alt, xs, ys, fi, fk, vnul1, year, elv, rd
  real :: x, y

  call premsg('GTRF : translate RF-residual between models')
  call premsg('        by 1) Adding OLD reference model field, and')
  call premsg('           2) Subtracting NEW reference model field')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam3 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input source data filename ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  read(1,'(a)') buf
  do while (buf(1:1) == '#')
    read(1,'(a)') buf
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
    if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
    nc = ncc
    if (ncc < 200) then
      nc = ncc; ncc = nc + 800
    else
      ncc = ncc - 200; nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  endif
  if ((nc >= 1) .and. (nc <= 62)) then
    iorg = 0; ispa = 0; ispb = 0
    korg = (nc-30)*360 - 180
    if (nc > 60) korg = 0
    if (nc == 61) iorg = +5400
    if (nc == 62) iorg = -5400
  endif
  read(1,'(a)') buf
  read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
  write(salt,'(1x,f7.0)') alt
  if (alt == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  xs = float(ixs)/1000.
  ys = float(iys)/1000.

  allocate(f(mx,my))

  write(anul,'(f9.2)') vnul
  if ((anul(9:9) /= '0') .or. (anul(1:1) /= ' ')) then
    call abendm('unable to handle vnul value')
  endif
  write(out,'(a12,a8,6x,a,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6,a8)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '---   source    ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  read(1,*) ((f(i,j),i=1,mx),j=1,my)

  if (alt > 0.) then
    kn = 0; elv = alt
  else if (alt < 0.) then
    kn = 0; elv = 0.
  else
    kn = 1
    read(1,'(a)') buf
    do while (buf(1:1) == '#')
      read(1,'(a)') buf
    enddo
    if (buf(13:16) == '    ') then        !-- v2005
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
      nc1 = ncc1
      if (ncc1 < 200) then
        nc1 = ncc1; ncc1 = nc1 + 800
      else
        ncc1 = ncc1 - 200; nc1 = ncc1
      endif
    else                                  !-- v2018
      read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
    endif
    if ((nc1 >= 1) .and. (nc1 <= 62)) then
      iorg1 = 0;  ispa1 = 0;  ispb1 = 0
      korg1 = (nc1-30)*360 - 180
      if (nc1 > 60) korg1 = 0
      if (nc1 == 61) iorg1 = +5400
      if (nc1 == 62) iorg1 = -5400
    endif
    if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg) .or.   &
        (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
      call abendm('Header1 inconsistent')
    read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
    if ((ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx) .or.   &
        (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
      call abendm('Header2 inconsistent')
    write(anul,'(f10.3)') vnul1
    if ((anul(9:10) /= '00') .or. (anul(1:1) /= ' '))   &
       call abendm('unable to handle vnul1 value')
    allocate(g(mx,my))
    read(1,*) ((g(i,j),i=1,mx),j=1,my)
  endif
  close(1)

  call gparmf('Enter Survey-Year  ==> ', year)
  call premsg('Select OLD RF model to be added')
  call gparmi('  1:IGRF(tentative) 2:DGRF 3:PGFR or 0:byGen. ==> ', mrf0)
  if ((mrf0 < 0) .or. (mrf0 > 3)) call abendm('invalid value')
  srf0 = 'IGRF-       '
  if (mrf0 == 0) then
    call premsg('  Enter Generation Number [1-]  or  0')
    call premsg('   (0: [Input data must be TotalForces, not Residuals]')
    call gparmi('       No RF model will be Added)      ==> ', mdl0)
    if ((mdl0 < 0) .or. (mdl0 > 99)) call abendm('invalid value')
    write(srf0(6:7),'(i2)') mdl0
    if (mdl0 == 0) srf0 = '0  (InputTF)'
  else
    call gparmi('  RF model base year  ==> ', mdl0)
    if (mdl0 < 1900) call abendm('invalid value')
    write(srf0(6:9),'(i4)') mdl0
    if (mrf0 == 2) srf0(1:1) = 'D'
    if (mrf0 == 3) srf0(1:1) = 'P'
  endif
  call premsg('Select NEW RF model to be subtracted')
  call gparmi('  1:IGRF(tentative) 2:DGRF 3:PGFR or 0:byGen. ==> ', mrf1)
  if ((mrf1 < 0) .or. (mrf1 > 3)) call abendm('invalid value')
  srf1 = 'IGRF-       '
  if (mrf1 == 0) then
    call premsg('  Enter Generation Number [1-]  or  0')
    call premsg('   (0: [Output data is TotalForces]')
    call gparmi('       No RF model will be Subtracted) ==> ', mdl1)
    if ((mdl1 < 0) .or. (mdl1 > 99)) call abendm('invalid value')
    write(srf1(6:7),'(i2)') mdl1
    if (mdl1 == 0) srf1 = '0 (OutputTF)'
  else
    call gparmi('  RF model base year  ==> ', mdl1)
    if (mdl1 < 1900) call abendm('invalid value')
    write(srf1(6:9),'(i4)') mdl1
    if (mrf1 == 2) srf0(1:1) = 'D'
    if (mrf1 == 3) srf0(1:1) = 'P'
  endif
  call gparma('Enter Name-label for new file  ==> ', 8, narea)
  call gparma('Enter output new data filename ==> ', 81-ldr, fnam3(ldr:80))
  call clspin()
  open(3,file=fnam3,status='new')
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'GTRF : translate RF-residual between models'
  write(6,'(a)') '        by 1) Adding OLD reference model field, and'
  write(6,'(a)') '           2) Subtracting NEW reference model field'
  write(6,'(a,a)') ' Source infile  :   ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') '    areaname :        ', area
  write(6,'(a,f8.2)') '    year of survey  : ', year
  write(6,'(a,a8)')   '    flying alt. (m) : ', salt
  write(6,'(a,a)') '    Reference Field : ', srf0
  write(6,'(a,a)') ' Output new file :  ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') '    areaname :        ', narea
  write(6,'(a,a)') '    Reference Field : ', srf1
  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', narea, 'Coord.', ncc, iorg, korg, ispa, ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'grid out=in ', xs, ys, mszx, mszy, mx, my, salt
  write(3,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
  write(3,'(a2,a8,4x,i4,4i8)') '# ', narea, ncc, iorg, korg, ispa, ispb
  write(3,'(a,f8.2)') '#    year of survey  : ', year
  write(3,'(a,a8)')   '#    flying alt. (m) : ', salt
  write(3,'(a,a)')    '#    Reference Field : ', srf0
  write(3,'(a)')      '# converted into another RF residual'
  write(3,'(a,a)')    '#    Reference Field : ', srf1
  write(3,'(a8,4x,i4,4i8)') narea, ncc, iorg,korg, ispa,ispb
  write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy, mx,my, vnul, alt

  call cvinit(ncc, float(iorg), float(korg), float(ispa), float(ispb))
  call cvenik(ys, xs, fi, fk)
  if (mdl1 > 0) then
    if (mrf0 == 0) then
      call gigrf(mdl0, year)
    else if (mrf0 == 1) then
      call sigrf(float(mdl0))
      call tyear(year)
    else if (mrf0 == 2) then
      call sdgrf(float(mdl0))
      call tyear(year)
    else if (mrf0 == 3) then
      call spgrf(float(mdl0))
      call tyear(year)
    endif
    call dpcini('... adding Old RF model : ')
    elv = alt
    do j=1,my
      y = ys + (j-1)*mszy/1000.
      do i=1,mx
        x = xs + (i-1)*mszx/1000.
        if (f(i,j) /= vnul) then
          call cvenik(y, x, fi, fk)
          if (kn == 1) elv = g(i,j)
          call igrfc(fi/60., fk/60., elv, rd)
          f(i,j) = f(i,j) + rd
        endif
        call dpcent((j-1)*mx+i, mx*my)
      enddo
    enddo
  endif
  if (mdl1 > 0) then
    if (mrf1 == 0) then
      call gigrf(mdl1, year)
    else if (mrf1 == 1) then
      call sigrf(float(mdl1))
      call tyear(year)
    else if (mrf1 == 2) then
      call sdgrf(float(mdl1))
      call tyear(year)
    else if (mrf1 == 3) then
      call spgrf(float(mdl1))
      call tyear(year)
    endif
    call dpcini('... subtracting New RF  : ')
    do j=1,my
      y = ys + (j-1)*mszy/1000.
      do i=1,mx
        x = xs + (i-1)*mszx/1000.
        if (f(i,j) /= vnul) then
          call cvenik(y, x, fi, fk)
          if (kn == 1) elv = g(i,j)
          call igrfc(fi/60., fk/60., elv, rd)
          f(i,j) = f(i,j) - rd
        endif
        call dpcent((j-1)*mx+i, mx*my)
      enddo
    enddo
  endif
  do j=1,my
    write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo

  if (kn == 1) then
    write(3,'(a)') '# Var.Obs.Altitude data are copied from Source'
    write(3,'(a8,4x,i4,4i8)') narea, ncc, iorg, korg, ispa, ispb
    write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy, mx,my, vnul1,-1.
    do j=1,my
      write(3,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'Var.Obs.Altitude data are copied from Source'
  endif
  close(3)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'GTRF : translate RF-residual between models'
    write(99,'(a)') '        by 1) Adding OLD reference model field, and'
    write(99,'(a)') '           2) Subtracting NEW reference model field'
    write(99,'(a,a)') ' Source infile  :  ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') '    areaname :        ', area
    write(99,'(a,f8.2)') '    year of survey  : ', year
    write(99,'(a,a8)')   '    flying alt. (m) : ', salt
    write(99,'(a,a)') '    Reference Field : ', srf0
    write(99,'(a,a)') ' Output new file : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') '    areaname :        ', narea
    write(99,'(a,a)') '    Reference Field : ', srf1
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', narea, 'Coord.', ncc, iorg, korg, ispa, ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'grid out=in ', xs, ys, mszx, mszy, mx, my, salt
    if (kn == 1)   &
      write(99,'(a)') 'Var.Obs.Altitude data are copied from Source'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end