!--------( EMAGF : preparation for <amag> inversion )--------
  implicit none
  integer,parameter :: KVER = 2018
  integer :: lwkdir, lrtrim;  real :: calma
  real,allocatable :: f(:,:), h(:,:), t(:,:), g(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, buf
  character(72) :: out;  character(25) :: strt;  character(20) :: strs, sdtm
  character(8) :: area, tarea, area1, salt
  integer :: ldr, i, j, nhx, nhy, ks, nmsx, nmsy, ixc, iyc, ia, ja
  integer :: kic, kjc, ki, kj, ir, jr, is, it, js, jt, ixt, iyt
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nct, ncct, iorgt, korgt, ispat, ispbt
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixst, iyst, mszxt, mszyt, mxt, myt
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: ixl, iyl, mszxm, mszym, mxm, mym
  real :: xs, ys, smx, smy, xst, yst, smxt, smyt, xt, yt, xtt, ytt
  real :: vnul, alt, tnul, hnul, xl, yl, xh, yh, smxm, smym
  real :: hw, botm, dip, dec, the, phi, vm, am
  real :: x, y, hc, za, sz, hnmsx, hnmsy

  call premsg('EMAGF : preparation for <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')
  do
    read(12,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
    if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
    nc = ncc
    if (ncc < 200) then
      nc = ncc; ncc = nc + 800
    else
      ncc = ncc - 200; nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  endif
  if ((nc >= 1) .and. (nc <= 62)) then
    iorg = 0;  ispa = 0;  ispb = 0
    korg = (nc-30)*360 - 180
    if (nc > 60) korg = 0
    if (nc == 61) iorg = +5400
    if (nc == 62) iorg = -5400
  endif
  read(12,'(a)') buf
  read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
  write(salt,'(1x,f7.0)') alt
  if (alt == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  xs = float(ixs)/1000.;  ys = float(iys)/1000.
  smx = float(mszx);  smy = float(mszy)

  allocate(f(mx,my));  allocate(h(mx, my))

  read(12,*) ((f(i,j),i=1,mx),j=1,my)

  call gparma('Enter Src.Altitude data filename  ==> ', 81-ldr, fnam2(ldr:80))
  open(10,file=fnam2,status='old')
  do
    read(10,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    if ((ncct<0) .or. (ncct>=400)) call abendm('invalid coord.no.')
    nct = ncct
    if (ncct < 200) then
      nct = ncct; ncct = nct + 800
    else
      ncct = ncct - 200; nct = ncct
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    nct = ncct
    if ((ncct >= 800) .and. (ncct < 1000)) nct = ncct - 800
  endif
  if ((nct >= 1) .and. (nct <= 62)) then
    iorgt = 0;  ispat = 0;  ispbt = 0
    korgt = (nct-30)*360 - 180
    if (nct > 60) korgt = 0
    if (nct == 61) iorgt = +5400
    if (nct == 62) iorgt = -5400
  endif
  read(10,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
  xst = float(ixst)/1000.;  yst = float(iyst)/1000.
  smxt = float(mszxt);  smyt = float(mszyt)

  allocate(t(mxt,myt))

  read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
  close(10)

  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  if ((ncct /= ncc) .or. (iorgt /= iorg) .or. (korgt /= korg)   &
      .or. (ispat /= ispa) .or. (ispbt /= ispb)) then
    write(out,'(12x,a8,12x,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    call premsg(out)
    call abendm('Header inconsistent')
  endif
  if (tarea /= area) call premsg('            ' // tarea)
  write(out,'(a14,2a10,2x,2a8,2a6,a8)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '--- source alt. ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
  call premsg(out)
  xt = xs + float(mszx*(mx-1))/1000.
  yt = ys + float(mszy*(my-1))/1000.
  xtt = xst + float(mszxt*(mxt-1))/1000.
  ytt = yst + float(mszyt*(myt-1))/1000.
  if ((xst > xs) .or. (yst > ys) .or.   &
      (xtt < xt) .or. (ytt < yt)) call abendm('ranges conflict')

  if (alt == 0.) then
    do
      read(12,'(a)') buf
      if (buf(1:1) /= '#') exit
    enddo
    if (buf(13:16) == '    ') then        !-- v2005
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
      nc1 = ncc1
      if (ncc1 < 200) then
        nc1 = ncc1; ncc1 = nc1 + 800
      else
        ncc1 = ncc1 - 200; nc1 = ncc1
      endif
    else                                  !-- v2018
      read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
    endif
    if ((nc1 >= 1) .and. (nc1 <= 62)) then
      iorg1 = 0;  ispa1 = 0;  ispb1 = 0
      korg1 = (nc1-30)*360 - 180
      if (nc1 > 60) korg1 = 0
      if (nc1 == 61) iorg1 = +5400
      if (nc1 == 62) iorg1 = -5400
    endif
    if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg)   &
        .or. (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
      call abendm('Header1 inconsistent')
    read(12,*) ixs1, iys1, mszx1, mszy1, mx1, my1, hnul
    if ((ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx)   &
        .or. (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
      call abendm('Header2 inconsistent')
    read(12,*) ((h(i,j),i=1,mx),j=1,my)
    strt = ' given as 2nd set data   '
  else
    if (alt < 0.) call abendm('Undefined Altitude')
    write(strt,'(f7.0,a)') alt, ' m above sea level'
  endif
  close(12)
  call premsg('  obs.surface : ' // strt)

  call premsg('')
  call premsg('Specify Source distribution Parameters')
  call premsg('  Southwest corner Coord. in km ?')
  call gparmf('    Northing : ', xl)
  ixl = ixst + mszxt*nint((xl-xst)*1000./smxt);  xl = float(ixl)/1000.
  if ((xl < xst) .or. (xl > xs)) call abendm('ranges conflict')
  call gparmf('    Easting  : ', yl)
  iyl = iyst + mszyt*nint((yl-yst)*1000./smyt);  yl = float(iyl)/1000.
  if ((yl < yst) .or. (yl > ys)) call abendm('ranges conflict')
  call premsg('  Mesh size and numbers to N ?')
  call gparmi('    Unit Mesh size in m :  ', mszxm)
  if (mszxm <= 0) call abendm('invalid value')
  nmsx = mszxm / mszxt
  if((nmsx*mszxt) /= mszxm) call abendm('mesh size illegal')
  hnmsx = float(nmsx-1)/2.
  smxm = float(mszxm)
  call gparmi('    Number of Mesh to N :  ', mxm)
  xh = xl + float(mszxm*(mxm-1))/1000.
  if ((xh < xt) .or. (xh > xtt)) call abendm('ranges conflict')
  call premsg('  Mesh size and numbers to E ?')
  call gparmi('    Unit Mesh size in m :  ', mszym)
  if (mszym <= 0) call abendm('invalid value')
  nmsy = mszym / mszyt
  if((nmsy*mszyt) /= mszym) call abendm('mesh size illegal')
  hnmsy = float(nmsy-1)/2.
  smym = float(mszym)
  call gparmi('    Number of Mesh to E :  ', mym)
  yh = yl + float(mszym*(mym-1))/1000.
  if ((yh < yt) .or. (yh > ytt)) call abendm('ranges conflict')
  call premsg('  Truncation of Source effect ?')
  call gparmf('    Half width of calc.window (in km) :  ', hw)
  nhx = nint(hw*1000./smx);  nhy = nint(hw*1000./smy)
  if ((nhx < 10) .or. (nhy < 10)) call abendm('truncation too short')

  allocate(g(-nhx:+nhx, -nhy:+nhy))

  call premsg('')
  call premsg('Specify Source-Bottom configuration')
  call prompt('    Constant Thickness (1) or Flat Bottom (0) ')
  call gparmi(': ', ks)
  if ((ks < 0) .or. (ks > 1)) call abendm('invalid value')
  strs = ' below sea level    '
  if (ks == 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 < -90.) .or. (dip > 90.) .or.   &
      (dec < -180.) .or. (dec >= 360.)) call abendm('invalid value')
  call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
  if ((the < -90.) .or. (the > 90.) .or.   &
      (phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')
  call gparmf('Enter Initial source magnetization (A/m) ==> ', vm)

  call premsg('')
  call gparma('Enter COEF matrix output filename ==> ', 81-ldr, fnam3(ldr:80))
  open(20,file=fnam3,form='unformatted',status='new')
  call gparma('Enter AMAG model output filename  ==> ', 81-ldr, fnam4(ldr:80))
  call clspin()
  call strdtm(sdtm)
  open(21,file=fnam4,status='new')
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'EMAGF : preparation for <amag> inversion'
  write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Src. Alt. infile    : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' COEF matrix outfile : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' AMAG model outfile  : ', fnam4(1:lrtrim(fnam4))
  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'source alt. ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'source model', xl, yl, mszxm, mszym, mxm, mym, '    --- '
  write(6,'(a,f6.2,a,i4,a,i4,a)')   &
    ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mag.mesh)'
  write(6,'(3x,a,a25)') 'obs. surface  : ', strt
  write(6,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'ambient field : ', 'dip=', dip, 'dec=', dec
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'magnetization : ', 'dip=', the, 'dec=', phi
  write(6,'(a,f8.2)') ' Initial source magnetization (A/m) :', vm

  call dpcini('... gen. COEF matrix data : ')
  call magafd(the, phi, dip, dec)
  write(20) area, KVER, ncc, iorg,korg, ispa,ispb,   &
           ixs, iys, mszx,  mszy,  mx,  my,  alt,   &
           ixl, iyl, mszxm, mszym, mxm, mym, nhx, nhy
  am = smxt*smyt

  do j=1,mym
    iyc = iyl + (j-1)*mszym
    ja = nint(float(iyc-iyst)/smyt - hnmsy)
    kjc = nint(float(iyc-iys)/smy)
    do i=1,mxm
      ixc = ixl + (i-1)*mszxm
      ia = nint(float(ixc-ixst)/smxt - hnmsx)
      kic = nint(float(ixc-ixs)/smx)
      do kj=-nhy,+nhy; do ki=-nhx,+nhx;  g(ki,kj) = 0.;  enddo; enddo
      do js=1,nmsy
        jt = ja + js
        if ((jt <= 0) .or. (jt > myt)) cycle
        iyt = iyst + (jt-1)*mszyt
        do is=1,nmsx
          it = ia + is
          if ((it <= 0) .or. (it > mxt)) cycle
          ixt = ixst + (it-1)*mszxt
          if (t(it,jt) == tnul) cycle
          do kj=-nhy,nhy
            jr = kjc + kj
            if ((jr < 0) .or. (jr >= my)) cycle
            y = float(iyt - iys - jr*mszy)
            do ki=-nhx,nhx
              ir = kic + ki
              x = float(ixt - ixs - ir*mszx)
              if ((ir < 0) .or. (ir >= mx)) cycle
              if (f(ir+1,jr+1) == vnul) cycle
              hc = alt
              if (alt == 0.) hc = h(ir+1,jr+1)
              za = hc - t(it,jt);  sz = botm
              if (ks == 0) sz = botm + t(it,jt)
              if ((iabs(ki) <= 10) .and. (iabs(kj) <= 10)) then
                call mprism(0.01, x, smxt, y, smyt, za, sz)
              else
                call mvline(0.01, x, y, za, sz, am)
              endif
              g(ki,kj) = g(ki,kj) + calma(0., 0., 0.)
            enddo
          enddo
        enddo
      enddo
      write(20) i,j, kic,kjc, ((g(ki,kj),ki=-nhx,nhx),kj=-nhy,nhy)
      call dpcent((j-1)*mxm+i, mxm*mym)
    enddo
  enddo
  close(20)

  ixl = nint(xl*1000.)
  iyl = nint(yl*1000.)
  write(21,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)')   &
    ixl, iyl, mszxm, mszym, mxm, mym, 9999.9, -1., 0
  do j=1,mym
    write(21,'((f7.1,9(1x,f7.1)))') ((vm*100.),i=1,mxm)
  enddo
  close(21)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'EMAGF : preparation for <amag> inversion'
    write(99,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Src. Alt. infile    : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' COEF matrix outfile : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' AMAG model outfile  : ', fnam4(1:lrtrim(fnam4))
    write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'source alt. ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'source model', xl, yl, mszxm, mszym, mxm, mym, '    --- '
    write(99,'(a,f6.2,a,i4,a,i4,a)') ' Half width of calc.window :',   &
      hw, ' km (', nhx, ' *', nhy, '  mag.mesh)'
    write(99,'(3x,a,a25)') 'obs. surface  : ', strt
    write(99,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'ambient field : ', 'dip=', dip, 'dec=', dec
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'magnetization : ', 'dip=', the, 'dec=', phi
    write(99,'(a,f8.2)') ' Initial source magnetization (A/m) :', vm
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end