c--------( EDEQ : preparation for <adeq> inversion )--------
      parameter (LHW=100, LXO=1201, LYO=1201, MXPT=500000)
      parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW)
      parameter (DPI=3.141593*2.)
c
      dimension f(3,MXPT), t(LXO,LYO), tt(LXM,LYM)
      dimension g(-LHW:+LHW, -LHW:+LHW)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character buf*80, out*72, area*8, sdtm*20
c
      call premsg('EDEQ : preparation for <adeq> inversion')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      fnam3 = wdr
      fnam4 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter StdLIN data filename     ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(1,file=fnam1,status='old')
      ndpt = 0
    1 read(1,'(a)',end=2) buf
      if(buf(1:1).eq.'&'.or.buf(1:1).eq.'%'.or.buf(1:1).eq.'#') goto 1
      ndpt = ndpt + 1
      if(ndpt.gt.MXPT) call abendm('too many data points')
      read(buf,'(f10.4,1x,f11.4,1x,f8.1)') (f(i,ndpt),i=1,3)
      goto 1
    2 close(1)
c
      call gparma('Enter Reduction Alt. filename  ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(2,file=fnam2,status='old')
    3 read(2,'(a)') buf
      if(buf(1:1).eq.'#') then
        call premsg('  ' // buf(1:lrtrim(buf)))
        goto 3
      endif
      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(2,'(a)') buf
      read(buf,*) ixl, iyl, mszx, mszy, mx, my, tnul
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      read(2,*) ((t(i,j),i=1,mx),j=1,my)
      close(2)
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)')
     *     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6)')
     *     '--- reduc. alt. ', xl, yl, mszx, mszy, mx, my
      call premsg(out)
      smx = float(mszx)
      smy = float(mszy)
      amxy = smx * smy
      call premsg('')
      call premsg('Specify Equivalent Source Parameter')
      call premsg('  Vertical distance from reduction surface ?')
      call gparmf('    Downward (in m)    : ', vd)
      if(vd.le.0.) call abendm('invalid value')
      call premsg('  Truncation of Source effect ?')
      call gparmf('    Half width (in km) : ', hw)
      hwm = hw * 1000.
      nhx = nint(hwm/smx)
      nhy = nint(hwm/smy)
      if(nhx.gt.LHW.or.nhy.gt.LHW) call abendm('short of work area')
      if(nhx.lt.10.or.nhy.lt.10) call abendm('truncation too short')
      nfx = nhx + nhx
      nfy = nhy + nhy
      mxm = mx + nfx
      mym = my + nfy
      ixlm = ixl - mszx*nhx
      iylm = iyl - mszy*nhy
      xlm = float(ixlm)/1000.
      ylm = float(iylm)/1000.
