!--------( AXDEQ : Gen.Mis-tie control <axdeq> )--------
  implicit none
  integer,parameter :: KVER = 2018
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:), r(:), q(:), s(:,:), p(:,:), v(:,:), t(:,:)
  real,allocatable :: g(:,:), c(:), d(:), e(:), b(:)
  integer,allocatable :: mv(:)
  character(8),allocatable :: sl(:), sl2(:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, sfnam, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1
  integer :: ldr, i, j, n, nhx, nhy, k, nver, IER
  integer :: ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixl, iyl, mszx, mszy, mx, my
  integer :: ixl1, iyl1, mszx1, mszy1, mxm1, mym1
  integer :: ixlm, iylm, mxm, mym, ic, jc
  integer :: mdpt, ndpt, npv, kc, kk, klp, l, nloop, nlp, npv1
  real :: tnul, xl, yl, vd, hw, ss, xlm, ylm, x, y, gn
  real :: a, anul, dmy, ds, dv, h, pqq, pvv, pvvp

  call premsg('AXDEQ : Gen.Mis-tie control <axdeq>')
  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 StdLIN data filename        ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  mdpt = 0; ndpt = 0
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) cycle
    if (buf(1:1) /= '#') mdpt = mdpt + 1
  enddo
  rewind(1)

  allocate(f(mdpt));  allocate(r(mdpt));  allocate(q(mdpt));

  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%') .or. (buf(1:1) == '#')) cycle
    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) = ' '
    i = index(buf,'nT');  if (i /= 0) buf(i:i+1) = '  '
    read(buf,*) x, y, h, f(ndpt)
  enddo
  close(1)

  call gparma('Enter CXFUP matrix data filename  ==> ', 81-ldr, fnam2(ldr:80))
  open(20,file=fnam2,form='unformatted',status='old')
  read(20) area, nver, ncc, iorg, korg, ispa, ispb,   &
         ixl, iyl, mszx, mszy, mx, my,   &
         nhx, nhy, npv, hw, vd, gn
  if (nver /= KVER) call abendm('Version mismatch')
  if (npv <= 0) call abendm('illegal number of control points')
  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size spec')
  if ((nhx < 10) .or. (nhy < 10)) call abendm('truncation too short')
  xl = float(ixl)/1000.;  ixlm = ixl - mszx*nhx;  xlm = float(ixlm)/1000.
  yl = float(iyl)/1000.;  iylm = iyl - mszy*nhy;  ylm = float(iylm)/1000.
  mxm = mx + nhx + nhx;  mym = my + nhy + nhy

  allocate(g(-nhx:+nhx,-nhy:+nhy))
  allocate(s(mxm,mym));  allocate(p(mxm,mym))
  allocate(v(mxm,mym));  allocate(t(mxm,mym))
  allocate(c(npv));  allocate(d(npv));  allocate(e(npv));  allocate(b(npv))
  allocate(mv(npv));  allocate(sl(npv));  allocate(sl2(npv))

  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,a8)')   &
     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
     '--- reduc. surf.', xl, yl,  mszx,mszy, mx, my
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
     '--- eq.src.anom.', xlm,ylm, mszx,mszy, mxm,mym, '    var.'
  call premsg(out)

  call gparma('Enter AXDEQ model in/out filename ==> ', 81-ldr, fnam3(ldr:80))
  open(21,file=fnam3,status='old')
  read(21,'(a)') buf
  do while (buf(1:1) == '#')
    read(21,'(a)') buf
  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
  read(21,'(a)') buf
  read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, anul, dmy, klp
  if ((area1 /= area) .or. (ncc1 /= ncc) .or.   &
      (iorg1 /= iorg) .or. (korg1 /= korg) .or.   &
      (ispa1 /= ispa) .or. (ispb1 /= ispb) .or.   &
      (ixl1 /= ixlm) .or. (iyl1 /= iylm) .or.   &
      (mszx1 /= mszx) .or. (mszy1 /= mszy) .or.   &
      (mxm1 /= mxm) .or. (mym1 /= mym))   &
    call abendm('AXDEQ Header inconsistent')
  read(21,*) ((s(i,j),i=1,mxm),j=1,mym)
  read(21,'(a)') buf
  do while (buf(1:1) == '#')
    read(21,'(a)') buf
  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
  read(21,'(a)') buf
  read(buf,*) ixl1, iyl1, mszx1, mszy1, mxm1, mym1, tnul
  if ((area1 /= area) .or. (ncc1 /= ncc) .or.   &
      (iorg1 /= iorg) .or. (korg1 /= korg) .or.   &
      (ispa1 /= ispa) .or. (ispb1 /= ispb) .or.   &
      (ixl1 /= ixlm) .or. (iyl1 /= iylm) .or.   &
      (mszx1 /= mszx) .or. (mszy1 /= mszy) .or.   &
      (mxm1 /= mxm) .or. (mym1 /= mym))   &
    call abendm('AXDEQ Header inconsistent')
  read(21,*) ((t(i,j),i=1,mxm),j=1,mym)
  rewind(21)

  call gparma('Enter AXOFF estim in/out filename ==> ', 81-ldr, fnam4(ldr:80))
  open(22,file=fnam4,status='old')
  read(22,'(a)') buf
  do while (buf(1:1) == '#')
    read(22,'(a)') buf
  enddo
  read(buf,'(a8,i6)') area1, npv1
  read(22,'(a)') sfnam
  if ((area1 /= area) .or. (npv1 /= npv) .or. (sfnam /= fnam2))   &
    call abendm('AXOFF file broken')
  do n=1,npv
    read(22,'(a8,i12,2x,a8,f10.2)') sl(n), mv(n), sl2(n), d(n)
  enddo
  rewind(22)

  call premsg('How many loops to execute ?')
  call gparmi('   (minus to restart from loop-0) ==> ', nloop)
  nlp = iabs(nloop)
  call clspin()
  call strdtm(sdtm)

  call dpcini('... read CXFUP matrix (0) : ')
  do kc=1,ndpt
    read(20) k, ic, jc, ((g(i,j),i=-nhx,+nhx),j=-nhy,+nhy), (c(kk),kk=1,npv)
    if (k /= kc) call abendm('CXFUP file broken')
    call dpcent(kc, ndpt)
  enddo
  read(20,iostat=IER) k
  if (IER >= 0) call abendm('CXFUP file broken')
  rewind(20)

  write(6,'(a)')   '----------------------------------------'
  write(6,'(a)')   'AXDEQ : Gen.Mis-tie control <axdeq>'
  write(6,'(a,a)') ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' CXFUP matrix infile   : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' AXDEQ model i/o file  : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' AXOFF estim i/o file  : ', 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)')   &
    '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
  write(6,'(a,i4)') ' Max. Loop count :', nloop

  if (nloop < 0) then
    klp = 0
    do i=1,mxm; do j=1,mym;  s(i,j) = 0.; p(i,j) = 0.;  enddo; enddo
    ds = 0.
  else
    ss = 0.
    do j=1,mym; do i=1,mxm
      if (t(i,j) /= tnul) ss = ss + s(i,j)*s(i,j)
      p(i,j) = 0.
    enddo; enddo
    do n=1,npv;  e(n) = 0.;  enddo
    ds = sqrt(ss / float(mxm*mym)) / gn
  endif

  call conv1(g,nhx,nhy, c,npv, s,mxm,mym, d, r,ndpt)
  ss = 0.
  do k=1,ndpt
    r(k) = f(k) - r(k);  ss = ss + r(k)*r(k)
  enddo
  dv = sqrt(ss / float(ndpt))
  write(6,'(a5,i4,3x,a10,f8.2,3x,a14,f7.2)')   &
    ' loop', klp, 'rms.Step :', ds, 'rms.Residual :', dv
  open(9,status='scratch')
  write(9,*) klp, ds, dv

  pvv = 1.
  do l=1,nlp
    klp = klp + 1
    pvvp = pvv
    call conv2(g,nhx,nhy, c,npv, r,ndpt, v,mxm,mym, b)
    pvv = 0.
    do j=1,mym; do i=1,mxm;  pvv = pvv + v(i,j)*v(i,j);  enddo; enddo
    do n=1,npv;  pvv = pvv + b(n)*b(n);  enddo
    a = pvv / pvvp
    do j=1,mym; do i=1,mxm;  p(i,j) = v(i,j) + a*p(i,j);  enddo; enddo
    do n=1,npv;  e(n) = b(n) + a*e(n);  enddo
    call conv1(g,nhx,nhy, c,npv, p,mxm,mym, e, q,ndpt)
    pqq = 0.
    do k=1,ndpt;  pqq = pqq + q(k)*q(k); enddo
    a = pvv / pqq;  ss = 0.
    do j=1,mym; do i=1,mxm
      s(i,j) = s(i,j) + a*p(i,j)
      if (t(i,j) /= tnul) ss = ss + p(i,j)*p(i,j)
    enddo; enddo
    ds = a * sqrt(ss / float(mxm*mym)) / gn
    do n=1,npv;  d(n) = d(n) + a*e(n);  enddo
    ss = 0.
    do k=1,ndpt
      r(k) = r(k) - a*q(k);  ss = ss + r(k)*r(k)
    enddo
    dv = sqrt(ss / float(ndpt))
    write(6,'(a5,i4,3x,a10,f8.2,3x,a14,f7.2)')   &
      ' loop', klp, 'rms.Step :', ds, 'rms.Residual :', dv
    write(9,*) klp, ds, dv
  enddo
  close(20)
  write(6,'(1x)')
  endfile(9)
  rewind(9)

  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., klp
  do j=1,mym
    do i=1,mxm
      if (s(i,j) < -9999.8) s(i,j) = -9999.8
      if (s(i,j) > +9999.8) s(i,j) = +9999.8
    enddo
    write(21,'((f7.1,9(1x,f7.1)))') (s(i,j),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)))') (t(i,j),i=1,mxm)
  enddo
  close(21)

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

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'AXDEQ : Gen.Mis-tie control <axdeq>'
    write(99,'(a,a)')   &
       ' StdLIN data infile    : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)')   &
       ' CXFUP matrix infile   : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)')   &
       ' AXDEQ model i/o file  : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)')   &
       ' AXOFF estim i/o file  : ', 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)')   &
      '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,i4)') ' Max. Loop count :', nloop
    do
      read(9,*,iostat=IER) klp, ds, dv
      if (IER < 0) exit
      write(99,'(a5,i4,3x,a10,f8.2,3x,a14,f7.2)')   &
        ' loop', klp, 'rms.Step :', ds, 'rms.Residual :', dv
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end
!---
  subroutine conv1(g,nhx,nhy, c,npv, s,mxm,mym, d, r,ndpt)
  implicit none
  integer,intent(in) :: nhx, nhy, npv, mxm, mym, ndpt
  real,intent(inout) :: g(-nhx:+nhx,-nhy:+nhy), c(npv)
  real,intent(in) :: s(mxm,mym), d(npv)
  real,intent(out) :: r(ndpt)
  character(8) :: area
  integer :: n, i, j, np, k, ic, jc, is, js
  real :: rp
  call dpcini('... read CXFUP matrix (1) : ')
  read(20) area, (i,n=1,18)
  do np=1,ndpt
    read(20) k, ic, jc, g, c
    rp = 0.
    do j=-nhy,+nhy
      js = jc + j
      if ((js <= 0) .or. (js > mym)) cycle
      do i=-nhx,+nhx
        is = ic + i
        if ((is <= 0) .or. (is > mxm)) cycle
        rp = rp + g(i,j)*s(is,js)
      enddo
    enddo
    do n=1,npv;  rp = rp + c(n)*d(n);  enddo
    r(np) = rp
    call dpcent(np, ndpt)
  enddo
  rewind(20)
  return
  end subroutine
