c--------( TMCORR : terr.effect correction )--------
      parameter (LHW=100, LXO=1201, LYO=1201)
      parameter (LXT=LXO+LHW+LHW, LYT=LXO+LHW+LHW)
      parameter (LXTT=LXT*3+6, LYTT=LYT*3+6)
c
      dimension f(LXO,LYO), g(LXO,LYO), h(LXO,LYO)
      dimension fr(LXO,LYO), gr(LXO,LYO)
      dimension t(LXT,LYT), tt(LXTT,LYTT), vs(2), coef(3,2), e(6)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character*8 area, tarea, area1, salt
      character buf*80, strt*25, strs*20, out*72
c
      call premsg('TMCORR : terr.effect correction')
      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 Obs.Mag.Anom. data filename ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(12,file=fnam1,status='old')
    1 read(12,'(a)') buf
      if(buf(1:1).eq.'#') goto 1
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.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.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      xt = xs + float(mszx*(mx-1))/1000.
      yt = ys + float(mszy*(my-1))/1000.
      fmx = float(mszx)
      fmy = float(mszy)
      read(12,*) ((f(i,j),i=1,mx),j=1,my)
c
      call gparma('Enter Topography data filename    ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(10,file=fnam2,status='old')
    2 read(10,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      read(buf,'(a8,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
      nct = ncct
      if(ncct.ge.200) nct = ncct - 200
      if(nct.ge.1.and.nct.le.62) then
        iorgt = 0
        ispat = 0
        ispbt = 0
        korgt = (nct-30)*360 - 180
        if(nct.gt.60) korgt = 0
        if(nct.eq.61) iorgt = +5400
        if(nct.eq.62) iorgt = -5400
      endif
      read(10,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
      if(mxt.gt.LXT.or.myt.gt.LYT) call abendm('too large array size')
      xst = float(ixst)/1000.
      yst = float(iyst)/1000.
      xtr = xst + float(mszxt*(mxt-1))/1000.
      ytr = yst + float(mszyt*(myt-1))/1000.
      fmxt = float(mszxt)
      fmyt = float(mszyt)
      f3xt = fmxt / 3.
      f3yt = fmyt / 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.ne.ncc.or.iorgt.ne.iorg.or.korgt.ne.korg
     *            .or.ispat.ne.ispa.or.ispbt.ne.ispb) then
        write(out,'(11x,a8,6x,i4,4x,4i8)')
     *     tarea, ncct, iorgt,korgt, ispat,ispbt
        call premsg(out)
        call abendm('Header inconsistent')
      endif
      if(tarea.ne.area) call premsg('           ' // tarea)
      write(out,'(a,2a10,2x,2a8,2a6,a7)')
     *     '---   (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)')
     *     '--- topography  ', xst, yst, mszxt, mszyt, mxt, myt
      call premsg(out)
      if(xst.gt.xs.or.yst.gt.ys .or.
     *   xtr.lt.xt.or.ytr.lt.yt) call abendm('ranges conflict')
c
      if(alt.eq.0.) goto 3
      if(alt.lt.0.) then
        call gparmf('    Mag. Obs. Altitude [in m]  ==> ', alt)
        if (alt.le.0.) call abendm('invalid value')
        write(salt,'(1x,f7.0)') alt
      endif
      write(strt,'(f7.0,a)') alt, ' m above sea level'
      goto 4
    3 read(12,'(a)') buf
      if(buf(1:1).eq.'#') goto 3
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(ncc1.ge.200) nc1 = ncc1 - 200
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if(area1.ne.area.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header1 inconsistent')
      read(12,*) ix1, iy1, mszx1, mszy1, mx1, my1, hnul
      if(ix1.ne.ixs.or.iy1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my)
     *     call abendm('Header2 inconsistent')
      read(12,*) ((h(i,j),i=1,mx),j=1,my)
      strt = ' given as 2nd set data   '
    4 close(12)
      call premsg('  obs.surface : ' // strt)
c
      call premsg('Truncation of Source effect ?')
      call gparmf('    Half width (in km) :  ', hw)
      nhx = nint(hw*1000./fmx)
      nhy = nint(hw*1000./fmy)
      if(nhx.lt.10.or.nhy.lt.10) call abendm('truncation too short')
c
      call premsg('Specify Source-Bottom configuration')
      call prompt('    Constant Thickness (1) or Flat Bottom (0) ')
      call gparmi(': ', ks)
      if(ks.lt.0.or.ks.gt.1) call abendm('invalid value')
      strs = ' below sea level    '
      if(ks.eq.1) strs = ' below terr.surface '
      call gparmf('  Bottom-depth (in m)' // strs // ': ', botm)
c
      call premsg('Specify Direction Parameters (in degrees)')
      call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90. .or.
     *   dec.lt.-180..or.dec.ge.360.) call abendm('invalid value')
      call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
      if(the.lt.-90..or.the.gt.90. .or.
     *   phi.lt.-180..or.phi.ge.360.) call abendm('invalid value')
c
      call premsg('')
      call gparma('Enter Ter.M.Corr. output filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(13,file=fnam3,status='new')
      call premsg('Enter Auxiary(1) output filename')
      call gparma('     ( blank to suppress output ) ==> ',
     *            81-ldr, fnam4(ldr:80))
      call clspin()
      kout = 0
      if(fnam4(ldr:ldr).ne.' ') then
        open(8,file=fnam4,status='new')
        kout = 1
      endif
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'TMCORR : terr.effect correction'
      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))
      if(kout.eq.1) write(6,'(a,a)')
     *                 ' Auxiary(1) outfile : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs,  ys,  mszx,  mszy,  mx,  my, salt,
     *              xst, yst, mszxt, mszyt, mxt, myt
      write(6,2000) hw, nhx, nhy
      write(6,3000) strt, botm, strs, dip, dec, the, phi
c
      do 100 j=1,myt*3+6
      do 100 i=1,mxt*3+6
        tt(i,j) = tnul
  100 continue
      do 110 i=1,mxt
        li = i*3 + 2
        tn = t(i,1)
        tt(li,5) = tn
        do 110 j=2,myt
          tp = tn
          lj = j*3 + 2
          tn = t(i,j)
          tt(li,lj) = tn
          if(tp.eq.tnul.or.tn.eq.tnul) goto 110
          tt(li,lj-2) = (tp + tp + tn) / 3.
          tt(li,lj-1) = (tp + tn + tn) / 3.
  110 continue
      do 120 lj=5,myt*3+2
        tn = tt(5,lj)
        do 120 i=2,mxt
          tp = tn
          li = i*3 + 2
          tn = tt(li,lj)
          if(tp.eq.tnul.or.tn.eq.tnul) goto 120
          tt(li-2,lj) = (tp + tp + tn) / 3.
          tt(li-1,lj) = (tp + tn + tn) / 3.
  120 continue
      call prompt('... calc. unit mag. effect: ')
      call dpcent(0, 0)
      xsm = (xs-xst) * 1000.
      ysm = (ys-yst) * 1000.
      rad = 3.141593 / 180.
      ex = cos(dip*rad) * cos(dec*rad)
      ey = cos(dip*rad) * sin(dec*rad)
      ez = sin(dip*rad)
      px = cos(the*rad) * cos(phi*rad)
      py = cos(the*rad) * sin(phi*rad)
      pz = sin(the*rad)
      cx = ex * px
      cy = ey * py
      cz = ez * pz
      cxy = (ex*py + ey*px) / 2.
      cyz = (ey*pz + ez*py) / 2.
      czx = (ez*px + ex*pz) / 2.
      c3 = cx - cz
      c4 = cy - cz
      c5 = cx - cy
      c1 = c3 + c5
      c2 = c4 - c5
      am = fmxt*fmyt
      do 200 j=1,my
        yc = ysm + fmy*float(j-1)
        jc = nint(yc/fmyt) + 1
        do 200 i=1,mx
          xc = xsm + fmx*float(i-1)
          ic = nint(xc/fmxt) + 1
          if(alt.ne.0.) h(i,j) = alt
          if(h(i,j).eq.hnul) goto 7
          hc = h(i,j)
          if(hc.lt.(t(ic,jc)+1.)) goto 7
          fv = 0.
          do 210 kj=-nhy,nhy
            jt = jc + kj
            if(jt.le.0.or.jt.gt.myt) goto 210
            y = float(jt-1)*fmyt - yc
            do 220 ki=-nhx,nhx
              it = ic + ki
              if(it.le.0.or.it.gt.mxt) goto 220
              x = float(it-1)*fmxt - xc
              if(t(it,jt).eq.tnul) goto 220
              if(ki.eq.0.and.kj.eq.0) goto 5
              if(iabs(ki).le.1.and.iabs(kj).le.1) goto 220
              za = hc - t(it,jt)
              zb = hc + botm
              if(ks.eq.1) zb = za + botm
              if(abs(za).lt.1.) za = 1.
              if(abs(zb).lt.1.) zb = 1.
              if(iabs(ki).le.5.and.iabs(kj).le.5) then
                call prism(x, y, fmxt, fmyt, za, zb, e)
                v = cyz*e(1) + czx*e(2) + cxy*e(3)
     *              + cx*e(4) + cy*e(5) + cz*e(6)
              else
                zza = za * za
                zzb = zb * zb
                xx = x * x
                xy = x * y
                yy = y * y
                ss = xx + yy
                ssss = ss * ss
                rra = ss + zza
                rrb = ss + zzb
                ra = sqrt(rra)
                rb = sqrt(rrb)
                dza = za / ra
                dzb = zb / rb
                u1a = dza / ssss
                u1b = dzb / ssss
                u2a = u1a * zza / rra
                u2b = u1b * zzb / rrb
                u3a = 2. / ra / rra
                u3b = 2. / rb / rrb
                v1 = c1*xx + c2*yy + cxy*xy*6.
                v2 = c3*xx + c4*yy + cxy*xy*2.
                v3 = czx*x + cyz*y
                v = ((u1b-u1a)*v1 - (u2b-u2a)*v2 - (u3b-u3a)*v3) * am
              endif
              goto 6
    5         v = 0.
              li = it*3 + 2
              lj = jt*3 + 2
              do 230 ljs=-4,4
                yp = float(ljs) * f3yt
                do 230 lis=-4,4
                  xp = float(lis) * f3xt
                  tp = tt(li+lis,lj+ljs)
                  if(tp.eq.tnul) goto 230
                  zat = hc - tp
                  zbt = hc + botm
                  if(ks.eq.1) zbt = zat + botm
                  if(abs(zat).lt.1.) zat = 1.
                  if(abs(zbt).lt.1.) zbt = 1.
                  call prism(xp, yp, f3xt, f3yt, zat, zbt, e)
                  w = cyz*e(1) + czx*e(2) + cxy*e(3)
     *                + cx*e(4) + cy*e(5) + cz*e(6)
                  v = v + w
  230         continue
    6         fv = fv + v
  220       continue
  210     continue
          g(i,j) = fv * 100.
          goto 8
    7     if(f(i,j).ne.vnul) call abendm('obs. alt. undefined')
          g(i,j) = vnul
    8     call dpcent((j-1)*mx+i, mx*my)
  200 continue
c
      if(kout.eq.1) then
        write(8,'(a8,i4,4x,4i8)') 'Tef-1A/m', ncc, iorg,korg,ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 810 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mx)
  810   continue
        write(6,'(a)') 'AUX1 <== Tef-1A/m : Ter.Efct of Unit magnetiz.'
      endif
c
      ft = 0.
      gt = 0.
      n = 0
      do 300 j=1,my
      do 300 i=1,mx
        fr(i,j) = vnul
        gr(i,j) = vnul
        if(f(i,j).eq.vnul) goto 300
        ft = ft + f(i,j)
        gt = gt + g(i,j)
        n = n + 1
  300 continue
      an = float(n)
      ft = ft / an
      gt = gt / an
c
      call sm2opn(1, 2)
      cmx = float(mx+1) / 2.
      cmy = float(my+1) / 2.
      do 400 j=1,my
        y = float(j) - cmy
        do 400 i=1,mx
          x = float(i) - cmx
          if(f(i,j).eq.vnul) goto 400
          vs(1) = f(i,j)-ft
          vs(2) = g(i,j)-gt
          call sm2ex(x, y, vs)
  400 continue
      call sm2cls(coef, 3, 2)
      write(6,4000)
      write(6,5000) 'observation', coef(1,1)+ft, coef(2,1), coef(3,1)
      write(6,5000) '1A/m effect', coef(1,2)+gt, coef(2,2), coef(3,2)
c
      fg = 0.
      gg = 0.
      do 500 j=1,my
        y = float(j) - cmy
        do 500 i=1,mx
          x = float(i) - cmx
          call sm2rv(x, y, vs)
          if(g(i,j).ne.vnul) gr(i,j) = g(i,j) - (vs(2) + gt)
          if(f(i,j).eq.vnul) goto 500
          fr(i,j) = f(i,j) - (vs(1) + ft)
          fg = fg + fr(i,j)*gr(i,j)
          gg = gg + gr(i,j)*gr(i,j)
  500 continue
      gm = fg / gg
      write(6,6000) gm
c
      do 600 j=1,my
      do 600 i=1,mx
        if(g(i,j).ne.vnul) then
          gr(i,j) = gr(i,j) * gm
          g(i,j) = g(i,j) * gm
        endif
        if(f(i,j).ne.vnul) f(i,j) = f(i,j) - g(i,j)
  600 continue
      write(13,'(a)') '# terr.effect corrected magnetic anomalies'
      write(13,'(a8,i4,4x,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 700 j=1,my
        write(13,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  700 continue
      write(13,'(a8,i4,4x,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 710 j=1,my
        write(13,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  710 continue
      close(13)
c
      if(kout.eq.1) then
        write(8,'(a8,i4,4x,4i8)') 'Obs-Trnd', ncc, iorg,korg,ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 820 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (fr(i,j),i=1,mx)
  820   continue
        write(6,'(a)') 'AUX1 <== Obs-Trnd : Trend-removed Observation'
        write(8,'(a8,i4,4x,4i8)') 'TerEfect', ncc, iorg,korg,ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 830 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mx)
  830   continue
        write(6,'(a)') 'AUX1 <== TerEfect : Simulated Terrain effect'
        write(8,'(a8,i4,4x,4i8)') 'Tef-Trnd', ncc, iorg,korg,ispa,ispb
        write(8,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, vnul, -1.
        do 840 j=1,my
          write(8,'((f7.1,9(1x,f7.1)))') (gr(i,j),i=1,mx)
  840   continue
        write(6,'(a)') 'AUX1 <== Tef-Trnd : Trend-removed Terr.effect'
        close(8)
      endif
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        write(99,'(//a)') 'TMCORR : terr.effect correction'
        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))
        if(kout.eq.1) write(99,'(a,a)')
     *       ' Auxiary(1) outfile : ', fnam4(1:lrtrim(fnam4))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xs,  ys,  mszx,  mszy,  mx,  my, salt,
     *                 xst, yst, mszxt, mszyt, mxt, myt
        write(99,2000) hw, nhx, nhy
        write(99,3000) strt, botm, strs, dip, dec, the, phi
        write(99,4000)
        write(99,5000) 'observation', coef(1,1)+ft, coef(2,1), coef(3,1)
        write(99,5000) '1A/m effect', coef(1,2)+gt, coef(2,2), coef(3,2)
        write(99,6000) gm
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'obs. mag.   ', 2f10.2, 2i8, 2i6, a8/
     *       3x, 'topography  ', 2f10.2, 2i8, 2i6)
 2000 format('  Half width of calc.window :', f6.2, ' km (', i4,
     *       ' *', i4,'  terr.mesh )')
 3000 format('  obs. surface  : ', a25/
     *       '  bottom depth  : ', f8.1, ' m', a20/
     *       '  ambient field :  dip=', f6.1, 3x, 'dec=', f7.1/
     *       '  magnetization :  dip=', f6.1, 3x, 'dec=', f7.1)
 4000 format(' Trend surface in nT'/ 5x, '( x/y : N/E coord. ',
     *       'in mesh unit referred to the center of area )')
 5000 format('  for ', a11, ' : ', f8.2, '  +  ',
     *       f7.3, ' * x  +  ', f7.3, ' * y')
 6000 format(' Magnetization of terrain model : ', f6.2, ' A/m')
      end
c---
      subroutine prism(xc, yc, fmx, fmy, za, zb, e)
      dimension e(6)
      do 1 i=1,6
        e(i) = 0.
    1 continue
      xa = xc - fmx/2.
      xb = xc + fmx/2.
      ya = yc - fmy/2.
      yb = yc + fmy/2.
      do 2 l=0,7
        x = xb
        y = yb
        z = zb
        if(mod(l,2).ne.0) x = xa
        if(mod(l/2,2).ne.0) y = ya
        if(l.ge.4) z = za
        rr = x*x + y*y + z*z
        r = sqrt(rr)
        e(1) = -e(1) + alog((r-x)/(r+x))
        e(2) = -e(2) + alog((r-y)/(r+y))
        e(3) = -e(3) - alog(r+z)
        e(4) = -e(4) - atan2((x*y),(rr+r*z-y*y))
        e(5) = -e(5) - atan2((x*y),(rr+r*z-x*x))
        e(6) = -e(6) + atan((x*y)/(r*z))
        if(l.ne.1.and.l.ne.5) goto 2
        do 3 i=1,6
          e(i) = -e(i)
    3   continue
    2 continue
      return
      end
