c--------( AMAG : <amag> inversion          )--------
      parameter (LHW=100, LXO=1201, LYO=1201)
      parameter (LXM=LXO+LHW+LHW, LYM=LYO+LHW+LHW)
c
      dimension g((LHW*2+1)*(LHW*2+1)), coef(3)
      dimension f(LXO,LYO), r(LXO,LYO), q(LXO,LYO)
      dimension s(LXM,LYM), p(LXM,LYM), v(LXM,LYM)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character buf*80, area*8, area1*8, out*72, salt*8, yn*1, sdtm*20
c
      call premsg('AMAG : <amag> 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 Mag.Anomaly data filename   ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(12,file=fnam1,status='old')
    1 read(12,'(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(12,'(a)') buf
      read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
      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(12,*) ((f(i,j),i=1,mx),j=1,my)
      close(12)
      call prompt('  Remove Linear trend ?  ')
      call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
      if(yn.ne.'y'.and.yn.ne.'n') call abendm('invalid parameter')
c
      call gparma('Enter COEF 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,  mszxm, mszym, mxm, mym, nhx, nhy
      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(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('Header inconsistent')
      if(alt.lt.0.) alt = alt1
      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.
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      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.   ', xs, ys, mszx, mszy, mx, my, salt
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- source model', xl,yl, mszxm,mszym, mxm,mym, '    --- '
      call premsg(out)
      if(mxm.gt.LXM.or.mym.gt.LYM) call abendm('too large array size')
      hw = float(mszx*nhx) / 1000.
      nfx = nhx + 1 + nhx
      nfy = nhy + 1 + nhy
      nnfw = nfx * nfy
c
      call gparma('Enter AMAG model in/out filename  ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(21,file=fnam3,status='old')
    2 read(21,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      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,*) ixl1, iyl1, mszxm1, mszym1, 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.ixl.or.iyl1.ne.iyl .or.
     *   mszxm1.ne.mszxm.or.mszym1.ne.mszym .or.
     *   mxm1.ne.mxm.or.mym1.ne.mym)
     *     call abendm('Header inconsistent')
      read(21,*) ((s(i,j),i=1,mxm),j=1,mym)
      rewind(21)
c
      call premsg('Enter Auxiliary(2) output filename')
      call gparma('      ( blank to suppress output ) ==> ',
     *            81-ldr, fnam4(ldr:80))
      kout = 0
      if(fnam4(ldr:ldr).ne.' ') then
        open(8,file=fnam4,status='new')
        kout = 1
      endif
      call premsg('How many loops to excute ?')
      call gparmi('   (minus to restart from loop-0) ==> ', nloop)
      nlp = iabs(nloop)
      if (nloop.lt.0) then
        call gparmf(' Initial value of magnetiz. (in A/m) ==> ', vm)
      endif
      call clspin()
      call strdtm(sdtm)
c
      call dpcini('.... read COEF matrix (0) : ')
      do 100 j=1,mym
      do 100 i=1,mxm
        read(20) ip, jp, k, k, (c,k=1,nnfw)
        if(ip.ne.i.or.jp.ne.j) call abendm('COEF file broken')
        p(i,j) = 0.
        call dpcent((j-1)*mxm+i, mxm*mym)
  100 continue
      rewind(20)
c
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'AMAG : <amag> inversion'
      write(6,'(a,a)')
     *     ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' COEF matrix infile  : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' AMAG model i/o file : ', fnam3(1:lrtrim(fnam3))
      if(kout.eq.1) write(6,'(a,a)')
     *              ' Aux(2) output file  : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs, ys, mszx,  mszy,  mx,  my,  salt,
     *              xl, yl, mszxm, mszym, mxm, mym, '    --- '
      write(6,2000) hw, nhx, nhy
      if (nloop.lt.0) then
        write(6,'(a,f6.2,a)')
     *       ' Starting from Initial magnetization :', vm, ' A/m'
      endif
c
      n = 0
      ss = 0.
      do 200 j=1,my
      do 200 i=1,mx
        if(f(i,j).eq.vnul) goto 200
        ss = ss + f(i,j)
        n = n + 1
  200 continue
      ft = ss / float(n)
      ss = 0.
      do 220 j=1,mym
      do 220 i=1,mxm
        ss = ss + s(i,j)*s(i,j)
  220 continue
      ds = sqrt(ss / float(mxm*mym))
c
      if(nloop.lt.0) then
        klp = 0
        do 300 j=1,mym
        do 300 i=1,mxm
          s(i,j) = vm * 100.
  300   continue
        ds = 0.
      endif
c
      if (yn.eq.'y') then
        cmx = float(mx+1) / 2.
        cmy = float(my+1) / 2.
        call sm2opn(1, 1)
        do 400 j=1,my
          y = float(j) - cmy
          do 400 i=1,mx
            x = float(i) - cmx
            if(f(i,j).ne.vnul) call sm2ex(x, y, f(i,j)-ft)
  400   continue
        call sm2cls(coef, 3, 1)
        cf = coef(1) + ft
        af = coef(2)
        bf = coef(3)
        write(6,3000) cf, af, bf
        do 500 j=1,my
          y = float(j) - cmy
          do 500 i=1,mx
            x = float(i) - cmx
            if(f(i,j).ne.vnul) f(i,j) = f(i,j) - (af*x + bf*y + cf)
  500   continue
        if(kout.eq.1) then
          write(8,'(a8,i4,4x,4i8)')
     *         'Obs-Trnd', ncc, iorg,korg, ispa,ispb
          write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *         ixs,iys, mszx,mszy, mx,my, vnul, -1.
          do 550 j=1,my
            write(8,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  550     continue
          write(6,'(a)')
     *      'AUX2 <== Obs-Trnd : Trend-removed Observation'
        endif
      else
        af = 0.
        bf = 0.
        cf = 0.
      endif
c
      call conv1(g, nfx, nfy, LXM, s, mxm, mym, LXO, r, mx, my)
      n = 0
      ss = 0.
      do 10 j=1,my
      do 10 i=1,mx
        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))
      write(6,4000) klp, ds, dv
      open(9,status='scratch')
      write(9,*) klp, ds, dv
c
      pvv = 1.
      do 600 l=1,nlp
        klp = klp + 1
        pvvp = pvv
        call conv2(g, nfx, nfy, LXO, r, mx, my, 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, mx, my)
        pqq = 0.
        do 40 j=1,my
        do 40 i=1,mx
          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)
          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,my
        do 60 i=1,mx
          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))
        write(6,4000) klp, ds, dv
        write(9,*) klp, ds, dv
  600 continue
      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)')
     *     ixl, iyl, mszxm, mszym, mxm, mym, anul, -1., klp
      do 700 j=1,mym
        do 710 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
  710   continue
        write(21,'((f7.1,9(1x,f7.1)))') (s(i,j),i=1,mxm)
  700 continue
      close(21)
      if(kout.eq.1) then
        call conv1(g, nfx, nfy, LXM, s, mxm, mym, LXO, r, mx, my)
        do 800 j=1,my
          y = float(j) - cmy
          do 800 i=1,mx
            x = float(i) - cmx
            if(f(i,j).eq.vnul) then
              r(i,j) = vnul
              q(i,j) = vnul
            else
              q(i,j) = r(i,j) + (af*x + bf*y + cf)
              f(i,j) = f(i,j) - r(i,j)
            endif
  800   continue
        write(8,'(a8,i4,4x,4i8)') 'SynthMdl', ncc, iorg,korg, ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 810 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (r(i,j),i=1,mx)
  810   continue
        write(6,'(a)') 'AUX2 <== SynthMdl : Synthetic model anomaly '
        if (yn.eq.'y') then
          write(8,'(a8,i4,4x,4i8)')
     *         'Mdl+Trnd', ncc, iorg,korg, ispa,ispb
          write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *         ixs,iys, mszx,mszy, mx,my, vnul, -1.
          do 820 j=1,my
            write(8,'((f7.1,9(1x,f7.1)))') (q(i,j),i=1,mx)
  820     continue
          write(6,'(a)') 'AUX2 <== Mdl+Trnd : Trend-added Synth. model'
        endif
        write(8,'(a8,i4,4x,4i8)') 'Residual', ncc, iorg,korg, ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 830 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  830   continue
        write(6,'(a)')
     *       'AUX2 <== Residual : Residual <(Obs.) - (Model)>'
        close(8)
      endif
      close(20)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)') 'AMAG : <amag> inversion'
        write(99,'(a,a)')
     *       ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' COEF matrix infile  : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' AMAG model i/o file : ', fnam3(1:lrtrim(fnam3))
        if(kout.eq.1) write(99,'(a,a)')
     *       ' Aux(2) output file  : ', fnam4(1:lrtrim(fnam4))
        write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xs, ys, mszx,  mszy,  mx,  my,  salt,
     *                 xl, yl, mszxm, mszym, mxm, mym, '    --- '
        write(99,2000) hw, nhx, nhy
        if (nloop.lt.0) then
          write(99,'(a,f6.2,a)')
     *         ' Starting from Initial magnetization :', vm, ' A/m'
        endif
        write(99,3000) cf, af, bf
   92   read(9,*,end=93) klp, ds, dv
        write(99,4000) klp, 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, 'source model', 2f10.2, 2i8, 2i6, a8)
 2000 format(' Half width of calc.window :', f6.2, ' km (', i4,
     *       ' *', i4, '  mag.mesh )')
 3000 format(' Trend surface in nT :', f8.2, ' + (', f7.3, ' * x) + (',
     *       f7.3, ' * y)'/ 5x, '( x/y : N/E coord. in mesh unit ',
     *       'referred to the center of area )')
 4000 format(' loop', i4, '   rms.Step :', f9.3,
     *       '   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,*), r(lrx,*)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call dpcini('.... read COEF matrix (1) : ')
      read(20) area, (i,k=1,20)
      do 1 jp=1,my
      do 1 ip=1,mx
        r(ip,jp) = 0.
    1 continue
      do 2 js=1,mym
      do 2 is=1,mxm
        read(20) k, k, kic, kjc, g
        kis = kic - nhx
        kjs = kjc - nhy
        do 3 j=1,nfy
          jp = kjs + j
          if(jp.le.0.or.jp.gt.my) goto 3
          do 4 i=1,nfx
            ip = kis + i
            if(ip.le.0.or.ip.gt.mx) goto 4
            r(ip,jp) = r(ip,jp) + g(i,j)*s(is,js)
    4     continue
    3   continue
        call dpcent((js-1)*mxm+is, mxm*mym)
    2 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,*), r(lrx,*)
      character*8 area
      nhx = (nfx-1) / 2
      nhy = (nfy-1) / 2
      call dpcini('.... read COEF matrix (2) : ')
      read(20) area, (i,k=1,20)
      do 1 js=1,mym
      do 1 is=1,mxm
        read(20) k, k, kic, kjc, g
        kis = kic - nhx
        kjs = kjc - nhy
        ss = 0.
        do 2 j=1,nfy
        jp = kjs + j
        if(jp.le.0.or.jp.gt.my) goto 2
          do 3 i=1,nfx
            ip = kis + i
            if(ip.le.0.or.ip.gt.mx) goto 3
            ss = ss + g(i,j)*r(ip,jp)
    3     continue
    2   continue
        s(is,js) = ss
        call dpcent((js-1)*mxm+is, mxm*mym)
    1 continue
      rewind(20)
      return
      end
