c--------( GTOPO : generate topography grid-data from DEM )--------
      parameter (LXO=2401, LYO=2401, LXN=4001, LYN=4001)
      character*26 DEMDIR
      parameter (DEMDIR = '/home/SHARE/data/DEM/wm40/')
c
      dimension f1(LXO,LYO), f2(LXO,LYO), f3(LXO,LYO), g(LXN,LYN)
      character*80 wdr, flog, fnam1, fnam2
      character area*8, narea*8, buf*80, sdtm*20
      logical exist
c
      call premsg('GTOPO : generate topography grid-data from DEM')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call premsg('')
      call gparma('Enter Name-label for topo 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 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')
      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')
      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)') 'GTOPO : generate topography grid-data from DEM'
      write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
      write(6,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *              xl, yl, mszxn, mszyn, mxn, myn
c
      call cvinit(nccn, float(iorgn), float(korgn),
     *                  float(ispan), float(ispbn))
      nik1 = 0
      nik2 = 0
      nik3 = 0
      ksel = 3
      do 10 j=1,myn
        iy = iyl + (j-1)*mszyn
        ye = float(iy) / 1000.
        do 10 i=1,mxn
          ix = ixl + (i-1)*mszxn
          xn = float(ix) / 1000.
          if (nccn .lt. 200) then
            call cvenik(ye, xn, fit, fkt)
            call xtw84(fit, fkt, 0., fi, fk, alt)
          else
            call cvenik(ye, xn, fi, fk)
          endif
          ido = ifix(fi/60.)
          fif = fi - float(ido*60)
          kdo = ifix(fk/60.)
          fkf = fk - float(kdo*60)
          nik = ido*1000 + kdo
          key = 0
          if (nik .eq. nik1) then
            ksel = 1
          else if (nik .eq. nik2) then
            ksel = 2
          else if (nik .eq. nik3) then
            ksel = 3
          else
            write(fnam1,'(a,i5,a)') DEMDIR, nik, '.alt'
            inquire(file=fnam1,exist=exist)
            if (exist) then
              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
              if (ncc.ne.399) call abendm('Illegal DEM file spec.')
              if (iorg.ne.(ido*60) .or. korg.ne.(kdo*60))
     *            call abendm('DEM file broken')
              read(1,*) ixs,iys, mszx,mszy, mx,my, vnul, alt
              if (ixs.ne.0.or.iys.ne.0) call abendm('DEM file broken')
              if (mx.gt.LXO.or.my.gt.LYO)
     *            call abendm('too large array size')
              if (ksel.eq.3) then
                ksel = 1
                nik1 = nik
                uif1 = float(mszx) / 1000.
                ukf1 = float(mszy) / 1000.
                mmi1 = mx - 1
                mmk1 = my - 1
                read(1,*) ((f1(ii,jj),ii=1,mx),jj=1,my)
              else if (ksel.eq.1) then
                ksel = 2
                nik2 = nik
                uif2 = float(mszx) / 1000.
                ukf2 = float(mszy) / 1000.
                mmi2 = mx - 1
                mmk2 = my - 1
                read(1,*) ((f2(ii,jj),ii=1,mx),jj=1,my)
              else
                ksel = 3
                nik3 = nik
                uif3 = float(mszx) / 1000.
                ukf3 = float(mszy) / 1000.
                mmi3 = mx - 1
                mmk3 = my - 1
                read(1,*) ((f3(ii,jj),ii=1,mx),jj=1,my)
              endif
              close(1)
            else
              key = 1
              g(i,j) = -1.
            endif
          endif
          if (key.eq.0) then
            if (ksel .eq. 1) then
              si = fif / uif1
              mi = ifix(si)
              if (mi.eq.mmi1) mi = mmi1 - 1
              dif = si - float(mi)
              sk = fkf / ukf1
              mk = ifix(sk)
              if (mk.eq.mmk1) mk = mmk1 - 1
              dkf = sk - float(mk)
              a0 = f1(mi+1,mk+1)
              a1 = f1(mi+1,mk+2)
              a2 = f1(mi+2,mk+1)
              a3 = f1(mi+2,mk+2)
            else if (ksel.eq.2) then
              si = fif / uif2
              mi = ifix(si)
              if (mi.eq.mmi2) mi = mmi2 - 1
              dif = si - float(mi)
              sk = fkf / ukf2
              mk = ifix(sk)
              if (mk.eq.mmk2) mk = mmk2 - 1
              dkf = sk - float(mk)
              a0 = f2(mi+1,mk+1)
              a1 = f2(mi+1,mk+2)
              a2 = f2(mi+2,mk+1)
              a3 = f2(mi+2,mk+2)
            else
              si = fif / uif3
              mi = ifix(si)
              if (mi.eq.mmi3) mi = mmi3 - 1
              dif = si - float(mi)
              sk = fkf / ukf3
              mk = ifix(sk)
              if (mk.eq.mmk3) mk = mmk3 - 1
              dkf = sk - float(mk)
              a0 = f3(mi+1,mk+1)
              a1 = f3(mi+1,mk+2)
              a2 = f3(mi+2,mk+1)
              a3 = f3(mi+2,mk+2)
            endif
            if (a0.lt.0..or.a1.lt.0..or.a2.lt.0..or.a3.lt.0.) then
              g(i,j) = -1.
            else
              g(i,j) = (a0*(1.-dif) + a2*dif) * (1.-dkf)
     *               + (a1*(1.-dif) + a3*dif) * dkf
            endif
          endif
   10 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, mszxn,mszyn, mxn,myn, 9999.9, -1.
      do 100 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn)
  100 continue
      close(2)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a,$)') 'GTOPO :'
        write(99,'(a)')   ' generate topography grid-data from DEM'
        write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
        write(99,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *                 xl, yl, mszxn, mszyn, mxn, myn
        call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)    ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my'/
     *       3x, 'gen.topo grid', 2f10.2, 2i8, 2i6)
      end
