!--------( EMEQ : preparation for <ameq> inversion )--------
  implicit none
  integer,parameter :: KVER = 2018
  integer :: lwkdir, lrtrim
  real,parameter :: DPI = 3.141593*2.
  real,allocatable :: f(:,:), h(:,:), t(:,:), tt(:,:), g(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, buf
  character(72) :: out;  character(25) :: strt;  character(20) :: sdtm
  character(8) :: area, tarea, area1, salt
  integer :: ldr, i, j, nhx, nhy, nfx, nfy, ic, jc, ip, jp, is, js
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixsf, iysf, mszxf, mszyf, mxf, myf
  integer :: ixl, iyl, mszx, mszy, mx, my
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: ixlm, iylm, mxm, mym, ki, kj
  real :: vnul, alt, xsf, ysf, tnul, xl, yl, vd, hw, hwm, xlm, ylm, hnul
  real :: fmx, fmy, xtf, ytf, smx, smy, xh, yh, amxy, hhwm4, ss
  real :: x, y, rr, w, s, xc, yc, xp, yp, hp, ht, z, r

  call premsg('EMEQ : preparation for <ameq> 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 AM mesh data filename    ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  do
    read(1,'(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(1,'(a)') buf
  read(buf,*) ixsf, iysf, mszxf, mszyf, mxf, myf, vnul, alt
  write(salt,'(1x,f7.0)') alt
  if (alt == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  xsf = float(ixsf)/1000.;  ysf = float(iysf)/1000.

  allocate(f(mxf,myf))

  read(1,*) ((f(i,j),i=1,mxf),j=1,myf)

  call gparma('Enter Reduction Alt. filename  ==> ', 81-ldr, fnam2(ldr:80))
  open(2,file=fnam2,status='old')
  do
    read(2,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') tarea, 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)') tarea, 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
  read(2,'(a)') buf
  read(buf,*) ixl, iyl, mszx, mszy, mx, my, tnul
  xl = float(ixl)/1000.;  yl = float(iyl)/1000.

  allocate(t(mx,my))

  read(2,*) ((t(i,j),i=1,mx),j=1,my)
  close(2)

  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  if ((ncc1 /= ncc) .or. (iorg1 /= iorg) .or. (korg1 /= korg)   &
      .or. (ispa1 /= ispa) .or. (ispb1 /= ispb)) then
    write(out,'(12x,a8,12x,i4,4x,4i8)') tarea, ncc1, iorg1,korg1, ispa1,ispb1
    call premsg(out)
    call abendm('Header inconsistent')
  endif
  if (area /= tarea) 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.   ', xsf, ysf, mszxf, mszyf, mxf, myf, salt
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '--- reduc. alt. ', xl,  yl,  mszx,  mszy,  mx,  my, '    --- '
  call premsg(out)
  fmx = float(mszxf)/1000.;  fmy = float(mszyf)/1000.
  xtf = xsf + float(mxf-1)*fmx;  ytf = ysf + float(myf-1)*fmy
  smx = float(mszx);  smy = float(mszy)
  xh = xl + float(mx-1)*smx/1000.;  yh = yl + float(my-1)*smy/1000.
  amxy = smx * smy
  if ((xl > xsf) .or. (yl > ysf) .or.   &
      (xh < xtf) .or. (yh < ytf)) call abendm('ranges conflict')

  if (alt == 0.) then
    do
      read(1,'(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(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, hnul
    if ((ixs1 /= ixsf) .or. (iys1 /= iysf) .or. (mszx1 /= mszxf)   &
        .or. (mszy1 /= mszyf) .or. (mx1 /= mxf) .or. (my1 /= myf))   &
      call abendm('Header2 inconsistent')

    allocate(h(mxf, myf))

    read(1,*) ((h(i,j),i=1,mxf),j=1,myf)
    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(1)
  call premsg('  obs.surface : ' // strt)

  call premsg('')
  call premsg('Specify Equivalent Source Parameter')
  call premsg('  Vertical distance from reduction surface ?')
  call gparmf('    Downward (in m)    : ', vd)
  if (vd <= 0.) call abendm('invalid value')
  call premsg('  Truncation of Source effect ?')
  call gparmf('    Half width of calc.window (in km) : ', hw)
  hwm = hw * 1000.;  nhx = nint(hwm/smx);  nhy = nint(hwm/smy)
  if ((nhx < 10) .or. (nhy < 10)) call abendm('truncation too short')
  nfx = nhx + nhx;  nfy = nhy + nhy
  mxm = mx + nhx + nhx;  mym = my + nhy + nhy
  ixlm = ixl - mszx*nhx;  iylm = iyl - mszy*nhy
  xlm = float(ixlm)/1000.;  ylm = float(iylm)/1000.

  allocate(tt(mxm,mym));  allocate(g(-nhx:+nhx,-nhy:+nhy))

  call gparma('Enter CMUP matrix out filename ==> ', 81-ldr, fnam3(ldr:80))
  open(20,file=fnam3,form='unformatted',status='new')
  call gparma('Enter AMEQ model out filename  ==> ', 81-ldr, fnam4(ldr:80))
  open(21,file=fnam4,status='new')
  call clspin()
  call strdtm(sdtm)
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'EMEQ : preparation for <ameq> inversion'
  write(6,'(a,a)') ' AM mesh data infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' CMUP matrix outfile   : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' AMEQ 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.   ', xsf, ysf, mszxf, mszyf, mxf, myf, salt
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'reduc. surf.', xl,  yl,  mszx,  mszy,  mx,  my, '    --- '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'eq.src.anom.', xlm, ylm, mszx,  mszy,  mxm, mym, '    var.'
  write(6,'(a,f6.2,a,i4,a,i4,a)')   &
    ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mesh )'
  write(6,'(a,f8.1,a)')   &
    ' Equiv. source anom. :', vd, ' m below Reduction alt.'

  call dpcini('... gen. Source alt. data : ')
  do jp=1,mym; do ip=1,mxm;  tt(ip,jp) = tnul;  enddo; enddo
  do j=1,my; do i=1,mx
    if (t(i,j) /= tnul) tt(i+nhx,j+nhy) = t(i,j) - vd
  enddo; enddo
  hhwm4 = hwm*hwm * 4.
  do jp=1,mym; do ip=1,mxm
    if (tt(ip,jp) == tnul) then
      ss = hhwm4
      do js=-nhy,+nhy
        j = jp - nhy + js
        if ((j <= 0) .or. (j > my)) cycle
        y = float(js)*smy
        do is=-nhx,+nhx
          i = ip - nhx + is
          if ((i <= 0) .or. (i > mx)) cycle
          if (t(i,j) == tnul) cycle
          x = float(is)*smx
          rr = x*x + y*y
          if (rr < ss) ss = rr
        enddo
      enddo
      if (ss < hhwm4) then
        ss = ss + ss;  w = 0.;  s = 0.
        do js=-nfy,+nfy
          j = jp - nhy + js
          if ((j <= 0) .or. (j > my)) cycle
          y = float(js)*smy
          do is=-nfx,+nfx
            i = ip - nhx + is
            if ((i <= 0) .or. (i > mx)) cycle
            if (t(i,j) == tnul) cycle
            x = float(is)*smx
            rr = x*x + y*y
            if (rr < ss) then
              s = s + t(i,j)*(ss-rr)
              w = w + (ss-rr)
            endif
          enddo
        enddo
        tt(ip,jp) = s/w - vd
      endif
    endif
    call dpcent((jp-1)*mxm+ip, mxm*mym)
  enddo; enddo

  call dpcini('... gen. CMUP matrix data : ')
  write(20) area, KVER, ncc, iorg,korg, ispa,ispb,   &
           ixsf, iysf, mszxf, mszyf, mxf, myf, alt,   &
           ixl,  iyl,  mszx,  mszy,  mx,  my,   &
           nhx, nhy, hw, vd
  do j=1,myf
    yc = ysf + fmy*float(j-1)
    yp = (yc-ylm)*1000.;  jc = nint(yp/smy)
    do i=1,mxf
      xc = xsf + fmx*float(i-1)
      xp = (xc-xlm)*1000.;  ic = nint(xp/smx)
      do kj=-nhy,+nhy; do ki=-nhx,+nhx;  g(ki,kj) = 0.;  enddo; enddo
      hp = alt
      if (alt == 0.) hp = h(i,j)
      if (hp /= hnul) then
        do kj=-nhy,+nhy
          js = jc + kj
          if ((js < 0) .or. (js >= mym)) cycle
          y = float(js)*smy - yp
          do ki=-nhx,+nhx
            is = ic + ki
            if ((is < 0) .or. (is >= mxm)) cycle
            ht = tt(is+1,js+1)
            if (ht /= tnul) then
              x = float(is)*smx - xp;  z = hp - ht
              rr = x*x + y*y + z*z; r = sqrt(rr)
              g(ki,kj) = z*amxy / (rr*r*DPI)
            endif
          enddo
        enddo
      endif
      write(20) i,j, ic,jc, ((g(ki,kj),ki=-nhx,+nhx),kj=-nhy,+nhy)
      call dpcent((j-1)*mxf+i, mxf*myf)
    enddo
  enddo
  close(20)

  write(21,'(a8,4x,i4,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., 0
  do j=1,mym
    write(21,'((f7.1,9(1x,f7.1)))') (0.,i=1,mxm)
  enddo
  write(21,'(a8,4x,i4,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 j=1,mym
    write(21,'((f7.1,9(1x,f7.1)))') (tt(i,j),i=1,mxm)
  enddo
  close(21)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'EMEQ : preparation for <ameq> inversion'
    write(99,'(a,a)') ' AM mesh data infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' CMUP matrix outfile   : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' AMEQ 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.   ', xsf, ysf, mszxf, mszyf, mxf, myf, salt
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'reduc. surf.', xl,  yl,  mszx,  mszy,  mx,  my, '    --- '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'eq.src.anom.', xlm, ylm, mszx,  mszy,  mxm, mym, '    var.'
    write(99,'(a,f6.2,a,i4,a,i4,a)')   &
      ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  mesh )'
    write(99,'(a,f8.1,a)')   &
      ' Equiv. source anom. :', vd, ' m below Reduction alt.'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end