c   program to extract a subdomain out of a sequence of fort.54 files
c
c                          rl  3/21/02
c                          cf  4/19/02 - changed 74 to 64
c                          cf  10/1/04 - changed 64 to 54
c
c   an input script file named extract_54.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.54)
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.54 files to be extracted from (e.g, 2)
c   name of 1st fort.54 file to be extracted from (e.g, 1-72ts_fort.74)
c   name of 2nd fort.54 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_54

      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    xamp(maxnodes,maxfreq),xdir(maxnodes,maxfreq)
	real    yamp(maxnodes,maxfreq),ydir(maxnodes,maxfreq)		   

c    Open and read 1st few lines of extract_instruction file

      open(1,file='extract_54.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.54 file		

      open(154,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.54 file 
c    For each time step, read in fort.64 output, then write time info
c      and subdomain elevations into subdomain fort.54 file

      do k=1,nfiles
        read(1,'(a80)') infilename
        write(*,*) 'Processing file ',infilename
        open(54,file=infilename)
        read(54,*) num_freq
	  write(154,*) num_freq
	  do i=1,num_freq
           read(54,100) freq,factor,eq_arg,name
	     write(154,100) freq,factor,eq_arg,name
	     enddo
	  read(54,*) np
	  write(154,*) nnp
 
100   format(5x,E16.10,2x,f9.7,3x,f10.8,2x,a10)

        do i=1,np
	     read(54,*) num
	     do j=1,num_freq
              read(54,101) xamp(i,j),xdir(i,j),yamp(i,j),ydir(i,j)
              end do
	     enddo

101   format(3x,E15.8,4x,f8.4,3x,E15.8,4x,f8.4)

        do i=1,nnp
	     write(154,*) i
	     do j=1,num_freq
              write(154,101) xamp(nnode(i),j),xdir(nnode(i),j),
     &                     yamp(nnode(i),j),ydir(nnode(i),j)
			end do

          enddo
        close(54)				          
        enddo

      close(1)
      close(141)
      close(154)
      stop
      end

