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.07 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',recl=512) ! 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