c----------( RPMEQS : reduction-to-pole anomaly by <ameqs> source )----------
      parameter (LHW=100, LXO=1201, LYO=1201)
      parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW)
      parameter (DPI=3.14159*2.)
c
      dimension s(LXM,LYM), t(LXM,LYM), f(LXO,LYO), h(LXO,LYO)
      character*80 wdr, flog, fnam1, fnam2, fnam3
      character buf*80, out*72, area*8, area1*8, area2*8
c
      call premsg('RPMEQS : reduc-to-pole anomaly by <ameqs> source')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      fnam3 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Reduction Alt. 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,*) ixl, iyl, mszx, mszy, mx, my, hnul
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      if(mx.le.0.or.my.le.0) call abendm('illegal array size')
      read(1,*) ((h(i,j),i=1,mx),j=1,my)
      close(1)
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      ixh = ixl + mszx*(mx-1)
      iyh = iyl + mszy*(my-1)
c
      call gparma('Enter AMEQS model input filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(2,file=fnam2,status='old')
    2 read(2,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(mcc1.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
      read(2,*) ixlm, iylm, mszx1, mszy1, mxm, mym, anul, dmy, klp
      if(ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy)
     *     call abendm('Header inconsistent')
      if(mxm.gt.LXM.or.mym.gt.LYM) call abendm('too large array size')
      if(mxm.lt.10.or.mym.lt.10) call abendm('illegal array size')
      read(2,*) ((s(i,j),i=1,mxm),j=1,mym)
    3 read(2,'(a)') buf
      if(buf(1:1).eq.'#') goto 3
      read(buf,'(a8,i4,4x,4i8)') area2, 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
      read(2,'(a)') buf
      read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, tnul
      if(area2.ne.area1.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header inconsistent')
      if(ixl1.ne.ixlm.or.iyl1.ne.iylm .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mxm1.ne.mxm.or.mym1.ne.mym)
     *     call abendm('Header inconsistent')
      read(2,*) ((t(i,j),i=1,mxm),j=1,mym)
      read(2,'(2x,4f8.1)') dip, dec, the, phi
      close(2)
      xlm = float(ixlm)/1000.
      ylm = float(iylm)/1000.
      ixhm = ixlm + mszx*(mxm-1)
      iyhm = iylm + mszy*(mym-1)
      if(ixl.le.ixlm.or.ixh.ge.ixhm .or.
     *   iyl.le.iylm.or.iyh.ge.iyhm) call abendm('ranges conflict')
      smx = float(mszx)
      smy = float(mszy)
      amxy = smx * smy
      vol = amxy * sqrt(amxy)
c
      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)')
     *     '--- reduc. alt. ', xl, yl,  mszx,mszy, mx, my
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- equiv.source', xlm,ylm, mszx,mszy, mxm,mym, '    var.'
      call premsg(out)
      call premsg('  Truncation of Source effect ?')
      call gparmf('    Half width (in km) : ', hw)
      hwm = hw * 1000.
      nhx = nint(hwm/smx)
      nhy = nint(hwm/smy)
c
      call gparma('Enter Reduc-to-Pole out filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      call clspin()
      open(3,file=fnam3,status='new')
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'RPMEQS : reduc-to-pole anomaly by <ameqs> source'
      write(6,'(a,a)')
     *     ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' AMEQS model infile    : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' Reduc-to-Pole outfile : ', fnam3(1:lrtrim(fnam3))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xl,  yl,  mszx, mszy, mx,  my,
     *              xlm, ylm, mszx, mszy, mxm, mym, '    var.'
      write(6,2000) hw, nhx, nhy
c
      call prompt('... calc. Reduc-Pole anom.: ')
      call dpcent(0, 0)
      do 200 j=1,my
      do 200 i=1,mx
        if(h(i,j).eq.hnul) then
          f(i,j) = 9999.9
        else
          xp = float(ixl - ixlm + mszx*(i-1))
          yp = float(iyl - iylm + mszy*(j-1))
          ib = nint(xp/smx) + 1
          jb = nint(yp/smy) + 1
          ss = 0.
          do 100 kj=-nhy,+nhy
            js = jb + kj
            y = float(kj)*smy
            do 100 ki=-nhx,+nhx
              is = ib + ki
              x = float(ki)*smx
              if(is.le.0.or.is.gt.mxm) goto 100
              if(js.le.0.or.js.gt.mym) goto 100
              if(t(is,js).eq.tnul) goto 100
              z = h(i,j) - t(is,js)
              rr = x*x + y*y + z*z
              r = sqrt(rr)
              v = z*z*2. - x*x - y*y
              ss = ss + s(is,js)*v*vol/(rr*rr*r)
  100     continue
          f(i,j) = ss
        endif
        call dpcent((j-1)*mx+i, mx*my)
  200 continue
c
      write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl, iyl, mszx, mszy, mx, my, 9999.9, 0.
      do 300 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  300 continue
      write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl, iyl, mszx, mszy, mx, my, hnul, -1.
      do 400 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  400 continue
      close(3)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        write(99,'(//a)')
     *       'RPMEQS : reduc-to-pole anomaly by <ameqs> source'
        write(99,'(a,a)')
     *       ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' AMEQS model infile    : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' Reduc-to-Pole outfile : ', fnam3(1:lrtrim(fnam3))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xl,  yl,  mszx, mszy, mx,  my,
     *                 xlm, ylm, mszx, mszy, mxm, mym, '    var.'
        write(99,2000) hw, nhx, nhy
        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, 'reduc. surf.', 2f10.2, 2i8, 2i6/
     *       3x, 'equiv.source', 2f10.2, 2i8, 2i6, a8)
 2000 format(' Half width of calc.window :', f6.2,
     *       ' km (', i4, ' *', i4, '  source-mesh )')
      end
