c--------( GOVLAY : gen. combined grid with slit )--------
      parameter (LXN=4501, LYN=4501, LXO=4501, LYO=4501)
c
      dimension f(LXO,LYO), a(LXO,LYO), g(LXN,LYN), h(LXN,LYN)
      character*80 wdr, flog, fnam1, fnam2
      character area*8, narea*8, area1*8, anul*10, salt*8
      character buf*80, sdtm*20
c
      call premsg('GOVLAY : gen. combined grid (with slit) '
     *            // 'for multiple areas')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call premsg('')
      call gparma('Enter Name-label for new grid ==> ', 8, narea)
      call gparmi('Select Projection-Number ==> ', nccn)
      ncn = nccn
      if(nccn.ge.200) ncn = nccn - 200
      ispan = 0
      ispbn = 0
      if(ncn.ge.1.and.ncn.le.62) then
        iorgn = 0
        korgn = (ncn-30)*360 - 180
        if(ncn.gt.60) korgn = 0
        if(ncn.eq.61) iorgn = +5400
        if(ncn.eq.62) iorgn = -5400
      else if(ncn.eq.65) then
        iorgn = 0
        call gparmi('  Central Longitude ?  (Deg.) : ', i1)
        korgn = i1*60
      else
        call gparmi2('  Origin Latitude ?   (Deg. Min.) : ', i1,i2)
        iorgn = i1*60 + i2
        call gparmi2('  Origin Longitude ?  (Deg. Min.) : ', i1,i2)
        korgn = i1*60 + i2
        if(ncn.eq.72) then
          call gparmi2('  St.Parallel-1 Lat. ?  (Deg. Min.) : ', i1,i2)
          ispan = i1*60 + i2
          call gparmi2('  St.Parallel-2 Lat. ?  (Deg. Min.) : ', i1,i2)
          ispbn = i1*60 + i2
        endif
      endif
      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 gparmi('  Unit Mesh size in m :  ', msz)
      if(msz.le.0) call abendm('invalid value')
      call gparmi('    Number of Mesh to N :  ', mxn)
      if(mxn.gt.LXN) call abendm('too large array size')
      call gparmi('    Number of Mesh to E :  ', myn)
      if(myn.gt.LYN) call abendm('too large array size')
      vnuln = 9999.9
      call gparma('Enter new grid output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(2,file=fnam2,status='new')
      do 100 jp=1,myn
      do 100 ip=1,mxn
        g(ip,jp) = vnuln
        h(ip,jp) = vnuln
  100 continue
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'GOVLAY : gen. combined grid (with slit) '
     *               // 'for multiple areas'
      write(6,'(a,a)') ' Output new file : ', fnam2
      write(6,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *              xl, yl, msz, msz, mxn, myn, '    var.'
      call strdtm(sdtm)
      open(9,status='scratch')
c
    1 call premsg('')
      call premsg('Enter Input source data filename')
      call gparma('             (Hit RETURN to End) ==> ',
     *            81-ldr, fnam1(ldr:80))
      if(fnam1(ldr:ldr).eq.' ') goto 8
      open(1,file=fnam1,status='old')
    2 read(1,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      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
      if(ncc.ne.nccn .or.
     *   iorg.ne.iorgn.or.korg.ne.korgn .or.
     *   ispa.ne.ispan.or.ispb.ne.ispbn)
     *     call abendm('Header inconsistent')
      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'
      if (mszx.ne.msz.or.mszy.ne.msz)
     *     call abendm('illegal mesh-size')
      is = (ixs-ixl)/mszx
      js = (iys-iyl)/mszy
      if (is*mszx+ixl.ne.ixs .or.
     *    js*mszy+iyl.ne.iys) call abendm('ranges conflict')
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      read(1,*) ((f(i,j),i=1,mx),j=1,my)
      if(alt.eq.0.) goto 3
      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
      goto 4
    3 read(1,'(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(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)
     *     call abendm('Header2 inconsistent')
      write(anul,'(f10.3)') vnul1
      if(anul(9:10).ne.'00'.or.anul(1:1).ne.' ')
     *     call abendm('unable to handle vnul1 value')
      read(1,*) ((a(i,j),i=1,mx),j=1,my)
    4 close(1)
      write(6,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      write(6,'(4x,a,2a10,2x,2a8,2a6,a7)')
     *     '  (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *     'alt'
      write(6,'(4x,a,2f10.2,2i8,2i6,a8)')
     *     'source data ', xs, ys, mszx, mszy, mx, my, salt
      write(9,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(9,*) xs, ys, mszx, mszy, mx, my, alt
      do 300 j=1,my
      do 300 i=1,mx
        if (f(i,j).eq.vnul) goto 300
        do 200 jp=js+j-1,js+j+1
        do 200 ip=is+i-1,is+i+1
          if (ip.le.0.or.ip.gt.mxn.or.jp.le.0.or.jp.gt.myn) goto 200
          g(ip,jp) = vnuln
          h(ip,jp) = vnuln
  200   continue
  300 continue
      do 400 j=1,my
        jp = js + j
        do 400 i=1,mx
          if (f(i,j).eq.vnul) goto 400
          ip = is + i
          if (ip.le.0.or.ip.gt.mxn.or.jp.le.0.or.jp.gt.myn) goto 400
          g(ip,jp) = f(i,j)
          if (alt.ne.0.) then
            h(ip,jp) = alt
          else if (a(i,j).ne.vnul1) then
            h(ip,jp) = a(i,j)
          else
            write(buf,'(a30,f8.1,a8,i4,a4,i4)') 
     *             '  a(i,j)=vnul  where  f(i,j) =', f(i,j),
     *             '  at  i=', i, '  j=', j
            call premsg(buf(1:58))
            call abendm('Aborted')
          endif
  400 continue
      goto 1
c
    8 call clspin()
      endfile(9)
      rewind(9)
      write(2,'(a)') '# prog.GOVLAY output'
      write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl,iyl, mszx,mszy, mxn,myn, vnuln, 0.
      do 500 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn)
  500 continue
      write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl,iyl, mszx,mszy, mxn,myn, vnuln, -1.
      do 600 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mxn)
  600 continue
      close(2)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a)') 'GOVLAY : gen. combined grid (with slit) '
     *                    // 'for multiple areas'
        write(99,'(a,a)') ' Output new file : ', fnam2
        write(99,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *                xl, yl, msz, msz, mxn, myn, '    var.'
        write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
   92   read(9,'(a8,i4,4x,4i8)',end=93)
     *       area, ncc, iorg, korg, ispa, ispb
        read(9,*) xs, ys, mszx, mszy, mx, my, alt
        write(salt,'(1x,f7.0)') alt
        if (alt.eq.0.) salt = '    var.'
        write(99,'(a,a,a8,a,i4,4x,4i8)') '--- ',
     *       'Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
        write(99,'(4x,a,2a10,2x,2a8,2a6,a7)')
     *       '  (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *       'alt'
        write(99,'(4x,a,2f10.2,2i8,2i6,a8)')
     *       'source data ', xs, ys, mszx, mszy, mx, my, salt
        goto 92
   93   call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      close(9)
c
      stop
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'gen.new grid', 2f10.2, 2i8, 2i6, a8)
      end
