real lai,lao,loi,loo
      parameter(ni=4800, nj=6000)
      parameter(lao=38.3400, lai=39.99583333333333)      
      parameter(loi=100.00416666666666, loo=102.010)

      real ulxmap, ulymap, xdim, ydim
      integer i, j, a, b, c, d, u, v, w, z
      integer*2 data(ni,nj), idata
      real height(ni, nj), x(ni), y(nj)
      character*2 cdata, ctmp
      equivalence (idata, cdata)

      ulxmap = 100.00416666666666
      ulymap = 39.99583333333333
      xdim = 0.00833333333333
      ydim = 0.00833333333333

C read DEM data

      open( unit=10, form='unformatted', access='direct', recl=ni*nj*2,
     +    file='E100N40.DEM', err=999 )

      read(unit=10, rec=1) ((data(i,j), i=1,ni), j=1,nj)

      close(10)

C  make lon. & lat.
      do i=1,ni
      x(i) = ulxmap + xdim*(i-1)
      end do
      
      do j=1,nj
      y(j) = ulymap - ydim*(j-1)
      end do


C   byte swapping

      do i=1,ni
         do j=1,nj
         idata = data(i,j)
         ctmp(1:1) = cdata(1:1)
         cdata(1:1) = cdata(2:2)
         cdata(2:2) = ctmp(1:1)
         height(i,j) = idata
         end do
      end do

	
C output
      u= (loi - ulxmap) / xdim + 1      
      v= (loo - ulxmap) / xdim + 1
      w= (ulymap - lai) / ydim + 1
      z= (ulymap - lao) / ydim + 1
C      write(*,*)u,v,w,z
      do i=u,v
         do j=w,z
      write(*,*) x(i), y(j), height(i,j)
         end do
       end do
       

      open( unit=20, form='unformatted', access='direct',
     +  recl=((v-u) + 1)*((z-w) + 1)*4, file='chE100N40.DEM', err=999 )

      write(unit=20, rec=1) ((height(i,j), i=u,v), j=w,z)

      close(20)

      stop

  999  write(*,*) 'can not open'
      end