!--------( CXDEQ : reduction of anom. by Gen.Mis-tie control <axdeq> )--------
  implicit none
  integer,parameter :: KVER = 2018
  real,parameter :: DPI = 3.141593 * 2.
  integer :: lwkdir, lrtrim
  real,allocatable :: s(:,:), t(:,:), f(:,:), h(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1
  integer :: ldr, i, j, nhx, nhy, nhx1, nhy1, nsx, nsy, nver
  integer :: nc, 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, klp, npv1
  integer :: is, js, ixh, ixhm, iyh, iyhm, ki, kj, mx1, my1
  real :: tnul, xl, yl, hw, ss, xlm, ylm, x, y, gn
  real :: a, anul, dmy, amxygn, hnul, hw1, hwm, rr, smx, smy, sr, vd1, z

  call premsg('CXDEQ : reduction of aanom. by 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 Reduction Alt. filename    ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  read(1,'(a)') buf
  do while (buf(1:1) == '#')
    read(1,'(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(1,'(a)') buf
  read(buf,*) ixl, iyl, mszx, mszy, mx, my, hnul
  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')

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

  read(1,*) ((h(i,j),i=1,mx),j=1,my)
  close(1)
  xl = float(ixl)/1000.;  ixh = ixl + mszx*(mx-1);  smx = float(mszx)
  yl = float(iyl)/1000.;  iyh = iyl + mszy*(my-1);  smy = float(mszy)

  call gparma('Enter CXFUP matrix data filename ==> ', 81-ldr, fnam2(ldr:80))
  open(20,file=fnam2,form='unformatted',status='old')
  read(20) area1, nver, ncc1, iorg1, korg1, ispa1, ispb1,   &
         ixl1, iyl1, mszx1, mszy1, mx1, my1,   &
         nhx1, nhy1, npv1, hw1, vd1, gn
  if (nver /= KVER) call abendm('Version mismatch')
  close(20)
  nc1 = ncc1
  if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
  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) .or.   &
      (ixl1 /= ixl) .or. (iyl1 /= iyl) .or.   &
      (mszx1 /= mszx) .or. (mszy1 /= mszy) .or.   &
      (mx1 /= mx) .or. (my1 /= my))   &
    call abendm('CXFUP Header inconsistent')
  amxygn = smx * smy / gn

  call gparma('Enter AXDEQ model input filename ==> ', 81-ldr, fnam3(ldr:80))
  open(2,file=fnam3,status='old')
  read(2,'(a)') buf
  do while (buf(1:1) == '#')
    read(2,'(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(2,'(a)') buf
  read(buf,*) ixlm, iylm, mszx1, mszy1, mxm, mym, anul, dmy, klp
  if ((area1 /= area) .or. (ncc1 /= ncc) .or.   &
      (iorg1 /= iorg) .or. (korg1 /= korg) .or.   &
      (ispa1 /= ispa) .or. (ispb1 /= ispb) .or.   &
      (mszx1 /= mszx) .or. (mszy1 /= mszy))   &
    call abendm('AXDEQ Header inconsistent')
  if ((mxm < 10) .or. (mym < 10)) call abendm('truncation too short')

  allocate(s(mxm,mym));  allocate(t(mxm,mym))

  read(2,*) ((s(i,j),i=1,mxm),j=1,mym)
  read(2,'(a)') buf
  do while (buf(1:1) == '#')
    read(2,'(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(2,'(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(2,*) ((t(i,j),i=1,mxm),j=1,mym)
  close(2)
  xlm = float(ixlm)/1000.;  ixhm = ixlm + mszx*(mxm-1)
  ylm = float(iylm)/1000.;  iyhm = iylm + mszy*(mym-1)
  if ((ixl <= ixlm) .or. (ixh >= ixhm) .or.   &
      (iyl <= iylm) .or. (iyh >= iyhm))   &
    call abendm('ranges conflict')
  nsx = (ixl-ixlm) / mszx;  nsy = (iyl-iylm) / mszy

  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. alt. ', 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 premsg('  Truncation of Source effect ?')
  call gparmf('    Half width (in km) : ', hw)
  hwm = hw * 1000.;  nhx = nint(hwm/smx);  nhy = nint(hwm/smy)

  call gparma('Enter Reduced Anom. out filename ==> ', 81-ldr, fnam4(ldr:80))
  call clspin()
  call strdtm(sdtm)
  open(3,file=fnam4,status='new')
  write(6,'(a)')   '----------------------------------------'
  write(6,'(a)')   'CXDEQ : reduction of anom. by Gen.Mis-tie control <axdeq>'
  write(6,'(a,a)') ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' CXFUP matrix infile   : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' AXDEQ model infile    : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' Reduced anom. 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)')   &
    '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 )'

  call dpcini('... calc. Reduced anomaly : ')
  do j=1,my; do i=1,mx
    if (h(i,j) == hnul) then
      f(i,j) = 9999.9
    else
      ss = 0.
      do kj=-nhy,+nhy
        js = j + nsy + kj;  y = float(kj)*smy
        do ki=-nhx,+nhx
          is = i + nsx + ki;  x = float(ki)*smx
          if ((is <= 0) .or. (is > mxm)) cycle
          if ((js <= 0) .or. (js > mym)) cycle
          if (t(is,js) == tnul) cycle
          z = h(i,j) - t(is,js);  rr = x*x + y*y + z*z;  sr = sqrt(rr)
          a = z*amxygn / (rr*sr*DPI);  ss = ss + a*s(is,js)
        enddo
      enddo
      f(i,j) = ss
    endif
    call dpcent((j-1)*mx+i, mx*my)
  enddo; enddo

  write(3,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
     ixl, iyl, mszx, mszy, mx, my, 9999.9, 0.
  do j=1,my
    write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo
  write(3,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
     ixl, iyl, mszx, mszy, mx, my, hnul, -1.
  do j=1,my
    write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  enddo
  close(3)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)')   &
      'CXDEQ : reduction of anom. by Gen.Mis-tie control <axdeq>'
    write(99,'(a,a)') ' Reduction alt. infile : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' CXFUP matrix infile   : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' AXDEQ model infile    : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' Reduced anom. 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)')   &
      '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 )'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end