c--------( LCECORR : loop-current effect correction )--------
c
      parameter (LXO=1201, LYO=1201, NPT=150)
      parameter (RAD=180./3.14159)
      parameter (LBUF=255)
      character*256 sbf
c
      dimension f(LXO,LYO), g(LXO,LYO), h(LXO,LYO)
      dimension fr(LXO,LYO), gr(LXO,LYO)
      dimension xpt(NPT+1), ypt(NPT+1), dpt(NPT+1)
      dimension vs(2), coef(3,2)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character area*8, area1*8, strt*25, pat*8, salt*8
      character buf*80, out*72
c
      sbf = ' '
      sbf(LBUF+1:LBUF+1) = char(0)
      call premsg('LCECORR : loop-current 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.
      read(12,*) ((f(i,j),i=1,mx),j=1,my)
      if(alt.lt.0.) then
        call gparmf('  Altitude [in m] for this data ==> ', 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'
c
      call gparma('Enter Loop-Location data filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(10,file=fnam2,status='old')
    2 read(10,'(a)') sbf(1:LBUF)
      if(sbf(1:1).eq.'#') goto 2
      i = 0
   10 i = i + 1
      if(sbf(i:i).eq.' ') goto 10
      pat = '        '
      j = 1
   11 n = ichar(sbf(i:i))
      if(n.ge.ichar('a').and.n.le.ichar('z'))
     *    n = n - ichar('a') + ichar('A')
      pat(j:j) = char(n)
      if(j.lt.8) j = j + 1
      i = i + 1
      if(sbf(i:i).ne.' ') goto 11
   12 i = i + 1
      if(sbf(i:i).eq.' ') goto 12
      if(pat.eq.'HPOLYG  ') then
        if(sbf(i:i).eq.'*') then
          ncct = -1
          read(sbf(i+1:LBUF),*) m, deep
        else
          read(sbf(i:LBUF),*) ncct, m, deep
          if (ncct.ne.ncc) call abendm('invalid parameter')
        endif
      else if(pat.eq.'SPOLYG  ') then
        if(sbf(i:i).eq.'*') then
          ncct = -1
          read(sbf(i+1:LBUF),*) m
        else
          read(sbf(i:LBUF),*) ncct, m
          if (ncct.ne.ncc) call abendm('invalid parameter')
        endif
      else
        call abendm('invalid parameter')
      endif
      if(m.le.2.or.m.gt.NPT) call abendm('invalid parameter')
      read(10,'(a)') sbf(1:LBUF)
      i = 1
      do 3 n=1,m
        goto 32
   31   read(10,'(a)') sbf(1:LBUF)
        i = 1
   32   if (ncct.eq.-1) then
          call scandm(i, sbf, xpt(n))
        else
          call scanfv(i, sbf, xpt(n))
        endif
        if (i.lt.0) call abendm('invalid parameter')
        if (i.eq.0) goto 31
        goto 34
   33   read(10,'(a)') sbf(1:LBUF)
        i = 1
   34   if (ncct.eq.-1) then
          call scandm(i, sbf, ypt(n))
        else
          call scanfv(i, sbf, ypt(n))
        endif
        if (i.lt.0) call abendm('invalid parameter')
        if (i.eq.0) goto 33
        if(pat.eq.'SPOLYG  ') then
          goto 36
   35     read(10,'(a)') sbf(1:LBUF)
          i = 1
   36     call scanfv(i, sbf, dpt(n))
          if (i.lt.0) call abendm('invalid parameter')
          if (i.eq.0) goto 35
        else
          dpt(n) = deep
        endif
    3 continue
      close(10)
      write(out,'(a,i4,a)') '  total ', m, ' points.'
      call premsg(out)
      if(ncct.eq.-1) then
        call cvinit(ncc,float(iorg),float(korg),
     *                  float(ispa),float(ispb))
        do 4 n=1,m
          call cviken(xpt(n), ypt(n), ye, xn)
          xpt(n) = xn
          ypt(n) = ye
    4   continue
      endif
      do 40 n=1,m
        xpt(n) = xpt(n) * 1000.
        ypt(n) = ypt(n) * 1000.
   40 continue
      xpt(m+1) = xpt(1)
      ypt(m+1) = ypt(1)
      dpt(m+1) = dpt(1)
c
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      if(ncct.ne.-1.and.ncct.ne.ncc) call abendm('invalid parameter')
      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)')
     *     '--- obs. mag.   ', xs, ys, mszx, mszy, mx, my, salt
      call premsg(out)
c
      if(alt.eq.0.) goto 5
      do 50 j=1,my
      do 50 i=1,mx
        h(i,j) = alt
   50 continue
      hnul = 9999.9
      goto 6
    5 read(12,'(a)') buf
      if(buf(1:1).eq.'#') goto 5
      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   '
    6 close(12)
      call premsg('  obs.surface : ' // strt)
c
      call premsg('Ambient field direction (in degrees) ?')
      call gparmf2('  Inclin., Declin.  : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90. .or.
     *   dec.lt.-180..or.dec.ge.360.) call abendm('invalid value')
      cx = cos(dip/RAD) * cos(dec/RAD)
      cy = cos(dip/RAD) * sin(dec/RAD)
      cz = sin(dip/RAD)
      call premsg('Range for loop-current estimation ?')
      call gparmf2('  SW corner offset to N, and N-S size (in km) ==> ',
     *             rxs, rxd)
      szmx = float(mszx)/1000.
      irxs = ifix(rxs/szmx + 1.5)
      nrx = ifix(rxd/szmx + 1.5)
      irxt = irxs + nrx - 1
      if (irxs.le.0.or.nrx.le.3.or.irxt.gt.mx)
     *     call abendm('ranges conflict')
      call gparmf2('  SW corner offset to E, and E-W size (in km) ==> ',
     *             rys, ryd)
      szmy = float(mszy)/1000.
      irys = ifix(rys/szmy + 1.5)
      nry = ifix(ryd/szmy + 1.5)
      iryt = irys + nry - 1
      if (irys.le.0.or.nry.le.3.or.iryt.gt.my)
     *     call abendm('ranges conflict')
      rxs = xs + rxs
      rys = ys + rys
c
      call gparma('Enter LCE.Corr. output filename   ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(13,file=fnam3,status='new')
      call premsg('Enter Auxiary(L) 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)') 'LCECORR : loop-current effect correction'
      write(6,'(a,a)') ' Obs. Mag. infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)') ' Loop-Loc data file : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)') ' LCE.Corr. outfile  : ', fnam3(1:lrtrim(fnam3))
      if(kout.eq.1) write(6,'(a,a)')
     *                 ' Auxiary(L) outfile : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs, ys, mszx, mszy, mx, my, salt
      write(6,1500) rxs, rys, mszx, mszy, nrx, nry, '    var.'
      write(6,2000) strt, dip, dec
c
      do 100 j=1,my
        y0 = float(iys + (j-1)*mszy)
        do 100 i=1,mx
          x0 = float(ixs + (i-1)*mszx)
          if(h(i,j).eq.hnul) then
            if(f(i,j).ne.vnul) call abendm('obs. alt. undefined')
            g(i,j) = vnul
          else
            v = 0.
            do 110 k=1,m
              x1 = xpt(k) - x0
              y1 = ypt(k) - y0
              z1 = h(i,j) + dpt(k)
              x2 = xpt(k+1) - x0
              y2 = ypt(k+1) - y0
              z2 = h(i,j) + dpt(k+1)
              call vbcal(x1,y1,z1, x2,y2,z2, x,y,z)
              v = v + (x*cx + y*cy + z*cz)*100.
  110       continue
            g(i,j) = v
          endif
  100 continue
c
      if(kout.eq.1) then
        write(8,'(a8,i4,4x,4i8)') 'LCE-100A', 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)') 'AUXL <== LCE-100A : loop-current (100A) effect'
      endif
c
      ft = 0.
      gt = 0.
      n = 0
      do 200 j=irys,iryt
      do 200 i=irxs,irxt
        if(f(i,j).eq.vnul) goto 200
        ft = ft + f(i,j)
        gt = gt + g(i,j)
        n = n + 1
  200 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 300 j=irys,iryt
        y = float(j) - cmy
        do 300 i=irxs,irxt
          x = float(i) - cmx
          if(f(i,j).eq.vnul) goto 300
          vs(1) = f(i,j) - ft
          vs(2) = g(i,j) - gt
          call sm2ex(x, y, vs)
  300 continue
      call sm2cls(coef, 3, 2)
      write(6,3000)
      write(6,4000) 'observation', coef(1,1)+ft, coef(2,1), coef(3,1)
      write(6,4000) '100A effect', coef(1,2)+gt, coef(2,2), coef(3,2)
c
      do 400 j=1,my
      do 400 i=1,mx
        fr(i,j) = vnul
        gr(i,j) = vnul
  400 continue
      fg = 0.
      gg = 0.
      do 440 j=irys,iryt
        y = float(j) - cmy
        do 440 i=irxs,irxt
          x = float(i) - cmx
          call sm2rv(x, y, vs)
          if(f(i,j).eq.vnul) goto 440
          fr(i,j) = f(i,j) - (vs(1) + ft)
          gr(i,j) = g(i,j) - (vs(2) + gt)
          fg = fg + fr(i,j)*gr(i,j)
          gg = gg + gr(i,j)*gr(i,j)
  440 continue
      gm = fg / gg
      write(6,5000) gm*100.
c
      do 500 j=irys,iryt
      do 500 i=irxs,irxt
        if(gr(i,j).ne.vnul) gr(i,j) = gr(i,j) * gm
  500 continue
      do 550 j=1,my
      do 550 i=1,mx
        if(g(i,j).ne.vnul) g(i,j) = g(i,j) * gm
        if(f(i,j).ne.vnul) f(i,j) = f(i,j) - g(i,j)
  550 continue
      write(13,'(a)') '# loop-current effect corrected mag. anom.'
      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 600 j=1,my
        write(13,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  600 continue
      if(alt.eq.0.) then
        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 700 j=1,my
          write(13,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  700   continue
      endif
      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)') 'AUXL <== Obs-Trnd : Trend-removed Observation'
        write(8,'(a8,i4,4x,4i8)') 'LC-Efect', 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)') 'AUXL <== LC-Efect : Simulated LC.effect'
        write(8,'(a8,i4,4x,4i8)') 'LCE-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)') 'AUXL <== LCE-Trnd : Trend-removed LC.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)') 'LCECORR : loop-current effect correction'
        write(99,'(a,a)')
     *       ' Obs. Mag. infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' Loop-Loc data file : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' LCE.Corr. outfile  : ', fnam3(1:lrtrim(fnam3))
        if(kout.eq.1) write(99,'(a,a)')
     *       ' Auxiary(L) outfile : ', fnam4(1:lrtrim(fnam4))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xs, ys, mszx, mszy, mx, my, salt
        write(99,1500) rxs, rys, mszx, mszy, nrx, nry, '    var.'
        write(99,2000) strt, dip, dec
        write(99,3000)
        write(99,4000) 'observation', coef(1,1)+ft,coef(2,1),coef(3,1)
        write(99,4000) '100A effect', coef(1,2)+gt,coef(2,2),coef(3,2)
        write(99,5000) gm*100.
        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)
 1500 format(3x, 'LCE estim.  ', 2f10.2, 2i8, 2i6, a8)
 2000 format('  obs. surface  : ', a25/
     *       '  ambient field :  dip=', f6.1, 3x, 'dec=', f7.1)
 3000 format(' Trend surface in nT'/ 5x, '( x/y : N/E coord. ',
     *       'in mesh unit referred to the center of area )')
 4000 format('  for ', a11, ' : ', f8.2, '  +  ',
     *       f7.3, ' * x  +  ', f7.3, ' * y')
 5000 format(' Intensity of loop-current : ', f7.1, ' A')
      end
