c   program to extract a subdomain out of a sequence of fort.53 files
c
c                          rl  3/21/02
c                          cf  4/19/02 - changed 74 to 64
c                          cf  10/1/04 - changed 64 to 53 and 53
c
c   an input script file named extract_53.in is required with the following
c   information:
c
c   name of subdomain grid (e.g, SAB_ec95d.grd)
c   name of output file from this program  (e.g., SAB_ec95d_extracted_fort.53)
c   y or n depending on whether it is desired to create a renumbered subdomain grid file 
c   only if last line answered y: name of renumbered subdomain grid file (e.g, SAB_ec95d_renum.grd)
c   total number of fort.53 files to be extracted from (e.g, 2)
c   name of 1st fort.53 file to be extracted from (e.g, 1-72ts_fort.74)
c   name of 2nd fort.53 file to be extracted from (e.g, 73-144ts_fort.74)
c                        .
c                        .
c                        .
c   NOTE: all file names must be no longer than 80 characters
c
c******************************************************************************

      program extract_53

      parameter (maxnodes=100000)
	parameter (maxfreq=20)

      character*80 inputa,inputb
      character*80 infilename,outfilename,gridfilename,renumgridfilename
	character*10 name
      character*1  ans
      integer i,j,k
      integer nts,np,ts_inc,ind
      integer nne,nnp
      integer ts
      integer node
      integer nnode(maxnodes)
      integer node_array(maxnodes)
      integer elenum,idum,n1,n2,n3
      integer nfiles,ntstotal
      real    t_inc
      real    x(maxnodes),y(maxnodes),z(maxnodes)		     
      real    time,dt
      real    amp(maxnodes,maxfreq),dir(maxnodes,maxfreq)

c    Open and read 1st few lines of extract_instruction file

      open(1,file='extract_53.in')
      read(1,'(a80)') gridfilename
      read(1,'(a80)') outfilename
      read(1,*) ans
      if((ans.eq.'Y').or.(ans.eq.'y')) then
        read(1,'(a80)') renumgridfilename
        endif
      read(1,*) nfiles


c    Open and read header info from subdomain grid

      open(141,file=gridfilename)
      read(141,'(a80)') inputb
      read(141,*) nne,nnp

	
c    Open subdomain fort.53 file		

      open(153,file=outfilename)


c    Read in nodes in subdomain grid

      do i=1,nnp
        read(141,*) nnode(i),x(i),y(i),z(i)
        end do

c    Create a renumbered grid if desired

      if ((ans.eq.'y').or.(ans.eq.'Y')) then
        write(*,*) 'Creating renumbered grid'
        open(241,file=renumgridfilename)
        write(241,*) ' ec95d SAB subdomain renumbered'
        write(241,*) nne,nnp
        do i=1,nnp
          write(241,*) i,x(i),y(i),z(i),nnode(i)
          node_array(nnode(i))=i
          end do
        do i=1,nne
          read(141,*) elenum,idum,n1,n2,n3
          if(node_array(n1).eq.0) write(*,*) 'warning node missing',
     &                                        i,n1
          if(node_array(n2).eq.0) write(*,*) 'warning node missing',
     &                                        i,n2
          if(node_array(n3).eq.0) write(*,*) 'warning node missing',
     &                                        i,n3
          write(241,*) i,3,node_array(n1),node_array(n2),node_array(n3)
          end do
        close(241)
        endif

c    For each fort.53 file 
c    For each time step, read in fort.53 output, then write time info
c      and subdomain elevations into subdomain fort.53 file

      do k=1,nfiles
        read(1,'(a80)') infilename
        write(*,*) 'Processing file ',infilename
        open(53,file=infilename)
        read(53,*) num_freq
	  write(153,*) num_freq
	  do i=1,num_freq
           read(53,100) freq,factor,eq_arg,name
	     write(153,100) freq,factor,eq_arg,name
	     enddo
	  read(53,*) np
	  write(153,*) nnp
 
100   format(5x,E16.10,2x,f9.7,3x,f10.8,2x,a10)

        do i=1,np
	     read(53,*) num
	     do j=1,num_freq
              read(53,101) amp(i,j),dir(i,j)
              end do
	     enddo

101   format(3x,E15.8,4x,f8.4)

        do i=1,nnp
	     write(153,*) i
	     do j=1,num_freq
              write(153,101) amp(nnode(i),j),dir(nnode(i),j)
			end do

          enddo
        close(53)				          
        enddo

      close(1)
      close(141)
      close(153)
      stop
      end