!---
  subroutine conv2(g,nhx,nhy, c,npv, r,ndpt, s,mxm,mym, d)
  implicit none
  integer,intent(in) :: nhx, nhy, npv, ndpt, mxm, mym
  real,intent(inout) :: g(-nhx:+nhx,-nhy:+nhy), c(npv)
  real,intent(in) :: r(ndpt)
  real,intent(out) :: s(mxm,mym), d(npv)
  character(8) :: area
  integer :: n, i, j, np, k, ic, jc, is, js
  real  :: ss

  call dpcini('... read CXFUP matrix (2) : ')
  read(20) area, (i,n=1,18)
  do js=1,mym; do is=1,mxm;  s(is,js) = 0.;  enddo; enddo
  do n=1,npv;  d(n) = 0.;  enddo
  do np=1,ndpt
    read(20) k, ic, jc, g, c
    ss = 0.
    do j=-nhy,+nhy
      js = jc + j
      if ((js <= 0) .or. (js > mym)) cycle
      do i=-nhx,+nhx
        is = ic + i
        if ((is <= 0) .or. (is > mxm)) cycle
        s(is,js) = s(is,js) + g(i,j)*r(np)
      enddo
    enddo
    do n=1,npv;  d(n) = d(n) + c(n)*r(np);  enddo
    call dpcent(np, ndpt)
  enddo
  rewind(20)
  return
  end subroutine