c--------( SEL : generate grid-data for selected area )--------
      parameter (LXO=4001, LYO=4001, LXN=4001, LYN=4001)
c
      dimension f(LXO,LYO), g(LXN,LYN)
      character*80 wdr, flog, fnam1, fnam2
      character area*8, area1*8, narea*8, salt*8, saltn*8
      character buf*80, out*72
c
      call premsg('SEL : generate grid-data for selected area')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Input source data filename ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(1,file=fnam1,status='old')
    1 read(1,'(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(1,'(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'
      rewind(1)
c
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      ixt = ixs + (mx-1)*mszx
      iyt = iys + (my-1)*mszy
      smx = float(mszx)
      smy = float(mszy)
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      write(buf,'(f9.2)') vnul
      if(buf(9:9).ne.'0'.or.buf(1:1).ne.' ')
     *     call abendm('unable to handle vnul value')
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      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)')
     *     '--- source data ', xs, ys, mszx, mszy, mx, my, salt
      call premsg(out)
      saltn = salt
      if(alt.lt.0.) then
        call gparmf('  Altitude [in m] for this data ==> ', alt)
        if (alt.le.0.) call abendm('invalid value')
        write(saltn,'(1x,f7.0)') alt
      endif
      call premsg('')
      call gparma('Enter Name-label for new grid ==> ', 8, narea)
      call premsg('Specify gridding parameters')
      call premsg('  Southwest corner Coord. in km ?')
      call gparmf('    Northing : ', xl)
      ixl = nint(xl*1000.)
      call gparmf('    Easting  : ', yl)
      iyl = nint(yl*1000.)
      call premsg('  Mesh size and numbers to N ?')
      call gparmi('    Unit Mesh size in m :  ', mszxn)
      if(mszxn.le.0) call abendm('invalid value')
      call gparmi('    Number of Mesh to N :  ', mxn)
      if(mxn.gt.LXN) call abendm('too large array size')
      ixh = ixl + (mxn-1)*mszxn
      call premsg('  Mesh size and numbers to E ?')
      call gparmi('    Unit Mesh size in m :  ', mszyn)
      if(mszyn.le.0) call abendm('invalid value')
      call gparmi('    Number of Mesh to E :  ', myn)
      if(myn.gt.LYN) call abendm('too large array size')
      iyh = iyl + (myn-1)*mszyn
      call gparma('Enter new grid output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      call clspin()
      open(2,file=fnam2,status='new')
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'SEL : generate grid-data for selected area'
      write(6,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
      write(6,1000) narea, ncc, iorg,korg, ispa,ispb,
     *              xs, ys, mszx,  mszy,  mx,  my,  salt,
     *              xl, yl, mszxn, mszyn, mxn, myn, saltn
      n = 0
c
    2 read(1,'(a)') buf
      if(buf(1:1).ne.'#') goto 3
      write(2,'(a)') buf(1:lrtrim(buf))
      goto 2
    3 n = n + 1
      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(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
      if(ixs1.ne.ixs.or.iys1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my .or. vnul1.ne.vnul)
     *     call abendm('Header2 inconsistent')
      write(2,'(2a)') '# Source : ', fnam1(1:lrtrim(fnam1))
      write(2,'(a)')  '# prog.SEL applied to above file'
      write(2,'(a8,i4,4x,4i8)') narea, ncc, iorg,korg, ispa,ispb
      if(n.eq.2) alt = -1.
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl,iyl, mszxn,mszyn, mxn,myn, vnul, alt
c
      read(1,*) ((f(i,j),i=1,mx),j=1,my)
      do 100 j=1,myn
        iy = iyl + (j-1)*mszyn
        jp = (iy-iys)/mszy
        dy = float(iy-iys-jp*mszy)/smy
        do 100 i=1,mxn
          ix = ixl + (i-1)*mszxn
          ip = (ix-ixs)/mszx
          dx = float(ix-ixs-ip*mszx)/smx
          if(ix.lt.ixs.or.ix.gt.ixt.or.iy.lt.iys.or.iy.gt.iyt) then
            g(i,j) = vnul
          else if(f(ip+1,jp+1).eq.vnul) then
            g(i,j) = vnul
          else if(dy.eq.0.) then
            if(dx.eq.0.) then
              g(i,j) = f(ip+1,jp+1)
            else if(f(ip+2,jp+1).eq.vnul) then
              g(i,j) = vnul
            else
              g(i,j) = f(ip+1,jp+1) + (f(ip+2,jp+1)-f(ip+1,jp+1))*dx
            endif
          else if(f(ip+1,jp+2).eq.vnul) then
            g(i,j) = vnul
          else
            if(dx.eq.0.) then
              g(i,j) = f(ip+1,jp+1) + (f(ip+1,jp+2)-f(ip+1,jp+1))*dy
            else if(f(ip+2,jp+1).eq.vnul.or.f(ip+2,jp+2).eq.vnul) then
              g(i,j) = vnul
            else
              f1 = f(ip+1,jp+1) + (f(ip+2,jp+1)-f(ip+1,jp+1))*dx
              f2 = f(ip+1,jp+2) + (f(ip+2,jp+2)-f(ip+1,jp+2))*dx
              g(i,j) = f1 + (f2-f1)*dy
            endif
          endif
  100 continue
      do 200 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn)
  200 continue
      if(n.eq.1.and.alt.eq.0.) goto 2
c
    8 close(1)
      close(2)
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,$)') 'SEL :'
        write(99,'(a)')   ' generate grid-data for selected area'
        write(99,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
        write(99,1000) narea, ncc, iorg,korg, ispa,ispb,
     *                 xs, ys, mszx,  mszy,  mx,  my,  salt,
     *                 xl, yl, mszxn, mszyn, mxn, myn, saltn
        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, 'source grid ', 2f10.2, 2i8, 2i6, a8/
     *       3x, 'gen.new grid', 2f10.2, 2i8, 2i6, a8)
      end