c
c      frnfxy = float((nhx+1+nhx)*(nhy+1+nhy) + 4) * 4. / 1024.
c      frndpt = float(ndpt + 1) / 1024. / 1024.
c      if((frnfxy*frndpt).gt.2.) then
c        call abendm('CFUP file-size exceeds 2GB.')
c      endif
c
      call gparma('Enter CFUP matrix output filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(20,file=fnam3,form='unformatted',status='new')
      call gparma('Enter ADEQ model output filename  ==> ',
     *            81-ldr, fnam4(ldr:80))
      call clspin()
      call strdtm(sdtm)
      open(21,file=fnam4,status='new')
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'EDEQ : preparation for <adeq> inversion'
      write(6,'(a,a)')
     *     ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' CFUP matrix outfile   : ', fnam3(1:lrtrim(fnam3))
      write(6,'(a,a)')
     *     ' ADEQ model outfile    : ', fnam4(1:lrtrim(fnam4))
      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, vd
c
      call dpcini('... gen. Source alt. data : ')
      do 100 jp=1,mym
      do 100 ip=1,mxm
        tt(ip,jp) = tnul
  100 continue
      do 110 j=1,my
      do 110 i=1,mx
        if(t(i,j).ne.tnul) tt(i+nhx,j+nhy) = t(i,j) - vd
  110 continue
      hhwm4 = hwm*hwm * 4.
      do 160 jp=1,mym
      do 160 ip=1,mxm
        if(tt(ip,jp).ne.tnul) goto 4
        ss = hhwm4
        do 130 js=-nhy,+nhy
          j = jp - nhy + js
          if(j.le.0.or.j.gt.my) goto 130
          y = float(js)*smy
          do 120 is=-nhx,+nhx
            i = ip - nhx + is
            if(i.le.0.or.i.gt.mx) goto 120
            if(t(i,j).eq.tnul) goto 120
            x = float(is)*smx
            rr = x*x + y*y
            if(rr.ge.ss) goto 120
            ss = rr
  120     continue
  130   continue
        if(ss.eq.hhwm4) goto 4
        ss = ss + ss
        w = 0.
        s = 0.
        do 150 js=-nfy,+nfy
          j = jp - nhy + js
          if(j.le.0.or.j.gt.my) goto 150
          y = float(js)*smy
          do 140 is=-nfx,+nfx
            i = ip - nhx + is
            if(i.le.0.or.i.gt.mx) goto 140
            if(t(i,j).eq.tnul) goto 140
            x = float(is)*smx
            rr = x*x + y*y
            if(rr.ge.ss) goto 140
            s = s + t(i,j)*(ss-rr)
            w = w + (ss-rr)
  140     continue
  150   continue
        tt(ip,jp) = s/w - vd
    4   call dpcent((jp-1)*mxm+ip, mxm*mym)
  160 continue
c
      call dpcini('... gen. CFUP matrix data : ')
      call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
      write(20) area, ncc, iorg,korg, ispa,ispb,
     *          ixl, iyl, mszx, mszy, mx, my,
     *          nhx, nhy, hw, vd
      do 230 k=1,ndpt
        call cviken(f(1,k), f(2,k), yp, xp)
        xp = (xp-xlm)*1000.
        yp = (yp-ylm)*1000.
        hp = f(3,k)
        ic = nint(xp/smx)
        jc = nint(yp/smy)
        do 200 kj=-LHW,+LHW
        do 200 ki=-LHW,+LHW
          g(ki,kj) = 0.
  200   continue
        do 220 kj=-nhy,+nhy
          js = jc + kj
          if(js.lt.0.or.js.ge.mym) goto 220
          y = yp - float(js)*smy
          do 210 ki=-nhx,+nhx
            is = ic + ki
            if(is.lt.0.or.is.ge.mxm) goto 210
            ht = tt(is+1,js+1)
            if (ht.ne.tnul) then
              x = xp - float(is)*smx
              z = hp - ht
              rr = x*x + y*y + z*z
              r = sqrt(rr)
              g(ki,kj) = z*amxy / (rr*r*DPI)
            endif
  210     continue
  220   continue
        write(20) k, ic, jc, ((g(ki,kj),ki=-nhx,+nhx),kj=-nhy,+nhy)
        call dpcent(k, ndpt)
  230 continue
      close(20)
c
      write(21,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)')
     *     ixlm, iylm, mszx, mszy, mxm, mym, 9999.9, 0., 0
      do 40 j=1,mym
        write(21,'((f7.1,9(1x,f7.1)))') (0.,i=1,mxm)
   40 continue
      write(21,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(21,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixlm, iylm, mszx, mszy, mxm, mym, tnul, -1.
      do 50 j=1,mym
        write(21,'((f7.1,9(1x,f7.1)))') (tt(i,j),i=1,mxm)
   50 continue
      close(21)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)') 'EDEQ : preparation for <adeq> inversion'
        write(99,'(a,a)')
     *       ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' CFUP matrix outfile   : ', fnam3(1:lrtrim(fnam3))
        write(99,'(a,a)')
     *       ' ADEQ model outfile    : ', fnam4(1:lrtrim(fnam4))
        write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
        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, vd
        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, 'reduc. surf.', 2f10.2, 2i8, 2i6/
     *       3x, 'eq.src.anom.', 2f10.2, 2i8, 2i6, a8)
 2000 format(' Half width of calc.window :', f6.2, ' km (', i4,
     *       ' *', i4, '  mesh )'/ ' Equiv. source anom. :',
     *       f8.1, ' m below Reduction alt.')
      end
