c--------( AMEQC : <ameq> inversion until converge )--------
      parameter (LHW=100, LXO=1201, LYO=1201)
      parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW, LFW=LHW+LHW+1)
      parameter (LLFW=LFW*LFW)
c
      dimension g(LLFW)
      dimension f(LXO,LYO), r(LXO,LYO), q(LXO,LYO)
      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, salt*8, sdtm*20
c
      call premsg('AMEQC : <ameq> 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 AM mesh 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,*) ixsf, iysf, mszxf, mszyf, mxf, myf, vnul, alt
      if(mxf.gt.LXO.or.myf.gt.LYO) call abendm('too large array size')
      if(mxf.le.0.or.myf.le.0) call abendm('illegal array size')
      read(1,*) ((f(i,j),i=1,mxf),j=1,myf)
      close(1)
c
      call gparma('Enter CMUP matrix data filename  ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(20,file=fnam2,form='unformatted',status='old')
      read(20) area1, ncc1, iorg1,korg1, ispa1,ispb1,
     *         ixs1, iys1, mszx1, mszy1, mx1, my1, alt1,
     *         ixl,  iyl,  mszx,  mszy,  mx,  my,
     *         nhx, nhy, hw, vd
      if(alt.lt.0.) alt = alt1
      if(area1.ne.area.or.ncc1.ne.ncc.or.alt1.ne.alt .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header inconsistent')
      if(ixs1.ne.ixsf.or.iys1.ne.iysf .or.
     *   mszx1.ne.mszxf.or.mszy1.ne.mszyf .or.
     *   mx1.ne.mxf.or.my1.ne.myf)
     *     call abendm('Header inconsistent')
      write(salt,'(1x,f7.0)') alt
      if (alt.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      xsf = float(ixsf)/1000.
      ysf = float(iysf)/1000.
      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('short of work area')
      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,a8)')
     *     '--- mag. obs.   ', xsf,ysf, mszxf,mszyf, mxf,myf, salt
      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 AMEQ 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()
      call strdtm(sdtm)
c
      call dpcini('.... read CMUP matrix (0) : ')
      do 5 jp=1,myf
      do 5 ip=1,mxf
        read(20) i,j, ic,jc, (g(n),n=1,nnfw)
        if(i.ne.ip.or.j.ne.jp) call abendm('CMUP file broken')
        call dpcent((jp-1)*mxf+ip, mxf*myf)
    5 continue
      rewind(20)
c
      write(6,'(a)')   '----------------------------------------'
      write(6,'(a)')   'AMEQC : <ameq> inversion until converge'
      write(6,'(a,a)')
     *     ' AM mesh data infile : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' CMUP matrix infile  : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' AMEQ model i/o file : ', fnam3(1:lrtrim(fnam3))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xsf, ysf, mszxf, mszyf, mxf, myf, salt,
     *              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, LXO,r,mxf,myf)
      n = 0
      ss = 0.
      do 10 j=1,myf
      do 10 i=1,mxf
        if(f(i,j).eq.vnul) then
          r(i,j) = 0.
        else
          r(i,j) = f(i,j) - r(i,j)
          ss = ss + r(i,j)*r(i,j)
          n = n + 1
        endif
   10 continue
      dv = sqrt(ss / float(n))
      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, LXO,r,mxf,myf, 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, LXO,q,mxf,myf)
        pqq = 0.
        do 40 j=1,myf
        do 40 i=1,mxf
          pqq = pqq + q(i,j)*q(i,j)
   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))
        n = 0
        ss = 0.
        do 60 j=1,myf
        do 60 i=1,mxf
          if(f(i,j).ne.vnul) then
            r(i,j) = r(i,j) - a*q(i,j)
            ss = ss + r(i,j)*r(i,j)
            n = n + 1
          endif
   60   continue
        dv = sqrt(ss / float(n))
        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,access='append')
        write(99,'(//a)') 'AMEQC : <ameq> inversion until converge'
        write(99,'(a,a)')
     *       ' AM mesh data infile : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' CMUP matrix infile  : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' AMEQ model i/o file : ', fnam3(1:lrtrim(fnam3))
        write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xsf, ysf, mszxf, mszyf, mxf, myf, salt,
     *                 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   call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        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, 'mag. obs.   ', 2f10.2, 2i8, 2i6, a8/
     *       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, lrx,r,mx,my)
      dimension g(nfx,nfy)
      dimension s(lsx,mym), r(lrx,my)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call dpcini('.... read CMUP matrix (1) : ')
      read(20) area, (i,k=1,22)
      do 1 jp=1,my
      do 1 ip=1,mx
        read(20) k,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(ip,jp) = rp
        call dpcent((jp-1)*mx+ip, mx*my)
    1 continue
      rewind(20)
      return
      end
c---
      subroutine conv2(g,nfx,nfy, lrx,r,mx,my, lsx,s,mxm,mym)
      dimension g(nfx,nfy)
      dimension s(lsx,mym), r(lrx,my)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call dpcini('.... read CMUP matrix (2) : ')
      read(20) area, (i,k=1,22)
      do 1 js=1,mym
      do 1 is=1,mxm
        s(is,js) = 0.
    1 continue
      do 2 jp=1,my
      do 2 ip=1,mx
        read(20) k,k, ic,jc, g
        kis = ic - nhx
        kjs = jc - nhy
        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(ip,jp)
    4     continue
    3   continue
        call dpcent((jp-1)*mx+ip, mx*my)
    2 continue
      rewind(20)
      return
      end
