!--------( LCECORR : loop-current effect correction )--------
  implicit none
  integer :: lwkdir, lrtrim
  integer,parameter :: LBUF = 255
  real,parameter :: RAD = (180./3.141593)
  real,allocatable :: f(:,:), g(:,:), h(:,:), fr(:,:), gr(:,:)
  real,allocatable :: xpt(:), ypt(:), dpt(:)
  real :: coef(3,2), vs(2)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, buf
  character(LBUF+1) :: sbf;  character(72) :: out
  character(25) :: strt;  character(20) :: sdtm
  character(8) :: area, area1, pat, salt
  integer :: ldr, i, j, n, m, k, kout
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my, ncct
  integer :: ix1, iy1, mszx1, mszy1, mx1, my1
  integer :: irxs, irxt, nrx, irys, iryt, nry
  real :: xs, ys, alt, vnul, hnul, deep, xn, ye, dip, dec, cx, cy, cz
  real :: rxs, rxd, szmx, rys, ryd, szmy, x0, y0, v
  real :: x1, y1, z1, x2, y2, z2, x, y, z, ft, gt, fg, gg, gm
  real :: an, cmx, cmy

  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')
  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'

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

  xs = float(ixs)/1000.;  ys = float(iys)/1000.
  read(12,*) ((f(i,j),i=1,mx),j=1,my)
  if (alt < 0.) call abendm('Undefined Altitude')

  call gparma('Enter Loop-Location data filename ==> ', 81-ldr, fnam2(ldr:80))
  open(10,file=fnam2,status='old')
  do
    read(10,'(a)') sbf(1:LBUF)
    if (sbf(1:1) /= '#') exit
  enddo
  i = 0
  do
    i = i + 1;  if (sbf(i:i) /= ' ') exit
  enddo
  j = 1;  pat = '        '
  do
    n = ichar(sbf(i:i))
    if (n == 0) call abendm('invalid parameter (0)')
    if ((n >= ichar('a')) .and. (n <= ichar('z')))   &
      n = n - ichar('a') + ichar('A')
    pat(j:j) = char(n)
    if (j < 8) j = j + 1
    i = i + 1;  if (sbf(i:i) == ' ') exit
  enddo
  do
    i = i + 1;  if (sbf(i:i) /= ' ') exit
  enddo
  if (pat == 'HPOLYG  ') then
    if (sbf(i:i) == '*') then
      ncct = -1
      read(sbf(i+1:LBUF),*) m, deep
    else
      read(sbf(i:LBUF),*) ncct, m, deep
      if (ncct /= ncc) call abendm('invalid parameter (ncc)')
    endif
  else if (pat == 'SPOLYG  ') then
    if (sbf(i:i) == '*') then
      ncct = -1
      read(sbf(i+1:LBUF),*) m
    else
      read(sbf(i:LBUF),*) ncct, m
      if (ncct /= ncc) call abendm('invalid parameter (ncc)')
    endif
  else
    call abendm('invalid parameter (pat)')
  endif
  if (m <= 2) call abendm('invalid parameter (m)')

  allocate(xpt(m));  allocate(ypt(m));  allocate(dpt(m))

  i = 0
  do n=1,m
    do
      if (i == 0) then;  read(10,'(a)') sbf(1:LBUF); i = 1;  endif
      if (ncct == -1) then
        call scandm(i, sbf, xpt(n))
      else
        call scanfv(i, sbf, xpt(n))
      endif
      if (i < 0) call abendm('invalid parameter')
      if (i > 0) exit
    enddo
    do
      if (i == 0) then;  read(10,'(a)') sbf(1:LBUF); i = 1;  endif
      if (ncct == -1) then
        call scandm(i, sbf, ypt(n))
      else
        call scanfv(i, sbf, ypt(n))
      endif
      if (i < 0) call abendm('invalid parameter')
      if (i > 0) exit
    enddo
    if (pat == 'SPOLYG  ') then
      do
        if (i == 0) then;  read(10,'(a)') sbf(1:LBUF); i = 1;  endif
        call scanfv(i, sbf, dpt(n))
        if (i < 0) call abendm('invalid parameter')
        if (i > 0) exit
      enddo
    else
      dpt(n) = deep
    endif
  enddo
  close(10)
  write(out,'(a,i4,a)') '  total ', m, ' points.'
  call premsg(out)
  if (ncct == -1) then
    call cvinit(ncc,float(iorg),float(korg), float(ispa),float(ispb))
    do n=1,m
      call cviken(xpt(n), ypt(n), ye, xn)
      xpt(n) = xn;  ypt(n) = ye
    enddo
  endif
  do n=1,m
    xpt(n) = xpt(n) * 1000.; ypt(n) = ypt(n) * 1000.
  enddo
  xpt(m+1) = xpt(1);  ypt(m+1) = ypt(1);  dpt(m+1) = dpt(1)

  write(out,'(a,a8,a,i4,4x,4i8)')   &
    '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  if ((ncct /= -1) .and. (ncct /= ncc)) call abendm('invalid parameter')
  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)')   &
    '--- obs. mag.   ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)

  if (alt /= 0.) then
    do j=1,my; do i=1,mx;  h(i,j) = alt;  enddo; enddo
    hnul = 9999.9
    write(strt,'(f7.0,a)') alt, ' m above sea level'
  else
    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; nc = 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   '
  endif
  close(12)
  call premsg('  obs.surface : ' // strt)

  call premsg('Ambient field direction (in degrees) ?')
  call gparmf2('  Inclin., Declin.  : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.) .or.   &
      (dec < -180.) .or. (dec >= 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 <= 0) .or. (nrx <= 3) .or. (irxt > 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 <= 0) .or. (nry <= 3) .or. (iryt > my))   &
    call abendm('ranges conflict')

  call gparma('Enter LCE.Corr. output filename   ==> ', 81-ldr, fnam3(ldr:80))
  open(13,file=fnam3,status='new')
  call premsg('Enter Auxiliary(L) output filename')
  call gparma('     ( blank to suppress output ) ==> ', 81-ldr, fnam4(ldr:80))
  call clspin()
  call strdtm(sdtm)
  kout = 0
  if (fnam4(ldr:ldr) /= ' ') 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 == 1)   &
    write(6,'(a,a)') ' Aux(L) output file : ', 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,'(3x,a,a25)') 'obs. surface  : ', strt
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'ambient field : ', 'dip=', dip, 'dec=', dec
  write(6,'(3x,a)') 'range for LCE estimation : '
  write(6,'(5x,a,2f8.1)') 'SW-corner offsetN, and N-S size =', rxs, rxd
  write(6,'(5x,a,2f8.1)') 'SW-corner offsetE, and E-W size =', rys, ryd

  do j=1,my
    y0 = float(iys + (j-1)*mszy)
    do i=1,mx
      x0 = float(ixs + (i-1)*mszx)
      if (h(i,j) == hnul) then
        if (f(i,j) /= vnul) call abendm('obs. alt. undefined')
        g(i,j) = vnul
      else
        v = 0.
        do 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.
        enddo
        g(i,j) = v
      endif
    enddo
  enddo

  if (kout == 1) then
    write(8,'(a8,4x,i4,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 j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUXL <== LCE-100A : loop-current (100A) effect'
  endif

  ft = 0.;  gt = 0.;  n = 0
  do j=irys,iryt; do i=irxs,irxt
    if (f(i,j) == vnul) cycle
    ft = ft + f(i,j);  gt = gt + g(i,j)
    n = n + 1
  enddo; enddo
  an = float(n);  ft = ft / an;  gt = gt / an

  call sm2opn(1, 2)
  cmx = float(mx+1) / 2.
  cmy = float(my+1) / 2.
  do j=irys,iryt
    y = float(j) - cmy
    do i=irxs,irxt
      x = float(i) - cmx
      if (f(i,j) == vnul) cycle
      vs(1) = f(i,j) - ft
      vs(2) = g(i,j) - gt
      call sm2ex(x, y, vs)
    enddo
  enddo
  call sm2cls(coef, 3, 2)
  write(6,'(1x,a)') 'Trend surface in nT'
  write(6,'(5x,a)')   &
    '( x/y : N/E coord. in mesh unit referred to the center of area )'
  write(6,'(1x,a,f8.2,a,f7.3,a,f7.3,a)') ' for observation : ',   &
    coef(1,1)+ft, '  +  ', coef(2,1), ' * x  +  ', coef(3,1), ' * y'
  write(6,'(1x,a,f8.2,a,f7.3,a,f7.3,a)') ' for 100A effect : ',   &
    coef(1,2)+gt, '  +  ', coef(2,2), ' * x  +  ', coef(3,2), ' * y'

  do j=1,my; do i=1,mx;  fr(i,j) = vnul; gr(i,j) = vnul;  enddo; enddo
  fg = 0.;  gg = 0.
  do j=irys,iryt
    y = float(j) - cmy
    do i=irxs,irxt
      x = float(i) - cmx
      call sm2rv(x, y, vs)
      if (f(i,j) == vnul) cycle
      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)
    enddo
  enddo
  gm = fg / gg
  write(6,'(1x,a,f7.1,a)') 'Intensity of loop-current : ', gm*100., ' A'

  do j=irys,iryt; do i=irxs,irxt
    if (gr(i,j) /= vnul) gr(i,j) = gr(i,j) * gm
  enddo; enddo
  do j=1,my; do i=1,mx
    if (g(i,j) /= vnul) g(i,j) = g(i,j) * gm
    if (f(i,j) /= vnul) f(i,j) = f(i,j) - g(i,j)
  enddo; enddo
  write(13,'(a)') '# loop-current effect corrected mag. anom.'
  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
  if (alt == 0.) then
    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
  endif
  close(13)

  if (kout == 1) then
    write(8,'(a8,4x,i4,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 j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (fr(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUXL <== Obs-Trnd : Trend-removed Observation'
    write(8,'(a8,4x,i4,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 j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUXL <== LC-Efect : Simulated LC.effect'
    write(8,'(a8,4x,i4,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 j=1,my
      write(8,'((f7.1,9(1x,f7.1)))') (gr(i,j),i=1,mx)
    enddo
    write(6,'(a)') 'AUXL <== LCE-Trnd : Trend-removed LC.effect'
    close(8)
  endif

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    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 == 1) write(99,'(a,a)')   &
      ' Aux(L) output file : ', 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,'(3x,a,a25)') 'obs. surface  : ', strt
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'ambient field : ', 'dip=', dip, 'dec=', dec
    write(99,'(3x,a)') 'range for LCE estimation : '
    write(99,'(5x,a,2f8.1)') 'SW-corner offsetN, and N-S size =', rxs, rxd
    write(99,'(5x,a,2f8.1)') 'SW-corner offsetE, and E-W size =', rys, ryd
    write(99,'(1x,a)') 'Trend surface in nT'
    write(99,'(5x,a)')   &
      '( x/y : N/E coord. in mesh unit referred to the center of area )'
    write(99,'(1x,a,f8.2,a,f7.3,a,f7.3,a)') ' for observation : ',   &
      coef(1,1)+ft, '  +  ', coef(2,1), ' * x  +  ', coef(3,1), ' * y'
    write(99,'(1x,a,f8.2,a,f7.3,a,f7.3,a)') ' for 100A effect : ',   &
      coef(1,2)+gt, '  +  ', coef(2,2), ' * x  +  ', coef(3,2), ' * y'
    write(99,'(1x,a,f7.1,a)') 'Intensity of loop-current : ', gm*100., ' A'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end
!---
  subroutine vbcal(x1,y1,z1, x2,y2,z2, x,y,z)
  implicit none
  real,intent(in) :: x1, y1, z1, x2, y2, z2
  real,intent(out) :: x, y, z
  real :: rrx, rry, rrz, rr1, rr2, p, r1, r2, rr, f

  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 subroutine
!---
  subroutine scandm(l, sbf, f)
  implicit none
  integer,parameter :: LBUF = 255
  character(LBUF+1),intent(in) :: sbf
  integer,intent(inout) :: l
  real,intent(out) :: f
  integer :: i, is, it, j, ndv
  real :: fv

  i = l
  do
    if ((sbf(i:i) /= ' ') .and. (sbf(i:i) /= char(8))) exit
    i = i + 1
    if (i > LBUF) then;  l = 0; return;  endif
  enddo
  is = i;  j = index(':+-0123456789', sbf(i:i))
  do
    if (j == 0) then;  l = -1; return;  endif
    if (j == 1) exit
    i = i + 1;  j = index(':0123456789', sbf(i:i))
  enddo
  it = i - 1
  if (it < is) then
    ndv = 0
  else
    read(sbf(is:it),*) ndv
  endif
  do
    i = i + 1
    if ((sbf(i:i) /= ' ') .and. (sbf(i:i) /= char(8))) exit
    if (i == LBUF) then;  l = -1; return;  endif
  enddo
  is = i;  j = index('.+-0123456789', sbf(i:i))
  if (j > 0) then
    do while (j > 1)
      i = i + 1;  j = index('.0123456789', sbf(i:i))
    enddo
    do while (j /= 0)
      i = i + 1;  j = index('0123456789', sbf(i:i))
    enddo
  endif
  it = i - 1
  if (it < is) then;  l = -1; return;  endif
  read(sbf(is:it),*) fv
  f = fv + float(ndv*60)
  l = i
  return
  end subroutine
!---
  subroutine scanfv(l, sbf, f)
  implicit none
  integer,parameter :: LBUF = 255
  character(LBUF+1),intent(in) :: sbf
  integer,intent(inout) :: l
  real,intent(out) :: f
  integer :: i, is, it, j

  i = l
  do
    if ((sbf(i:i) /= ' ') .and. (sbf(i:i) /= char(8))) exit
    i = i + 1
    if (i > LBUF) then;  l = 0; return;  endif
  enddo
  is = i;  j = index('.+-0123456789', sbf(i:i))
  if (j > 0) then
    do while (j > 1)
      i = i + 1;  j = index('.0123456789', sbf(i:i))
    enddo
    do while (j /= 0)
      i = i + 1;  j = index('0123456789', sbf(i:i))
    enddo
  endif
  it = i - 1
  if (it < is) then;  l = -1; return;  endif
  read(sbf(is:it),*) f
  l = i
  return
  end subroutine