c--------( CALMAS : calc. model anomaly on a given surface )--------
      parameter (LXO=1201, LYO=1201, MNM=20)
      parameter (DPI=3.14159*2.)
c
      dimension f(LXO,LYO), h(LXO,LYO)
      dimension pv(7,MNM), pam(MNM), pth(MNM), pph(MNM), ipv(MNM)
      character*80 wdr, flog, fnam1, fnam2
      character buf*80, out*72, area*8, area1*8, area2*8, sdtm*20
      character*5 nmdl(4)/ 'Block', 'hRect', 'vLine', 'Point'/
c
      call premsg('CALMAS : calc. model anomaly on a given surface')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Calc.Surface Alt. 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,*) ixl, iyl, 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')
      read(1,*) ((h(i,j),i=1,mx),j=1,my)
      close(1)
      xl = float(ixl)/1000.
      yl = float(iyl)/1000.
      do 100 j=1,my
      do 100 i=1,mx
        f(i,j) = 0.
  100 continue
      ixh = ixl + mszx*(mx-1)
      iyh = iyl + mszy*(my-1)
      call gparma('Enter Calc. MA output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(2,file=fnam2,status='new')
c
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'CALMAS : calc. model anomaly on a given surface'
      write(6,'(a,a)')
     *     ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)')
     *     ' Model anomaly outfile  : ', fnam2(1:lrtrim(fnam2))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xl, yl, mszx, mszy, mx, my,
     *              xl, yl, mszx, mszy, mx, my, '    var.'
c
      call premsg('')
      call premsg('How to specify the source location')
      call premsg('  0: Coordinates values   or')
      call gparmi('  1: Distance from SW corner ==> ', kl)
      if(kl.lt.0.or.kl.gt.1) call abendm('invalid selection')
      call premsg('Specify Ambient field Direction')
      call gparmf2('  Inclin., Declin. (in degrees) : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90. .or.
     *   dec.lt.-180..or.dec.ge.360.) call abendm('invalid value')
c
      nm = 0
    2 call premsg('Select source model type:')
      call prompt('  1:Block 2:hRect 3:vLine 4:Point or 0:End')
      call gparmi(' ==> ', ks)
      if(ks.lt.0.or.ks.gt.4) call abendm('invalid selection')
      if(ks.eq.0) goto 3
      nm = nm + 1
      if(nm.gt.MNM) call abendm('too many models')
      ipv(nm) = ks
      call premsg('  Horizontal location:')
      if(ks.eq.1.or.ks.eq.2) then
        call prompt('   Northing Center (cx), N-S size (rx) (in km)')
        call gparmf2(' ==> ', cx, rx)
        if(kl.eq.1) cx = cx + xl
        pv(1,nm) = cx
        pv(2,nm) = rx
        call prompt('   Easting Center (cy), E-W size (ry) (in km) ')
        call gparmf2(' ==> ', cy, ry)
        if(kl.eq.1) cy = cy + yl
        pv(3,nm) = cy
        pv(4,nm) = ry
      else
        call gparmf('   Northing (xc) (in km) ==> ', xc)
        if(kl.eq.1) xc = xc + xl
        pv(1,nm) = xc
        pv(2,nm) = 0.
        call gparmf('   Easting (yc) (in km)  ==> ', yc)
        if(kl.eq.1) yc = yc + yl
        pv(3,nm) = yc
        pv(4,nm) = 0.
      endif
      call premsg('  Depth below Sea-level:')
      if(ks.eq.1.or.ks.eq.3) then
        call prompt('   Top depth (za), and Thickness (rz) (in km)')
        call gparmf2(' ==> ', za, rz)
        pv(5,nm) = za
        pv(6,nm) = rz
      else
        call gparmf('    (in km) ==> ', zc)
        pv(5,nm) = zc
        pv(6,nm) = 0.
      endif
      if(ks.eq.4) then
        call prompt('  Equivalent source volume (in cubic-km)')
        call gparmf(' ==> ', vol)
      else if(ks.eq.3) then
        call prompt('  Equivalent source area (in square-km)')
        call gparmf(' ==> ', vol)
      else if(ks.eq.2) then
        call prompt('  Equivalent source thickness (in km)')
        call gparmf(' ==> ', vol)
      else
        vol = 0.
      endif
      pv(7,nm) = vol
      call gparmf('  Magnetization Intensity (A/m) : ', am)
      pam(nm) = am
      call premsg('  Magnetization Direction ?')
      call gparmf2('    Inclin., Declin. (in degrees) : ', the, phi)
      if(the.lt.-90..or.the.gt.90. .or.
     *   phi.lt.-180..or.phi.ge.360.) call abendm('invalid value')
      pth(nm) = the
      pph(nm) = phi
      write(6,2000) nm, nmdl(ipv(nm)), pam(nm), pth(nm), pph(nm),
     *              (pv(i,nm),i=1,7)
c
      call magafd(the, phi, dip, dec)
      if(ks.eq.4) then
        call mpoint(am, xc, yc, zc, vol)
      else if(ks.eq.3) then
        call mvline(am, xc, yc, za, rz, vol)
      else if(ks.eq.2) then
        call mhrect(am, cx, rx, cy, ry, zc, vol)
      else
        call mprism(am, cx, rx, cy, ry, za, rz)
      endif
      do 200 j=1,my
      do 200 i=1,mx
        if(h(i,j).eq.hnul) then
          f(i,j) = 9999.9
        else
          xp = float(ixl + mszx*(i-1)) / 1000.
          yp = float(iyl + mszy*(j-1)) / 1000.
          zp = - h(i,j) / 1000.
          f(i,j) = f(i,j) + calma(xp, yp, zp)
        endif
  200 continue
      goto 2
c
    3 call clspin()
      write(2,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl, iyl, mszx, mszy, mx, my, 9999.9, 0.
      do 300 j=1,my
        write(2,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  300 continue
      write(2,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl, iyl, mszx, mszy, mx, my, hnul, -1.
      do 400 j=1,my
        write(2,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  400 continue
      close(2)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)')
     *       'CALMAS : calc. model anomaly on a given surface'
        write(99,'(a,a)')
     *       ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' Model anomaly outfile  : ', fnam2(1:lrtrim(fnam2))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xl, yl, mszx, mszy, mx, my,
     *                 xl, yl, mszx, mszy, mx, my, '    var.'
        do 500 n=1,nm
          write(99,2000) n, nmdl(ipv(n)), pam(n), pth(n), pph(n),
     *                   (pv(i,n),i=1,7)
  500   continue
        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, 'calc. surf. ', 2f10.2, 2i8, 2i6/
     *       3x, 'model anom. ', 2f10.2, 2i8, 2i6, a8)
 2000 format(i4, 2x, a5, f8.2, 2f8.1/ 6x, 3(f9.2, f7.2), f9.4)
      end
