c--------( GTRF : translate RF-residual between models )--------
      parameter (LX=4001, LY=4001)
c
      dimension f(LX,LY)
      character*80 wdr, flog, fnam1, fnam3
      character area*8, narea*8, anul*10, area1*8, salt*8
      character srf0*9, srf1*9, buf*80, out*72, sdtm*20
c
      call premsg('GTRF : translate RF-residual between models')
      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')
    1 read(1,'(a)') buf
      if(buf(1:1).eq.'#') goto 1
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.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.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      if(mx.gt.LX.or.my.gt.LY) call abendm('too large array size')
      write(anul,'(f9.2)') vnul
      if(anul(9:9).ne.'0'.or.anul(1:1).ne.' ')
     *     call abendm('unable to handle vnul value')
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      write(out,'(a,2a10,2x,2a8,2a6,a7)')
     *     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *     'alt'
      call premsg(out)
      write(out,'(a,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)
      call gparmf('Enter Survey-Year  ==> ', year)
      call premsg('Select OLD RF model Spec. or by Generation')
      call gparmi('  1:IGRF(tentative) 2:DGRF 3:PGFR or 0:byGen. ==> ',
     *            mrf0)
      if (mrf0.lt.0.or.mrf0.gt.3) call abendm('invalid value')
      srf0 = 'IGRF     '
      if (mrf0.eq.0) then
        call gparmi('  Generation Number ==> ', mdl0)
        if (mdl0.lt.1.or.mdl0.gt.99) call abendm('invalid value')
        srf0(5:5) = '-'
        write(srf0(6:7),'(i2)') mdl0
      else
        call gparmi('  RF base year  ==> ', mdl0)
        if (mdl0.lt.1900) call abendm('invalid value')
        write(srf0(6:9),'(i4)') mdl0
        if (mrf0.eq.2) srf0(1:1) = 'D'
        if (mrf0.eq.3) srf0(1:1) = 'P'
      endif
      call premsg('Select NEW RF model Spec. or by Generation')
      call gparmi('  1:IGRF(tentative) 2:DGRF 3:PGFR or 0:byGen. ==> ',
     *            mrf1)
      if (mrf1.lt.0.or.mrf1.gt.3) call abendm('invalid value')
      srf1 = 'IGRF     '
      if (mrf1.eq.0) then
        call gparmi('  Generation Number ==> ', mdl1)
        if (mdl1.lt.1.or.mdl1.gt.99) call abendm('invalid value')
        srf1(5:5) = '-'
        write(srf1(6:7),'(i2)') mdl1
      else
        call gparmi('  RF base year  ==> ', mdl1)
        if (mdl1.lt.1900) call abendm('invalid value')
        write(srf1(6:9),'(i4)') mdl1
        if (mrf1.eq.2) srf0(1:1) = 'D'
        if (mrf1.eq.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,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,1000) narea, ncc, iorg,korg, ispa,ispb,
     *              xs, ys, mszx, mszy, mx, my, salt
c
      write(3,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
      write(3,'(a2,a8,i4,4x,4i8)')
     *     '# ', area, 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,i4,4x,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
c
      call cvinit(ncc, 0.,0.,0.,0.)
      call cvenik(ys, xs, fi, fk)
      if (mrf0.eq.0) then
        call gigrf(mdl0, year)
      else if (mrf0.eq.1) then
        call sigrf(float(mdl0))
        call tyear(year)
      else if (mrf0.eq.2) then
        call sdgrf(float(mdl0))
        call tyear(year)
      else if (mrf0.eq.3) then
        call spgrf(float(mdl0))
        call tyear(year)
      endif
      call igrfc(fi/60., fk/60., alt, rdb)
      call dpcini('... adding old model value: ')
      do 100 j=1,my
        y = ys + (j-1)*mszy/1000.
        do 100 i=1,mx
          x = xs + (i-1)*mszx/1000.
          if(f(i,j).ne.vnul) then
            call cvenik(y, x, fi, fk)
            call igrfc(fi/60., fk/60., alt, rd)
            f(i,j) = f(i,j) + (rd-rdb)
          endif
          call dpcent((j-1)*mx+i, mx*my)
  100 continue
      if (mrf1.eq.0) then
        call gigrf(mdl1, year)
      else if (mrf1.eq.1) then
        call sigrf(float(mdl1))
        call tyear(year)
      else if (mrf1.eq.2) then
        call sdgrf(float(mdl1))
        call tyear(year)
      else if (mrf1.eq.3) then
        call spgrf(float(mdl1))
        call tyear(year)
      endif
      call dpcini('... subtracting new model : ')
      do 200 j=1,my
        y = ys + (j-1)*mszy/1000.
        do 200 i=1,mx
          x = xs + (i-1)*mszx/1000.
          if(f(i,j).ne.vnul) then
            call cvenik(y, x, fi, fk)
            call igrfc(fi/60., fk/60., alt, rd)
            f(i,j) = f(i,j) - (rd-rdb)
          endif
          call dpcent((j-1)*mx+i, mx*my)
  200 continue
      do 300 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  300 continue
      kn = 0
c
    3 read(1,'(a)',end=8) buf
      if(buf(1:1).eq.'#') goto 3
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(ncc1.ge.200) nc1 = ncc1 - 200
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if(area1.ne.area.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header1 inconsistent')
      read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
      if(ixs1.ne.ixs.or.iys1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my)
     *     call abendm('Header2 inconsistent')
      write(anul,'(f10.3)') vnul1
      if(anul(9:10).ne.'00'.or.anul(1:1).ne.' ')
     *     call abendm('unable to handle vnul1 value')
      read(1,*) ((f(i,j),i=1,mx),j=1,my)
      close(1)
      write(3,'(a)') '# Var.Obs.Altitude data are copied from Source'
      write(3,'(a8,i4,4x,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 400 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  400 continue
      kn = 1
      write(6,'(a)') 'Var.Obs.Altitude data are copied from Source'
c
    8 close(1)
      close(3)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)') 'GTRF : translate RF-residual between models'
        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,1000) narea, ncc, iorg,korg, ispa,ispb,
     *                 xs, ys, mszx, mszy, mx, my, salt
        if (kn.eq.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
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'grid out=in ', 2f10.2, 2i8, 2i6, a8)
      end
