#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
c----------( ADEQC : <adeq> inversion until converge )----------
      parameter (LHW=100, LXO=1201, LYO=1201, MXPT=500000)
      parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW, LFW=LHW+LHW+1)
      parameter (LLFW=LFW*LFW)
c
      dimension g(LLFW)
      dimension f(MXPT), r(MXPT), q(MXPT)
      dimension s(LXM,LYM), p(LXM,LYM), v(LXM,LYM), t(LXM,LYM)
      character*80 wdr, flog, fnam1, fnam2, fnam3
      character buf*80, out*72, area*8, area1*8
c
      call premsg('ADEQC : <adeq> inversion until converge')
      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 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,1x,f8.1)') x, y, h, f(ndpt)
      goto 1
    2 close(1)
c
      call gparma('Enter CFUP matrix data filename  ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(20,file=fnam2,form='unformatted',status='old')
      read(20) area, ncc, iorg, korg, ispa, ispb,
     *         ixl, iyl, mszx, mszy, mx, my,
     *         nhx, nhy, hw, vd
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      ixlm = ixl - mszx*nhx
      iylm = iyl - mszy*nhy
      xlm = float(ixlm)/1000.
      ylm = float(iylm)/1000.
      mxm = mx + nhx + nhx
      mym = my + nhy + nhy
      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')
      if(nhx.gt.LHW.or.nhy.gt.LHW)
     *     call abendm('shortage of work area within truncation')
      if(nhx.lt.10.or.nhy.lt.10) call abendm('truncation too short')
      nfx = nhx + 1 + nhx
      nfy = nhy + 1 + nhy
      nnfw = nfx * nfy
      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. surf.', xl, yl,  mszx,mszy, mx, my
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- eq.src.anom.', xlm,ylm, mszx,mszy, mxm,mym, '    var.'
      call premsg(out)
c
      call gparma('Enter ADEQ model in/out filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(21,file=fnam3,status='old')
    3 read(21,'(a)') 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
      read(21,'(a)') buf
      read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, anul, dmy, klp
      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('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(21,*) ((s(i,j),i=1,mxm),j=1,mym)
    4 read(21,'(a)') buf
      if(buf(1:1).eq.'#') goto 4
      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
      read(21,'(a)') buf
      read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, tnul
      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('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(21,*) ((t(i,j),i=1,mxm),j=1,mym)
      rewind(21)
c
      call premsg('Threshold percentage to regard converged')
      call gparmf('   (cont. 5 loops) [default: 2%] ==> ', thrsh)
      if (thrsh.lt.0. .or. thrsh.ge.20.)
     *     call abendm('invalid parameter')
      if (thrsh.eq.0.) thrsh = 2.
      call clspin()
c
      call prompt('.... read CFUP matrix (0) : ')
      call dpcent(0, 0)
      kc = 0
    5 read(20,end=6) k, ic, jc, (g(i),i=1,nnfw)
      kc = kc + 1
      if(k.ne.kc) call abendm('CFUP file broken')
      call dpcent(kc, ndpt)
      goto 5
    6 rewind(20)
      if(kc.ne.ndpt) call abendm('CFUP file broken')
c
      write(6,'(a)')   '----------------------------------------'
      write(6,'(a)')   'ADEQC : <adeq> inversion until converge'
      write(6,'(a,a)')
     *     ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' CFUP matrix infile    : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' ADEQ model i/o file   : ', 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, vd
c
        ss = 0.
        do 200 j=1,mym
        do 200 i=1,mxm
          if(t(i,j).ne.tnul) ss = ss + s(i,j)*s(i,j)
          p(i,j) = 0.
  200   continue
        ds = sqrt(ss / float(mxm*mym))
c
      call conv1(g,nfx,nfy, LXM,s,mxm,mym, r,ndpt)
      ss = 0.
      do 10 k=1,ndpt
        r(k) = f(k) - r(k)
        ss = ss + r(k)*r(k)
   10 continue
      dv = sqrt(ss / float(ndpt))
      lp = 0
      write(6,3000) klp, lp, ds, dv
      open(9,status='scratch')
      write(9,*) klp, lp, ds, dv
c
      pvv = 1.
    7   klp = klp + 1
        rv = dv
        pvvp = pvv
        call conv2(g,nfx,nfy, r,ndpt, LXM,v,mxm,mym)
        pvv = 0.
        do 20 j=1,mym
        do 20 i=1,mxm
          pvv = pvv + v(i,j)*v(i,j)
   20   continue
        a = pvv / pvvp
        do 30 j=1,mym
        do 30 i=1,mxm
          p(i,j) = v(i,j) + a*p(i,j)
   30   continue
        call conv1(g,nfx,nfy, LXM,p,mxm,mym, q,ndpt)
        pqq = 0.
        do 40 k=1,ndpt
          pqq = pqq + q(k)*q(k)
   40   continue
        a = pvv / pqq
        ss = 0.
        do 50 j=1,mym
        do 50 i=1,mxm
          s(i,j) = s(i,j) + a*p(i,j)
          if(t(i,j).ne.tnul) ss = ss + p(i,j)*p(i,j)
   50   continue
        ds = a * sqrt(ss / float(mxm*mym))
        ss = 0.
        do 60 k=1,ndpt
          r(k) = r(k) - a*q(k)
          ss = ss + r(k)*r(k)
   60   continue
        dv = sqrt(ss / float(ndpt))
        if (((rv-dv)/rv)*100. .ge. thrsh) then
          lp = 0
        else
          lp = lp + 1
        endif
        write(6,3000) klp, lp, ds, dv
        write(9,*) klp, lp, ds, dv
        if (lp.lt.5.and.dv.ge.0.1) goto 7
c
      close(20)
      write(6,'(1x)')
      endfile(9)
      rewind(9)
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., klp
      do 410 j=1,mym
        do 400 i=1,mxm
          if(s(i,j).lt.-9999.8) s(i,j) = -9999.8
          if(s(i,j).gt.+9999.8) s(i,j) = +9999.8
  400   continue
        write(21,'((f7.1,9(1x,f7.1)))') (s(i,j),i=1,mxm)
  410 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 420 j=1,mym
        write(21,'((f7.1,9(1x,f7.1)))') (t(i,j),i=1,mxm)
  420 continue
      close(21)
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)') 'ADEQC : <adeq> inversion until converge'
        write(99,'(a,a)')
     *       ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' CFUP matrix infile    : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' ADEQ model i/o file   : ', 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, vd
   92   read(9,*,end=93) klp, lp, ds, dv
        write(99,3000) klp, lp, ds, dv
        goto 92
   93   close(99)
      endif
      close(9)
      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.')
 3000 format(' loop', i4, i3, '   rms.Step :', f8.2,
     *       '   rms.Residual :', f7.2)
      end
c---
      subroutine conv1(g,nfx,nfy, lsx,s,mxm,mym, r,ndpt)
      dimension g(nfx,nfy)
      dimension s(lsx,mym), r(ndpt)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call prompt('.... read CFUP matrix (1) : ')
      call dpcent(0, 0)
      read(20) area, (i,k=1,15)
      do 1 np=1,ndpt
        read(20) k, ic, jc, g
        rp = 0.
        kis = ic - nhx
        kjs = jc - nhy
        do 2 j=1,nfy
          js = kjs + j
          if(js.le.0.or.js.gt.mym) goto 2
          do 3 i=1,nfx
            is = kis + i
            if(is.le.0.or.is.gt.mxm) goto 3
            rp = rp + g(i,j)*s(is,js)
    3     continue
    2   continue
        r(np) = rp
        call dpcent(np, ndpt)
    1 continue
      rewind(20)
      return
      end
c---
      subroutine conv2(g,nfx,nfy, r,ndpt, lsx,s,mxm,mym)
      dimension g(nfx,nfy)
      dimension s(lsx,mym), r(ndpt)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call prompt('.... read CFUP matrix (2) : ')
      call dpcent(0, 0)
      read(20) area, (i,k=1,15)
      do 1 js=1,mym
      do 1 is=1,mxm
        s(is,js) = 0.
    1 continue
      do 2 np=1,ndpt
        read(20) k, ic, jc, g
        kis = ic - nhx
        kjs = jc - nhy
        ss = 0.
        do 3 j=1,nfy
          js = kjs + j
          if(js.le.0.or.js.gt.mym) goto 3
          do 4 i=1,nfx
            is = kis + i
            if(is.le.0.or.is.gt.mxm) goto 4
            s(is,js) = s(is,js) + g(i,j)*r(np)
    4     continue
    3   continue
        call dpcent(np, ndpt)
    2 continue
      rewind(20)
      return
      end
