c   program to extract a subdomain out of a sequence of fort.63 files
c
c                          rl  3/20/02
c
c   an input script file named extract_63.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.63)
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.63 files to be extracted from (e.g, 2)
c   name of 1st fort.63 file to be extracted from (e.g, 1-72ts_fort.63)
c   name of 2nd fort.63 file to be extracted from (e.g, 73-144ts_fort.63)
c                        .
c                        .
c                        .
c   NOTE: all file names must be no longer than 80 characters
c
c******************************************************************************

      program extract_63

      parameter (maxnodes=100000)

      character*80 inputa,inputb
      character*80 infilename,outfilename,gridfilename,renumgridfilename
      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    elev(maxnodes)		   

c    Open and read 1st few lines of extract_instruction file

      open(1,file='extract_63.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
      read(1,'(a80)') infilename

c    Open and read header info from all fort.63 elevation output files
c    to determine the total number of records that will be used.

      ntstotal=0
      do i=1,nfiles
        open(63,file=infilename)
        read(63,'(a80)') inputa
        read(63,*) nts,np,t_inc,ts_inc,ind
        ntstotal=ntstotal+nts
        close(63)
        end do

      close(1)
      open(1,file='extract_63.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 and write header info into subdomain fort.63 file		

      open(163,file=outfilename)
      write(163,'(a80)') inputa
      write(163,999) ntstotal,nnp,t_inc,ts_inc,ind
  999 format(i11,i11,e16.7,i6,i6)


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.63 file 
c    For each time step, read in fort.63 output, then write time info
c      and subdomain elevations into subdomain fort.63 file

      do k=1,nfiles
        read(1,'(a80)') infilename
        write(*,*) 'Processing file ',infilename
        open(63,file=infilename)
        read(63,'(a80)') inputa
        read(63,*) nts,np,t_inc,ts_inc,ind
 
        do i=1,nts
          if (i.eq.1) then
            read(63,*) time,ts
            dt = time/ts
            else
            read(63,*) time
            ts=time/dt
            endif

          do j=1,np
            read(63,*) node,elev(j)
            end do

          write(163,1000) time,ts
 1000     format(e22.10,i15)
          do j=1,nnp
            write(163,1001) j,elev(nnode(j))
 1001       format(i10,e17.8)
            end do

          enddo
        close(63)				          
        enddo

      close(1)
      close(141)
      close(163)
      stop
      end