c---
      subroutine vbcal(x1,y1,z1, x2,y2,z2, x,y,z)
      rrx = y1*z2 - y2*z1
      rry = z1*x2 - z2*x1
      rrz = x1*y2 - x2*y1
      rr1 = x1*x1 + y1*y1 + z1*z1
      rr2 = x2*x2 + y2*y2 + z2*z2
      p = x1*x2 + y1*y2 + z1*z2
      r1 = sqrt(rr1)
      r2 = sqrt(rr2)
	  rr = r1 * r2
	  f = (100./rr) * ((r1+r2)/(rr+p))
      x = rrx * f
      y = rry * f
      z = rrz * f
      return
      end
c---
      subroutine scandm(l, buf, fv)
      parameter (LBUF=255)
      character*(LBUF) buf
      i = l
   10 if (buf(i:i).ne.' ' .and. buf(i:i).ne.char(8)) goto 1
      i = i + 1
      if (i.gt.LBUF) then
        l = 0
        return
      endif
      goto 10
    1 is = i
      j = index(':+-0123456789', buf(i:i))
      if (j.eq.0) goto 9
      if (j.eq.1) goto 3
    2 i = i + 1
      j = index(':0123456789', buf(i:i))
      if (j.eq.0) goto 9
      if (j.gt.1) goto 2
    3 it = i - 1
      if (it.lt.is) then
        ndv = 0
      else
        read(buf(is:it),*) ndv
      endif
    4 i = i + 1
      if (buf(i:i).ne.' ' .and. buf(i:i).ne.char(8)) goto 20
      if (i.eq.LBUF) goto 9
      goto 4
   20 is = i
      j = index('.+-0123456789', buf(i:i))
      if (j.eq.0) goto 7
      if (j.eq.1) goto 6
    5 i = i + 1
      j = index('.0123456789', buf(i:i))
      if (j.eq.0) goto 7
      if (j.gt.1) goto 5
    6 i = i + 1
      j = index('0123456789', buf(i:i))
      if (j.ne.0) goto 6
    7 it = i - 1
      if (it.lt.is) goto 9
      read(buf(is:it),*) fv
      fv = fv + float(ndv*60)
      l = i
      return
    9 l = -1
      return
      end
c---
      subroutine scanfv(l, buf, fv)
      parameter (LBUF=255)
      character*(LBUF) buf
      i = l
    1 if (buf(i:i).ne.' ' .and. buf(i:i).ne.char(8)) goto 2
      i = i + 1
      if (i.gt.LBUF) then
        l = 0
        return
      endif
      goto 1
    2 is = i
      j = index('.+-0123456789', buf(i:i))
      if (j.eq.0) goto 5
      if (j.eq.1) goto 4
    3 i = i + 1
      j = index('.0123456789', buf(i:i))
      if (j.eq.0) goto 5
      if (j.gt.1) goto 3
    4 i = i + 1
      j = index('0123456789', buf(i:i))
      if (j.ne.0) goto 4
    5 it = i - 1
      if (it.lt.is) then
        l = -1
        return
      endif
      read(buf(is:it),*) fv
      l = i
      return
      end
