c--------( CMAGF : calc. anomaly of <amag> fine model )--------
      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), e(6)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character buf*80, out*72, area*8, area1*8, area2*8
      character strs*20, sdtm*20
c
      call premsg('CMAGF : calc. anomaly of <amag> fine model')
      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))
c
      call gparma('Enter Calc.Alt. 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,*) ixs, iys, 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')
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      ixt = ixs + mszx*(mx-1)
      iyt = iys + mszy*(my-1)
      read(1,*) ((h(i,j),i=1,mx),j=1,my)
      close(1)
      smx = float(mszx)
      smy = float(mszy)
c
      call gparma('Enter AMAG 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(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,*) ixl, iyl, mszxm, mszym, 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)
     *     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')
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      ixh = ixl + mszxm*(mxm-1)
      iyh = iyl + mszym*(mym-1)
      if(ixl.gt.ixs.or.iyl.gt.iys .or.
     *   ixh.lt.ixt.or.iyh.lt.iyt) call abendm('ranges conflict')
      read(2,*) ((s(i,j),i=1,mxm),j=1,mym)
      close(2)
      smxm = float(mszxm)
      smym = float(mszym)
c
      call gparma('Enter Source Alt. data filename  ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(10,file=fnam3,status='old')
    3 read(10,'(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(10,'(a)') buf
      read(buf,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
      if(ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header inconsistent')
      if(mxt.gt.LXM.or.myt.gt.LYM) call abendm('too large array size')
      if(mxt.lt.10.or.myt.lt.10) call abendm('illegal array size')
      nmsx = mszxm / mszxt
      if((nmsx*mszxt).ne.mszxm) call abendm('mesh size illegal')
      nmsy = mszym / mszyt
      if((nmsy*mszyt).ne.mszym) call abendm('mesh size illegal')
      hnmsx = float(nmsx-1)/2.
      hnmsy = float(nmsy-1)/2.
      xst = float(ixst)/1000.
      yst = float(iyst)/1000.
      ixtt = ixst + mszxt*(mxt-1)
      iytt = iyst + mszyt*(myt-1)
      if(ixst.gt.ixl.or.iyst.gt.iyl .or.
     *   ixtt.lt.ixh.or.iytt.lt.iyh) call abendm('ranges conflict')
      read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
      close(10)
      smxt = float(mszxt)
      smyt = float(mszyt)
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)')
     *     '--- Calc. Alt.  ', xs, ys,  mszx,mszy, mx,my
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- AMAG model  ', xl, yl, mszxm,mszym, mxm,mym
      write(out,'(a,2f10.2,2i8,2i6)')
     *     '--- Source Alt. ', xst, yst, mszxt,mszyt, mxt,myt
      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 premsg('')
      call premsg('Specify Source-Bottom configuration')
      call prompt('    Constant Thickness (1) or Flat Bottom (0) ')
      call gparmi(': ', ks)
      if(ks.lt.0.or.ks.gt.1) call abendm('invalid value')
      strs = ' below sea level    '
      if(ks.eq.1) strs = ' below terr.surface '
      call gparmf('  Bottom-depth (in m)' // strs // ': ', botm)
      call premsg('Specify Direction Parameters (in degrees)')
      call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90. .or.
     *   dec.lt.-180..or.dec.ge.360.) call abendm('invalid value')
      call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
      if(the.lt.-90..or.the.gt.90. .or.
     *   phi.lt.-180..or.phi.ge.360.) call abendm('invalid value')
c
      call gparma('Enter Calc.Anom. output filename ==> ',
     *            81-ldr, fnam4(ldr:80))
      call clspin()
      call strdtm(sdtm)
      open(3,file=fnam4,status='new')
      write(6,'(a)')   '----------------------------------------'
      write(6,'(a)')   'CMAGF : calc. anomaly of <amag> fine model'
      write(6,'(a,a)')
     *     ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' AMAG model infile     : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)')
     *     ' Source Alt. infile    : ', fnam3(1:lrtrim(fnam3))
      write(6,'(a,a)')
     *     ' Calc.Anom. outfile    : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs, ys,  mszx, mszy,  mx, my,
     *              xl, yl,  mszxm,mszym, mxm,mym,
     *              xst,yst, mszxt,mszyt, mxt,myt
      write(6,2000) hw, nhx, nhy
c
      do 100 j=1,my
      do 100 i=1,mx
        f(i,j) = 0.
        if(h(i,j).eq.hnul) f(i,j) = 9999.9
  100 continue
      call dpcini('... calc. model anomaly   : ')
      call magafd(the, phi, dip, dec)
      am = smxt * smyt
      do 200 j=1,mym
        iyc = iyl + (j-1)*mszym
        ja = nint(float(iyc-iyst)/smyt) - hnmsy
        kjc = nint(float(iyc-iys)/smy)
        do 200 i=1,mxm
          ixc = ixl + (i-1)*mszxm
          ia = nint(float(ixc-ixst)/smxt) - hnmsx
          kic = nint(float(ixc-ixs)/smx)
          do 300 js=1,nmsy
            jt = ja + js
            if(jt.le.0.or.jt.gt.myt) goto 300
            iyt = iyst + (jt-1)*mszyt
            do 400 is=1,nmsx
              it = ia + is
              if(it.le.0.or.it.gt.mxt) goto 400
              ixt = ixst + (it-1)*mszxt
              ht = t(ia+1,ja+1)
              if(ht.eq.tnul) call abendm('source alt. undefined')
              do 500 kj=-nhy,nhy
                jr = kjc + kj
                if(jr.lt.0.or.jr.ge.my) goto 500
                y = float(iyt - iys - jr*mszy)
                do 600 ki=-nhx,nhx
                  ir = kic + ki
                  x = float(ixt - ixs - ir*mszx)
                  if(ir.lt.0.or.ir.ge.mx) goto 600
                  hc = h(ir+1,jr+1)
                  if(hc.eq.hnul) goto 600
                  za = hc - ht
                  sz = botm
                  if(ks.eq.0) sz = ht + botm
                  if(iabs(ki).le.10.and.iabs(kj).le.10) then
                    call mprism(s(i,j)/100., x, smxt, y, smyt, za, sz)
                  else
                    call mvline(s(i,j)/100., x, y, za, sz, am)
                  endif
                  f(ir+1,jr+1) = f(ir+1,jr+1) + calma(0.,0.,0.)
  600           continue
  500         continue
  400       continue
  300     continue
          call dpcent((j-1)*mxm+i, mxm*mym)
  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)')
     *     ixs, iys, mszx, mszy, mx, my, 9999.9, 0.
      do 700 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  700 continue
      write(3,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixs, iys, mszx, mszy, mx, my, hnul, -1.
      do 800 j=1,my
        write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  800 continue
      close(3)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)')
     *        'CMAGF : calc. anomaly of <amag> fine model'
        write(99,'(a,a)')
     *       ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' AMAG model infile     : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' Source Alt. infile    : ', fnam3(1:lrtrim(fnam3))
        write(99,'(a,a)')
     *       ' Calc.Anom. outfile    : ', 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,
     *                xl, yl,  mszxm,mszym, mxm,mym,
     *                xst,yst, mszxt,mszyt, mxt,myt
        write(99,2000) hw, nhx, nhy
        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'/
     *       3x, 'Calc. Alt.  ', 2f10.2, 2i8, 2i6/
     *       3x, 'AMAG model  ', 2f10.2, 2i8, 2i6/
     *       3x, 'Source Alt  ', 2f10.2, 2i8, 2i6)
 2000 format(' Half width of calc.window :', f6.2,
     *       ' km (', i4, ' *', i4, '  calc.-mesh )')
      end
