!--------( EXDEQ3 : preparation for Gen.Mis-tie control <axdeq> )--------
!--------(            to Control on all crossover points )--------
  implicit none
  integer,parameter :: KVER = 2018
  real,parameter :: DPI = 3.141593 * 2.
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), t(:,:), tt(:,:), g(:,:), c(:)
  integer,allocatable :: iq(:), mv(:)
  character(8),allocatable :: sl(:), sl1(:), sl2(:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, fnam5, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, sln, sla, slb;  character(1) :: cln
  integer :: ldr, i, j, n, nhx, nhy, nfx, nfy, ip, jp, k, IER
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ixl, iyl, mszx, mszy, mx, my
  integer :: ixlm, iylm, mxm, mym, is, js, ic, jc, ki, kj
  integer :: mln, mdpt, mcp, nln, ndpt, iln, iql, npv
  integer :: ka, kb, kcl, kp, mva, mvb, mvs, na, nb, kx
  real :: tnul, xl, yl, alt, smx, smy, amxy, vd, hw, hwm, hhwm4, ss
  real :: xlm, ylm, x, y, rr, w, s, z, ht, r, fg, gn, zd
  real :: f1a1, f1a2, f2a1, f2a2, f1b1, f1b2, f2b1, f2b2
  real(8) :: dfi, dfk, dxn, dye

  call premsg('EXDEQ3 : preparation for Gen.Mis-tie control <axdeq>')
  call premsg('           to Control on all crossover points')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr
  fnam3 = wdr;  fnam4 = wdr;  fnam5 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter StdLIN data filename     ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  mln = 1;  mdpt = 0
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) then
      mln = mln + 1
    else
      if (buf(1:1) /= '#') mdpt = mdpt + 1
    endif
  enddo
  mcp = mln*mln/2
  rewind(1)

  allocate(f(3,mdpt));  allocate(c(mcp))
  allocate(iq(mln));  allocate(mv(mcp))
  allocate(sl(mln));  allocate(sl1(mcp));  allocate(sl2(mcp))

  call gparma('Enter Reduction Alt. filename  ==> ', 81-ldr, fnam2(ldr:80))
  open(2,file=fnam2,status='old')
  read(2,'(a)') buf
  do while (buf(1:1) == '#')
    call premsg('  ' // buf(1:lrtrim(buf)))
    read(2,'(a)') buf
  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(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)

  call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
  nln = 1;  iq(nln) = 0;  sl(nln) = '........';  ndpt = 0
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) then
      if (iq(nln) /= 0) nln = nln + 1
      iq(nln) = 0;  sln = buf(2:9)
      if (sln == '        ') call abendm('line name missing')
      do while (sln(1:1) == ' ');  sln = sln(2:8);  enddo
      sl(nln) = sln
      do i=1,(nln-1)
        if (sl(i) == sln) call abendm('line name repeated')
      enddo
    else if (buf(1:1) /= '#') then
      ndpt = ndpt + 1
      i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
      i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
      i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
      read(buf,*) dfi, dfk, alt
      call cvdiken(dfi, dfk, dye, dxn)
      f(1,ndpt) = real(dxn);  f(2,ndpt) = real(dye);  f(3,ndpt) = alt
      iq(nln) = ndpt
    endif
  enddo
  close(1)
  if (iq(nln) == 0) nln = nln - 1
  if (nln < 2) call abendm('too few lines')

  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
     '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6)')   &
     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
     '--- reduc. alt. ', xl, yl, mszx, mszy, mx, my
  call premsg(out)
  smx = float(mszx);  smy = float(mszy);  amxy = smx * smy
  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 (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 + nfx;  mym = my + nfy
  ixlm = ixl - mszx*nhx;  xlm = float(ixlm)/1000.
  iylm = iyl - mszy*nhy;  ylm = float(iylm)/1000.

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

  call gparma('Enter CXFUP matrix output filename ==> ', 81-ldr, fnam3(ldr:80))
  open(20,file=fnam3,form='unformatted',status='new')
  call gparma('Enter AXDEQ model output filename  ==> ', 81-ldr, fnam4(ldr:80))
  open(21,file=fnam4,status='new')
  call gparma('Enter AXOFF estim output filename  ==> ', 81-ldr, fnam5(ldr:80))
  open(22,file=fnam5,status='new')
  call clspin()
  call strdtm(sdtm)

  call dpcini('... search for Crossovers : ')
  npv = 0
  do na=1,nln
    sla = sl(na);  cln = sla(1:1);  kcl = 0
    if ((cln == 'B') .or. (cln == 'C') .or. (cln == 'X') .or.   &
        (cln == 'b') .or. (cln == 'c') .or. (cln == 'x')) kcl = 1
    do nb=1,nln
      call dpcent((na-1)*nln+(nb-1), nln*nln)
      slb = sl(nb);  cln = slb(1:1)
      if ((cln == 'B') .or. (cln == 'C') .or. (cln == 'X') .or.   &
          (cln == 'b') .or. (cln == 'c') .or. (cln == 'x')) then
        if (kcl == 1) cycle
      else
        if (kcl == 0) cycle
      endif
      kx = 0
      ka = 1;  if (na /= 1) ka = iq(na-1) + 1
      f1a2 = f(1,ka);  f2a2 = f(2,ka)
      do
        ka = ka + 1;  if (ka > iq(na)) exit
        f1a1 = f1a2;  f1a2 = f(1,ka)
        f2a1 = f2a2;  f2a2 = f(2,ka)
        kb = 1;  if (nb /= 1) kb = iq(nb-1) + 1
        f1b2 = f(1,kb);  f2b2 = f(2,kb)
        do
          kb = kb + 1;  if (kb > iq(nb)) exit
          f1b1 = f1b2;  f1b2 = f(1,kb)
          f2b1 = f2b2;  f2b2 = f(2,kb)
          call xcheck(f1a1,f2a1, f1a2,f2a2, f1b1,f2b1, f1b2,f2b2, kx)
          if (kx == 0) cycle
          npv = npv + 1
          if (npv > mcp) call abendm('[ERROR] too many control points')
          mvs = ka - 1;  kp = npv
          do
            if ((kp == 1) .or. (mv(kp-1) <= mvs)) exit
            kp = kp - 1
            mv(kp+1) = mv(kp);  sl1(kp+1) = sl1(kp);  sl2(kp+1) = sl2(kp)
          enddo
          mv(kp) = mvs;  sl1(kp) = sla;  sl2(kp) = slb
          if (kx == 1) exit
        enddo
        if (kx == 1) exit
      enddo
    enddo
  enddo
  call dpcent(nln*nln, nln*nln)

  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'EXDEQ3 : preparation for Gen.Mis-tie control <axdeq>'
  write(6,'(a)') '           to Control on all crossover points'
  write(6,'(a,a)') ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' CXFUP matrix outfile  : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' AXDEQ model outfile   : ', fnam4(1:lrtrim(fnam4))
  write(6,'(a,a)') ' AXOFF estim outfile   : ', fnam5(1:lrtrim(fnam5))
  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)')   &
    '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.'
  write(6,'(a,i5)') ' Number of Control points (for Mag. offset) :', npv

  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

  zd = vd
  do k=1,ndpt
    f(1,k) = (f(1,k)-xlm)*1000.;  ic = nint(f(1,k)/smx)
    f(2,k) = (f(2,k)-ylm)*1000.;  jc = nint(f(2,k)/smy)
    ht = tt(ic+1,jc+1)
    if (ht /= tnul) then
      z = f(3,k) - ht;  if (z < zd) zd = z
    endif
  enddo
  write(6,'(a,f6.0)') ' Min. height of Obs. above Eqv.Src.Surface :', zd
  if (zd <= 0.) call abendm('Invalid Obs. below Eqv.Src.')

  gn = amxy / (zd*zd*DPI)
  fg = zd*zd                            ! [ = (amxy / DPI)  /  gn ]

  call dpcini('... gen. CXFUP matrix data: ')
  write(20) area, KVER, ncc, iorg,korg, ispa,ispb,   &
          ixl, iyl, mszx, mszy, mx, my,   &
          nhx, nhy, npv, hw, vd, gn
  iln = 1;  sln = sl(iln);  iql = iq(iln)
  sla = '        ';  mva = 0
  nb = 1;  slb = sl1(nb);  mvb = mv(nb)
  do k=1,ndpt
    if (k > iql) then
      iln = iln + 1;  sln = sl(iln);  iql = iq(iln)
    endif
    do while (k > mvb)
      sla = slb;  mva = mvb;  nb = nb + 1
      if (nb <= npv) then
        slb = sl1(nb);  mvb = mv(nb)
      else
        slb = '        ';  mvb = ndpt
      endif
    enddo
    do n=1,npv;  c(n) = 0.;  enddo
    if (sln == sla) then
      if (sln == slb) then
        if (mva == mvb) then
          w = 0.
        else
          w = float(k-mva) / float(mvb-mva)
        endif
        c(nb-1) = 1.-w
        c(nb) = w
      else
        c(nb-1) = 1.
      endif
    else
      if (sln == slb) then
        c(nb) = 1.
      endif
    endif
    ic = nint(f(1,k)/smx);  jc = nint(f(2,k)/smy)
    do kj=-nhy,+nhy; do ki=-nhx,+nhx;  g(ki,kj) = 0.;   enddo; enddo
    do kj=-nhy,+nhy
      js = jc + kj
      if ((js < 0) .or. (js >= mym)) cycle
      y = f(2,k) - float(js)*smy
      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 = f(1,k) - float(is)*smx;  z = f(3,k) - ht
          rr = x*x + y*y + z*z;  r = sqrt(rr)
          g(ki,kj) = z*fg / (rr*r)      ! [ = z*amxy / (rr*r*DPI)  /  gn ]
        endif
      enddo
    enddo
    write(20) k, ic,jc, ((g(ki,kj),ki=-nhx,+nhx),kj=-nhy,+nhy), (c(n),n=1,npv)
    call dpcent(k, ndpt)
  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)

  write(22,'(a8,i6)') area, npv
  write(22,'(a)') fnam3(1:lrtrim(fnam3))
  do n=1,npv
    write(22,'(a8,i12,2x,a8,f10.2)') sl1(n), mv(n), sl2(n), 0.
  enddo
  close(22)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'EXDEQ3 : preparation for Gen.Mis-tie control <axdeq>'
    write(99,'(a)')   '           to Control on all crossover points'
    write(99,'(a,a)') ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Reduc.Alt.data infile : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' CXFUP matrix outfile  : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' AXDEQ model outfile   : ', fnam4(1:lrtrim(fnam4))
    write(99,'(a,a)') ' AXOFF estim outfile   : ', fnam5(1:lrtrim(fnam5))
    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)')   &
      '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.'
    write(99,'(a,i5)') ' Number of Control points (for Mag. offset) :', npv
    write(99,'(a,f6.0)') ' Min. height of Obs. above Eqv.Src.Surface :', zd
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end
!---
  subroutine xcheck(x1,y1, x2,y2, x3,y3, x4,y4, kx)
  implicit none
  real,intent(in) :: x1,y1, x2,y2, x3,y3, x4,y4
  integer,intent(out) :: kx
  real :: xa, xb, xc, ya, yb, yc, a, b, c
  kx = 0
  xa = x2 - x1;  xb = x4 - x3;  xc = x3 - x1
  ya = y2 - y1;  yb = y4 - y3;  yc = y3 - y1
  c = xa*yb - ya*xb
  if (c /= 0.) then
    b = (ya*xc - xa*yc) / c
    if ((0. <= b) .and. (b <= 1.)) then
      a = (yb*xc - xb*yc) / c
      if ((0. <= a) .and. (a <= 1.)) kx = 1
    endif
  endif
  return
  end subroutine