!--------( TMCFIX : correct for fixed terr.magnetization )--------
!--------(               [do not fill deep sea with mass])--------
  implicit none
  integer :: lwkdir, lrtrim;  real :: calma
  real,allocatable :: f(:,:), g(:,:), h(:,:), t(:,:), tt(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, buf
  character(72) :: out;  character(25) :: strt
  character(20) :: strs, sdtm
  character(8) :: area, tarea, area1, salt
  integer :: ldr, i, j, nhx, nhy, ks, li, lj, ic, jc
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: nct, ncct, iorgt, korgt, ispat, ispbt
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ix1, iy1, mszx1, mszy1, mx1, my1
  integer :: ixst, iyst, mszxt, mszyt, mxt, myt
  integer :: ki, kj, it, jt, lis, ljs
  real :: alt, vnul, xs, ys, xt, yt, smx, smy, hnul, tnul
  real :: xst, yst, xtr, ytr, smxt, smyt, hw, botm, x, y
  real :: dip, dec, the, phi, gm, tn, tp, xsm, ysm, am, xc, yc, hc
  real :: xp, yp, f3xt, f3yt, zap, szp, fv, za, sz

  call premsg(' TMCFIX : correct for fixed terr.magnetization')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr;  fnam3 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Obs.Mag.Anom. data filename ==> ', 81-ldr, fnam1(ldr:80))
  open(12,file=fnam1,status='old')
  do
    read(12,'(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(12,'(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'
  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')

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

  xs = float(ixs)/1000.;  ys = float(iys)/1000.
  xt = xs + float(mszx*(mx-1))/1000.
  yt = ys + float(mszy*(my-1))/1000.
  smx = float(mszx);  smy = float(mszy)
  read(12,*) ((f(i,j),i=1,mx),j=1,my)

  if (alt == 0.) then
    do
      read(12,'(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(12,*) ix1, iy1, mszx1, mszy1, mx1, my1, hnul
    if ((ix1 /= ixs) .or. (iy1 /= iys) .or. (mszx1 /= mszx)   &
        .or. (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
      call abendm('Header2 inconsistent')
    read(12,*) ((h(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'
  endif
  close(12)

  call gparma('Enter Topography data filename    ==> ', 81-ldr, fnam2(ldr:80))
  open(10,file=fnam2,status='old')
  do
    read(10,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') tarea, ncct, iorg,korgt, ispat,ispbt
    if ((ncct<0) .or. (ncct>=400)) call abendm('invalid coord.no.')
    nct = ncct
    if (ncct < 200) then
      nct = ncct; ncct = nct + 800
    else
      ncct = ncct - 200; nct = ncct
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    nct = ncct
    if ((ncct >= 800) .and. (ncct < 1000)) nct = ncct - 800
  endif
  if ((nct >= 1) .and. (nct <= 62)) then
    iorgt = 0;  ispat = 0;  ispbt = 0
    korgt = (nct-30)*360 - 180
    if (nct > 60) korgt = 0
    if (nct == 61) iorgt = +5400
    if (nct == 62) iorgt = -5400
  endif
  read(10,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul

  allocate(t(mxt,myt));  allocate(tt(mxt*3+6,myt*3+6))

  xst = float(ixst)/1000.
  yst = float(iyst)/1000.
  xtr = xst + float(mszxt*(mxt-1))/1000.
  ytr = yst + float(mszyt*(myt-1))/1000.
  smxt = float(mszxt);  f3xt = smxt / 3.
  smyt = float(mszyt);  f3yt = smyt / 3.
  read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
  close(10)
  write(out,'(a,a8,a,i4,4x,4i8)')   &
    '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  if ((ncct /= ncc) .or. (iorgt /= iorg) .or. (korgt /= korg)   &
      .or. (ispat /= ispa) .or. (ispbt /= ispb)) then
    write(out,'(11x,a8,6x,4x,i4,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
    call premsg(out)
    call abendm('Header inconsistent')
  endif
  if (tarea /= area) call premsg('           ' // tarea)
  write(out,'(a,2a10,2x,2a8,2a6,a8)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a,2f10.2,2i8,2i6,a8)')   &
    '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  write(out,'(a,2f10.2,2i8,2i6,a8)')   &
    '--- topography  ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
  call premsg(out)
  if ((xst > xs) .or. (yst > ys) .or.   &
      (xtr < xt) .or. (ytr < yt)) call abendm('ranges conflict')
  call premsg('  obs.surface : ' // strt)

  call premsg('Truncation of Source effect ?')
  call gparmf('    Half width of calc.window (in km) :  ', hw)
  nhx = nint(hw*1000./smx)
  nhy = nint(hw*1000./smy)
  if ((nhx < 10) .or. (nhy < 10)) call abendm('truncation too short')

  call premsg('Specify Source-Bottom configuration')
  call gparmi('    Constant Thickness (1) or Flat Bottom (0) : ', ks)
  if ((ks < 0) .or. (ks > 1)) call abendm('invalid value')
  strs = ' below sea level    '
  if (ks == 1) strs = ' below terr.surface '
  call gparmf('  Bottom-depth (in m)' // strs // ': ', botm)

  call premsg('Specify Direction Parameters (in degrees)')
  call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.) .or.   &
      (dec < -180.) .or. (dec >= 360.)) call abendm('invalid value')
  call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
  if ((the < -90.) .or. (the > 90.) .or.   &
      (phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')
  call gparmf('Specify Terrain Magnetization (A/m) ==> ', gm)

  call premsg('')
  call gparma('Enter Ter.M.Corr. output filename ==> ', 81-ldr, fnam3(ldr:80))
  open(13,file=fnam3,status='new')
  call clspin()
  call strdtm(sdtm)
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'TMCFIX : correct for fixed terr.magnetization'
  write(6,'(a,a)') ' Obs. Mag. infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Topography infile  : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' T.M.Corr. outfile  : ', fnam3(1:lrtrim(fnam3))
  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,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'topography  ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
  write(6,'(3x,a,f6.2,a,i4,a,i4,a)') 'Half width of calc.window :',   &
    hw, ' km (', nhx, ' *', nhy, ' terr.mesh )'
  write(6,'(3x,a,a25)')        'obs. surface  : ', strt
  write(6,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'ambient field : ', 'dip=', dip, 'dec=', dec
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'magnetization : ', 'dip=', the, 'dec=', phi
  write(6,'(2x,a,f6.2,a)') 'Terrain magnetization : ', gm, ' A/m'

  do j=1,myt*3+6; do i=1,mxt*3+6
    tt(i,j) = tnul
  enddo; enddo
  do i=1,mxt
    li = i*3 + 2;  tn = t(i,1);  tt(li,5) = tn
    do j=2,myt
      tp = tn;  lj = j*3 + 2;  tn = t(i,j);  tt(li,lj) = tn
      if ((tp == tnul) .or. (tn == tnul)) cycle
      tt(li,lj-2) = (tp + tp + tn) / 3.
      tt(li,lj-1) = (tp + tn + tn) / 3.
    enddo
  enddo
  do lj=5,myt*3+2
    tn = tt(5,lj)
    do i=2,mxt
      tp = tn;  li = i*3 + 2;  tn = tt(li,lj)
      if ((tp == tnul) .or. (tn == tnul)) cycle
      tt(li-2,lj) = (tp + tp + tn) / 3.
      tt(li-1,lj) = (tp + tn + tn) / 3.
    enddo
  enddo
  call dpcini('... calc. unit mag. effect: ')
  xsm = (xs-xst) * 1000.
  ysm = (ys-yst) * 1000.
  call magafd(the, phi, dip, dec)
  am = smxt*smyt
  do j=1,my
    yc = ysm + smy*float(j-1)
    jc = nint(yc/smyt) + 1
    do i=1,mx
      call dpcent((j-1)*mx+i-1, mx*my)
      xc = xsm + smx*float(i-1)
      ic = nint(xc/smxt) + 1
      if (alt /= 0.) h(i,j) = alt
      if ((h(i,j) == hnul) .or. (h(i,j) < (t(ic,jc)+1.))) then
        if (f(i,j) /= vnul) call abendm('obs. alt. undefined')
        g(i,j) = vnul; cycle
      endif
      hc = h(i,j);  fv = 0.
      do kj=-nhy,nhy
        jt = jc + kj
        if ((jt <= 0) .or. (jt > myt)) cycle
        y = float(jt-1)*smyt - yc
        do ki=-nhx,nhx
          it = ic + ki
          if ((it <= 0) .or. (it > mxt)) cycle
          x = float(it-1)*smxt - xc
          if (t(it,jt) == tnul) cycle
          if ((ki == 0) .and. (kj == 0)) then
            li = ic*3 + 2;  lj = jc*3 + 2
            do ljs=-4,4
              yp = y + float(ljs)*f3yt
              do lis=-4,4
                xp = x + float(lis)*f3xt
                tp = tt(li+lis,lj+ljs)
                if (tp == tnul) cycle
                zap = hc - tp;  szp = botm
                if (ks == 0) szp = botm + tp
                if (szp > 0.) then
                  call mprism(1., xp, f3xt, yp, f3yt, zap, szp)
                  fv = fv + calma(0., 0., 0.)
                endif
              enddo
            enddo
          else
            if ((iabs(ki) <= 1) .and. (iabs(kj) <= 1)) cycle
            za = hc - t(it,jt);  sz = botm
            if (ks == 0) sz = botm + t(it,jt)
            if (sz > 0.) then
              if ((iabs(ki) <= 5) .and. (iabs(kj) <= 5)) then
                call mprism(1., x, smxt, y, smyt, za, sz)
              else
                call mvline(1., x, y, za, sz, am)
              endif
              fv = fv + calma(0., 0., 0.)
            endif
          endif
        enddo
      enddo
      g(i,j) = fv
    enddo
  enddo
  call dpcent(mx*my, mx*my)

  do j=1,my; do i=1,mx
    if (f(i,j) /= vnul) f(i,j) = f(i,j) - g(i,j)*gm
  enddo; enddo
  write(13,'(a)') '# terr.effect corrected magnetic anomalies'
  write(13,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(13,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy, mx,my, vnul, alt
  do j=1,my
    write(13,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo
  write(13,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(13,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixs,iys, mszx,mszy, mx,my, hnul, -1.
  do j=1,my
    write(13,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  enddo
  close(13)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)')   &
      'TMCFIX : correct for fixed terr.magnetization'
    write(99,'(a,a)')   &
      ' Obs. Mag. infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)')   &
      ' Topography infile  : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)')   &
      ' T.M.Corr. outfile  : ', fnam3(1:lrtrim(fnam3))
    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,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'topography  ', xst, yst, mszxt, mszyt, mxt, myt, '    --- '
    write(99,'(3x,a,f6.2,a,i4,a,i4,a)') 'Half width of calc.window :',   &
      hw, ' km (', nhx, ' *', nhy, '  terr.mesh )'
    write(99,'(3x,a,a25)')        'obs. surface  : ', strt
    write(99,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'ambient field : ', 'dip=', dip, 'dec=', dec
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'magnetization : ', 'dip=', the, 'dec=', phi
    write(99,'(2x,a,f6.2,a)') 'Terrain magnetization : ', gm, ' A/m'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end