      program EXTREG9
!
! Extracts selected regions from ICOOL file  FOR009.DAT
! Writes new FOR009 file only containing selected regions
!
      parameter ( mxgood = 100 )
      real*8 tp,xp(3),pp(3),bfld(3),evtwt,efld(3),sarc,pol(3)
	integer igd(mxgood),npar(mxgood)
	logical header
	character*80 line,fmt
	data igd/mxgood*0/,npar/mxgood*0/,idec/0/,iexp/0/
!
      version = 1.06
	write(*,'(a,f6.2)') 'Welcome to program EXTREG9, version ',version
!
      open(unit=1,file='extreg9.inp',status='old',iostat=ioc)
	if(ioc .ne. 0) then
	  write(*,*)'Couldnt open EXTREG9.INP'
	  go to 900
	end if
!
      open(unit=9,file='for009.dat',status='old',iostat=ioc)
	if(ioc .ne. 0) then
	  write(*,*)'Couldnt open FOR009.DAT'
	  go to 900
	end if
!
      open(unit=10,file='extreg9.dat',status='unknown')
!
      read(1,*,iostat=ioc) knt,header
	if(ioc .ne. 0) then
	  write(*,*)'Couldnt read  KNT,HEADER  in   EXTREG9.INP'
	  go to 900
	end if
      if( knt .gt. mxgood ) then
        write(2,'(a,2x,a,i4)') 'Too many regions',
     &    'Only using ',mxgood
      end if
!
      read(1,*) (igd(i),i=1,knt)
! ---------------------------------------------------------
! Write new file with selected regions only
!
      read(9,'(a80)') line 
      if(header) write(10,'(a80)') line          
      read(9,'(a80)') line 
      if(header) write(10,'(a80)') line 
! 
      read(9,*)   ! skip line of column labels
      read(9,'(a80)') line    ! read beginning of 1st data row as characters
      do i=1,80
      if( line(i:i) .eq. '.' ) idec = i
      if( line(i:i) .eq. 'E' ) iexp = i
      if( line(i:i).eq.' ' .and. iexp.gt.0 ) then
        if( i-iexp .lt. 5 ) then
          f9dp = iexp - idec - 1
        else
          f9dp = 17
        end if
        go to 250
      end if
      end do
250   continue
      if( f9dp .eq. 4 ) then
        if(header) write(10,254)
254      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t43,'x',t55,'y',t67,'z',t79,'Px',t91,'Py',
     &  t103,'Pz',t115,'Bx',t127,'By',t139,'Bz',t151,'wt',
     &  t163,'Ex',t175,'Ey',t187,'Ez',t199,'arclength',t211,'polX',
     &  t223,'polY',t235,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e12.4)'
      else if( f9dp .eq. 6 ) then      
        if(header) write(10,256)
256      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t45,'x',t59,'y',t73,'z',t87,'Px',t101,'Py',
     &  t115,'Pz',t129,'Bx',t145,'By',t157,'Bz',t171,'wt',
     &  t185,'Ex',t199,'Ey',t213,'Ez',t227,'arclength',t241,'polX',
     &  t255,'polY',t269,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e14.6)'
      else if( f9dp .eq. 8 ) then      
        if(header) write(10,258)
258      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t47,'x',t63,'y',t79,'z',t95,'Px',t111,'Py',
     &  t127,'Pz',t143,'Bx',t159,'By',t175,'Bz',t191,'wt',
     &  t207,'Ex',t223,'Ey',t239,'Ez',t255,'arclength',t271,'polX',
     &  t287,'polY',t303,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e16.8)'
      else if( f9dp .eq. 10 ) then      
        if(header) write(10,260)
260      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t49,'x',t67,'y',t85,'z',t103,'Px',t121,'Py',
     &  t139,'Pz',t157,'Bx',t175,'By',t193,'Bz',t211,'wt',
     &  t229,'Ex',t247,'Ey',t265,'Ez',t283,'arclength',t301,'polX',
     &  t319,'polY',t337,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e18.10)'
      else if( f9dp .eq. 12 ) then      
        if(header) write(10,262)
262      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t51,'x',t71,'y',t91,'z',t111,'Px',t131,'Py',
     &  t151,'Pz',t171,'Bx',t191,'By',t211,'Bz',t231,'wt',
     &  t251,'Ex',t271,'Ey',t291,'Ez',t311,'arclength',t331,'polX',
     &  t351,'polY',t371,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e20.12)'
      else if( f9dp .eq. 14 ) then      
        if(header) write(10,264)
264      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t53,'x',t75,'y',t97,'z',t119,'Px',t141,'Py',
     &  t163,'Pz',t185,'Bx',t207,'By',t229,'Bz',t251,'wt',
     &  t273,'Ex',t295,'Ey',t317,'Ez',t339,'arclength',t361,'polX',
     &  t383,'polY',t405,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e22.14)'
      else if( f9dp .eq. 16 ) then      
        if(header) write(10,266)
266      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t55,'x',t79,'y',t103,'z',t127,'Px',t151,'Py',
     &  t175,'Pz',t199,'Bx',t223,'By',t247,'Bz',t271,'wt',
     &  t295,'Ex',t319,'Ey',t343,'Ez',t367,'arclength',t391,'polX',
     &  t415,'polY',t439,'polZ') 
        fmt = '(i7,2i4,i5,i6,1x,1p,18e24.16)'
      else if( f9dp .eq. 17 ) then      
        write(10,267)
267      format(t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg',
     &  t31,'time',t56,'x',t81,'y',t106,'z',t131,'Px',t156,'Py',
     &  t181,'Pz',t206,'Bx',t231,'By',t256,'Bz',t281,'wt',
     &  t306,'Ex',t331,'Ey',t356,'Ez',t381,'arclength',t406,'polX',
     &  t431,'polY',t456,'polZ')
        fmt = '(i7, 2i4, i5, i6, 1x, 1p, 18e25.16E3)'
      else
        if(header) write(10,266)
        fmt = '(i7,2i4,i5,i6,1x,1p,18e24.16)'
      end if
      backspace(9)      ! rewind 1st data line
! ---------------------------------------------------------
300   continue
      read(9,*,end=350) ievt,ipnum,iptyp,ipflg,jsrg,
     &  tp,xp,pp,bfld,evtwt,efld,sarc,pol
      if( ipflg .ne. 0 ) go to 300
!
      do i=1,knt
	if( jsrg .eq. igd(i) ) then
        write(10,fmt) ievt,ipnum,iptyp,ipflg,jsrg,
     &   tp,xp,pp,bfld,evtwt,efld,sarc,pol
        npar(i) = npar(i) + 1
        go to 300     ! leave loop if you find good one
      end if
	end do
!
      go to 300
! ---------------------------------------------------------
350   continue
      do i=1,knt
      write(*,380) 'Extracted ',npar(i),' particles at region ',igd(i)
380   format(a,i6,a,i4)
      end do
!
      close(1)
      close(9)
	close(10)
!
900   continue
      end