!--------( AAPTDP : point-dipole model fitting )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,parameter :: RAD = (180./3.141593)
  real,allocatable :: v(:,:), t(:,:), w(:,:), em(:,:)
  real :: c(3), sol(3), sol1(3), f(3)
  character(80) :: wdr, flog, fnam1, fnam3, fnam4, buf
  character(72) :: out;  character(25) :: strt;  character(20) :: sdtm
  character(8) :: area, area1, mname, rname, salt
  character(1) :: ps
  integer :: ldr, i, j, is, js, ni, nj, n, ky, kyr, k, IER
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my, ixt, iyt
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  real :: xs, ys, smx, smy, dip, dec, av, step, ah, at, the, phi
  real :: vnul, tnul, alt, xsw, ysw, xsz, ysz, xis, yjs, x, y, z
  real :: xsol, ysol, zsol, xsol1, ysol1, zsol1, drms, drms1, ww

  call premsg('AAPTDP : point-dipole model fitting')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = 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(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,*) 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.
  ixt = ixs + mszx*(mx-1);  iyt = iys + mszy*(my-1)
  smx = float(mszx)/1000.;  smy = float(mszy)/1000.

  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')
  allocate(v(mx,my));  allocate(t(mx,my));  allocate(w(mx,my))

  read(1,*) ((v(i,j),i=1,mx),j=1,my)
  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,a8)')   &
    '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)

  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('Header inconsistent')
    read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, tnul
    if ((ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx)   &
        .or. (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
      call abendm('Header inconsistent')
    read(1,*) ((t(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'
    tnul = 9999.9
    do j=1,my; do i=1,mx;  t(i,j) = alt;  enddo; enddo
  endif
  call premsg('  obs.surface : ' // strt)
  close(1)

  do j=1,my; do i=1,mx
    if (t(i,j) == tnul) then
      if (v(i,j) /= vnul) call abendm('input file broken')
      w(i,j) = vnul
    else
      w(i,j) = 0.
    endif
  enddo; enddo

  call premsg('')
  call premsg('Specify Direction Parameters (in degrees)')
  call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.)) call abendm('invalid Inclination')
  if ((dec < -180.) .or. (dec > 360.)) call abendm('invalid Declination')
  c(1) = cos(dip/RAD)*cos(dec/RAD)
  c(2) = cos(dip/RAD)*sin(dec/RAD)
  c(3) = sin(dip/RAD)

  call premsg('')
  call gparma('Enter Name-label for Model anom.  ==> ', 8, mname)
  call gparma('Enter Model anom. output filename ==> ', 81-ldr, fnam3(ldr:80))
  open(11,file=fnam3,status='new')
  write(11,'(a)') '# AAPTDP model anomaly'
  write(11,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
  write(11,'(a2,a8,4x,i4,4i8)') '# ', area, ncc, iorg, korg, ispa, ispb
  call premsg('')
  call gparma('Enter Name-label for Residual     ==> ', 8, rname)
  call gparma('Enter Residual output filename    ==> ', 81-ldr, fnam4(ldr:80))
  open(12,file=fnam4,status='new')
  write(12,'(a)') '# AAPTDP residuals'
  write(12,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
  write(12,'(a2,a8,4x,i4,4i8)') '# ', area, ncc, iorg, korg, ispa, ispb
  write(12,'(a,a)') '# Model : ', fnam3(1:lrtrim(fnam3))
  write(12,'(a2,a8,4x,i4)') '# ', mname, ncc
  call strdtm(sdtm)
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'AAPTDP : point-dipole model fitting'
  write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Model Anom. outfile : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' Residual 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)')   &
    'obs. mag.   ', xs,ys, mszx,mszy, mx,my, salt
  write(6,'(a,a25)') '  obs. surface  : ', strt
  write(6,'(a,f6.2,3x,a,f7.1)') '  ambient field :  dip=', dip, 'dec=', dec

  open(9,status='scratch')
  call premsg(' ')
  call premsg('Which way to define window,')
  call premsg('    UTM coordinates in km ("u") or')
  call premsg('    mesh-count within the area ("m"),')
  call gparma('  Select  "u" or "m" : ', 1, ps)
  if ((ps /= 'u') .and. (ps /= 'm')) call abendm('invalid selection')

  do
    if (ps == 'u') then
      call premsg ('Select fitting window (in km)')
      call gparmf2('    SW-corner coord. (Xsw, Ysw) [both 0 to End] ==> ',   &
                   xsw, ysw)
      if ((xsw == 0.) .and. (ysw == 0.)) exit
      is = nint((xsw-xs)/smx);  js = nint((ysw-ys)/smy)
      call gparmf2('    window size (N-S, E-W) ==> ', xsz, ysz)
      ni = nint(xsz/smx);  nj = nint(ysz/smy)
    else
      call premsg ('Select fitting window (in mesh-unit starting from 0)')
      call gparmi2('    window size (N-S, E-W) [both 0 to End] ==> ', ni, nj)
      if ((ni == 0) .and. (nj == 0)) exit
      call gparmi2('    SW-corner location   (IXs, JYw)   ==> ', is, js)
    endif
    xsz = smx*float(ni);  ysz = smy*float(nj)
    ni = ni + 1;  nj = nj + 1
    if ((is < 0) .or. (js < 0) .or. ((is+ni) > mx) .or. ((js+nj) > my)) then
      call premsg('  out of area'); cycle
    endif
    if ((ni <= 5) .or. (nj <= 5)) then
      call premsg('  too narrow window'); cycle
    endif
    allocate(em(4,ni*nj))

    xis = smx*float(is);  yjs = smy*float(js)
    xsw = xs + xis;  ysw = ys + yjs
    xsol = xsz/2.;  ysol = ysz/2.;  zsol = xsz/4.
    n = 0;  av = 0.
    do j=1,nj
      y = ysol - smy*float(j-1)
      do i=1,ni
        if (v(is+i,js+j) == vnul) cycle
        n = n + 1
        x = xsol - smx*float(i-1)
        z = zsol + t(is+i,js+j)/1000.
        call adipol(x, y, z, c, em(1,n))
        em(4,n) = v(is+i,js+j)
        av = av + em(4,n)
      enddo
    enddo
    if (n <= 12) call abendm('too few available data')
    av = av / float(n)
    do k=1,n;  em(4,k) = em(4,k) - av;  enddo
    call lsqsol(n, em, sol, drms)
    step = 0.1;  ky = 0;  kyr = 0
    do
      xsol1 = xsol;  ysol1 = ysol;  zsol1 = zsol
      if (ky == 0) then;  zsol1 = zsol + step
      else if (ky == 1) then;  zsol1 = zsol - step
      else if (ky == 2) then;  xsol1 = xsol + step
      else if (ky == 3) then;  xsol1 = xsol - step
      else if (ky == 4) then;  ysol1 = ysol + step
      else if (ky == 5) then;  ysol1 = ysol - step
      endif
      n = 0;  av = 0.
      do j=1,nj
        y = ysol1 - smy*float(j-1)
        do i=1,ni
          if (v(is+i,js+j) == vnul) cycle
          n = n + 1
          x = xsol1 - smx*float(i-1)
          z = zsol1 + t(is+i,js+j)/1000.
          call adipol(x, y, z, c, em(1,n))
          em(4,n) = v(is+i,js+j)
          av = av + em(4,n)
        enddo
      enddo
      if (n <= 12) call abendm('too few available data')
      av = av / float(n)
      do k=1,n;  em(4,k) = em(4,k) - av; enddo
      call lsqsol(n, em, sol1, drms1)
      if (drms1 < drms) then
        kyr = ky;  drms = drms1;  xsol = xsol1;  ysol = ysol1;  zsol = zsol1
        sol(1) = sol1(1);  sol(2) = sol1(2);  sol(3) = sol1(3);  cycle
      endif
      ky = ky + 1
      if (ky == 6) ky = 0
      if (ky /= kyr) cycle
      step = step / 2.;  ky = 0;  kyr = 0
      if (step <= 0.001) exit
    enddo
    ah = sqrt(sol(1)*sol(1) + sol(2)*sol(2))
    at = sqrt(sol(1)*sol(1) + sol(2)*sol(2) + sol(3)*sol(3))
    the = atan2(sol(3),ah)*RAD
    phi = atan2(sol(2),sol(1))*RAD
    write(6,'(a,2f9.3,2f8.3,a,3f8.3)')   &
      '#W# ', xsw,ysw, xsz,ysz, '  : ', xsol,ysol,zsol
    write(6,'(a,4x,3f10.3,2x,f10.3,2f8.2)')   &
      '##', (sol(k),k=1,3), at, the,phi
    write(11,'(a,2f9.3,2f8.3,a,3f8.3)')   &
      '#W# ', xsw,ysw, xsz,ysz, '  : ', xsol,ysol,zsol
    write(11,'(a,4x,3f10.3,2x,f10.3,2f8.2)')   &
      '##', (sol(k),k=1,3), at, the,phi
    write(9,*) xsw,ysw, xsz,ysz, xsol,ysol,zsol, (sol(k),k=1,3), at, the,phi
    do j=1,my
      y = ysol + yjs - smy*float(j-1)
      do i=1,mx
        if (t(i,j) == tnul) cycle
        x = xsol + xis - smx*float(i-1)
        z = zsol + t(i,j)/1000.
        call adipol(x, y, z, c, f)
        ww = sol(1)*f(1) + sol(2)*f(2) + sol(3)*f(3)
        w(i,j) = w(i,j) + ww
      enddo
    enddo
    deallocate(em)
  enddo

  call clspin()
  do j=1,my; do i=1,mx
    if (v(i,j) == vnul) cycle
    v(i,j) = v(i,j) - w(i,j)
    if (w(i,j) == vnul) v(i,j) = vnul
  enddo; enddo

  write(11,'(a8,4x,i4,4i8)') mname, ncc, iorg, korg, ispa, ispb
  write(11,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy, mx,my, vnul, alt
  do j=1,my
    write(11,'((f7.1,9(1x,f7.1)))') (w(i,j),i=1,mx)
  enddo
  if (alt == 0.) then
    write(11,'(a)') '# Var.Obs.Alt. data are copied from Source'
    write(11,'(a8,4x,i4,4i8)') mname, ncc, iorg, korg, ispa, ispb
    write(11,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy,mx,my, tnul, -1.
    do j=1,my
      write(11,'((f7.1,9(1x,f7.1)))') (t(i,j),i=1,mx)
    enddo
  endif
  close(11)
  endfile(9);  rewind(9)

  write(12,'(a8,4x,i4,4i8)') rname, ncc, iorg, korg, ispa, ispb
  write(12,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy,mx,my, vnul, alt
  do j=1,my
    write(12,'((f7.1,9(1x,f7.1)))') (v(i,j),i=1,mx)
  enddo
  if (alt == 0.) then
    write(12,'(a)') '# Var.Obs.Alt. data are copied from Source'
    write(12,'(a8,4x,i4,4i8)') rname, ncc, iorg, korg, ispa, ispb
    write(12,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy,mx,my, tnul, -1.
    do j=1,my
      write(12,'((f7.1,9(1x,f7.1)))') (t(i,j),i=1,mx)
    enddo
  endif
  close(12)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'AAPTDP : point-dipole model fitting'
    write(99,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Model Anom. outfile : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' Residual 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)')   &
      'obs. mag.   ', xs,ys, mszx,mszy, mx,my, salt
    write(99,'(a,a25)') '  obs. surface  : ', strt
    write(99,'(a,f6.2,3x,a,f7.1)') '  ambient field :  dip=', dip, 'dec=', dec
    do
      read(9,*,iostat=IER)   &
        xsw, ysw, xsz, ysz, xsol, ysol, zsol, (sol(k),k=1,3), at, the, phi
      if (IER < 0) exit
      write(99,'(a,2f9.3,2f8.3,a,3f8.3)')   &
        '#W# ', xsw,ysw, xsz,ysz, '  : ', xsol,ysol,zsol
      write(99,'(a,4x,3f10.3,2x,f10.3,2f8.2)')   &
        '##', (sol(k),k=1,3), at, the,phi
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end
!---
  subroutine lsqsol(n, em, sol, drms)
  implicit none
  integer,intent(in) :: n
  real,intent(in) :: em(4,n)
  real,intent(out) :: sol(3), drms
  real :: q(4,3)
  integer :: i, ii, k
  real :: g, sum

  do i=1,3; do ii=1,4;  q(ii,i) = 0.;  enddo; enddo
  do k=1,n; do i=1,3; do ii=i,4
    q(ii,i) = q(ii,i) + em(i,k)*em(ii,k)
  enddo; enddo; enddo
  q(1,2) = q(2,1);  q(1,3) = q(3,1);  q(2,3) = q(3,2)
  do k=1,3; do ii=k+1,4
    q(ii,k) = q(ii,k) / q(k,k)
    do i=1,3
      if (i /= k) q(ii,i) = q(ii,i) - q(ii,k)*q(k,i)
    enddo
  enddo; enddo
  sol(1) = q(4,1);  sol(2) = q(4,2);  sol(3) = q(4,3);  sum = 0.
  do k=1,n
    g = em(1,k)*q(4,1) + em(2,k)*q(4,2) + em(3,k)*q(4,3)
    sum = sum + (em(4,k)-g)*(em(4,k)-g)
  enddo
  drms = sqrt(sum/float(n))
  return
  end subroutine
!---
  subroutine adipol(x, y, z, c, f)
  implicit none
  real,intent(in) :: x, y, z, c(3)
  real,intent(out) :: f(3)
  real :: xx, yy, zz, xy, yz, zx, rr, r, r5

  xx = x * x;  yy = y * y;  zz = z * z
  xy = x * y;  yz = y * z;  zx = z * x
  rr = xx + yy + zz
  r = sqrt(rr)
  r5 = rr * rr * r / 100.
  f(1) = ((xx*3.-rr)*c(1) + xy*c(2)*3. + zx*c(3)*3.) / r5
  f(2) = ((yy*3.-rr)*c(2) + yz*c(3)*3. + xy*c(1)*3.) / r5
  f(3) = ((zz*3.-rr)*c(3) + zx*c(1)*3. + yz*c(2)*3.) / r5
  return
  end subroutine