c
c     file: /export/software/puddle0/dems/demjoin.f
c     Last Modified: 11/25/96
c
c     Written by: Glenn E. Moglen
c
c   Program combines two dems into one, either in a North-South or East-West
c   joining line.  The parameter typefile controls the assumed direction of
c   joining.  If typefile = 0 then North-South joining, East-West otherwise.
c
c      This program uses subroutines from fileutil.c, the set of subroutines
c      for I/O of matrix files.
c 
c
      Program demjoin            
      PARAMETER(IGRIDY=2404, igridx=2404)          
      integer*4 nnx, snx, enx, wnx, nny, sny, eny, wny, nx, ny
      integer*2 dem1(igridy, igridx), 
     +          dem2(igridy, igridx),
     +          demtot(igridy, igridx)
      real*4 dx, dy
      character*80 northfile, southfile, eastfile, westfile, totfile
      open (11,file = 'demjoin.in', form = 'formatted', status = 'old')
      read (11, *) typefile
c    
c  The variable typefile controls whether joining is E-W or N-S
c
c  typefile:
c             = 0 ==> North - South
c             = 1 ==> East - West
      if (typefile .eq. 0) then
         read (11, 100) northfile, southfile, totfile
         call demread (dem1, northfile, nnx, nny, igridy, dx, dy)
         call demread (dem2, southfile, snx, sny, igridy, dx, dy)          
         if (nny .ne. sny) then
            print *, 'Error: unequal E-W extents of files'
            stop
         end if
         do i = 1, snx
            do j = 1, sny
               demtot(j, i) = dem2(j, i)
            end do
         end do
         do i = 1, nnx
            do j = 1, nny
               demtot(j, i + snx) = dem1(j, i)
            end do
         end do
         ny = sny
         nx = snx + nnx
      else
         read (11, 100) eastfile, westfile, totfile
         CALL demread (dem1, eastfile, enx, eny, igridy, dx, dy)
         call demread (dem2, westfile, wnx, wny, igridy, dx, dy)
         if (enx .ne. wnx) then
            print *, 'Error: unequal N-S extents of files'
            stop
         end if
         do i = 1, wnx
            do j = 1, wny
               demtot(j, i) = dem2(j, i)
            end do
         end do
         do i = 1, enx
            do j = 1, eny
               demtot(j + wny, i) = dem1(j, i)
            end do
         end do
         ny = eny + wny
         nx = enx 
      end if
      close (11, status = 'keep')
      call demwrite (demtot, totfile, nx, ny, igridy, dx, dy)
 100  format (a80/a80/a80)
      end