C****************************************************************************** 
C SELECTOR TASK:
C	wgain
C
C FILE:
C	wgain.f
C
C DESCRIPTION: 
C	Process a SAX-WFC FITS file of events containing In-Flight
C	Calibration (IFC) sources to calculate a gain table.
C
C AUTHOR:  
C	Gerrit J. Wiersma 
C
C MODIFICATION HISTORY: 
C $Log:	wgain.f,v $
c Revision 1.10  99/10/13  17:02:24  17:02:24  gerwie (Gerrit J. Wiersma)
c corrected T_SCALE keyword to 1.0
c drop producing gain record with ZERO counts
c 
c Revision 1.9  97/09/16  12:30:08  12:30:08  gerwie (Gerrit J. Wiersma)
c no longer insert TASK_ID in program version
c 
C Revision 1.8  1997/09/16 12:24:28  gerwie
C added extracting IFC regions from raw (rn0) event file
C added warning output if no IFC events found
C added header keywords ONTIME, CD/FE_MEAN and CD/FE_INSTY
C
C Revision 1.7  1997/05/28 08:48:17  gerwie
C changed comment '#' to 'c'
C
C Revision 1.6  1997/04/08 13:59:33  gerwie
C added test on empty input file
C fix bug with calculating (n * step) from start
C
c Revision 1.5  1996/12/17  09:31:51  gerwie
c commentary changes
c
c Revision 1.4  1996/12/04  13:05:08  gerwie
c revised Fe to Cd ADC peak factor
c
c Revision 1.3  1996/11/28  16:18:37  gerwie
c changed kev to pha routine
c
C
C NOTES:
C	- wgain supported in IRAF and HOST environments
C	- The last bin will may not contain a full bin-time
C	- caurrent gain factors for Cd : 2.841 / 2.958
C
C ARGUMENTS:
C      none
C
C PRIMARY LOCAL VARIABLES:
C      ifcfile	- input event file 
C      gtifile  - input GTI FITS
C      outfile  - output gain file
C      cal1file - calibration file for WFC1
C      cal2file - calibration file for WFC2
C      imask	- array holding selection mask
C      xmax	- maximam x-dimension of imask
C      ymax	- maximam y-dimension of imask
C      bintime  - time of event binning in seconds
C      blockio  - fast block io option
C      verbsty  - verbosity volume level
C
C CALLED ROUTINES:
C      subroutine grunpar - gets run parameters from environment
C      subroutine procifc - process input FITS file with IFC events
C
C****************************************************************************** 
	SUBROUTINE wgain
C
	CHARACTER*160  ifcfile, gtifile, outfile, cal1file, cal2file
	INTEGER  verbsty
	LOGICAL  blockio
	DOUBLE PRECISION  bintime
C
	INTEGER  xmax,ymax
	PARAMETER (xmax = 1023)			! 767
	PARAMETER (ymax = 1023)			! 767
	INTEGER*2 imask(0:xmax,0:ymax)
C
	CHARACTER RCS_ID*256
	DATA RCS_ID /'$Header: wgain.f,v 1.10 99/10/13 17:02:24 gerwie Exp $'/
	CHARACTER TASK_ID*128
	DATA TASK_ID /'$State: Exp $'/
C
	CHARACTER ctime*8,cdate*9,fdate*8
	COMMON /wtim/ ctime,cdate,fdate
	CHARACTER*40 taskname
	COMMON /task/ taskname
	INTEGER  i1,i2
C
	ctime = ' '
	CALL Wgtime(ctime,cdate,fdate)
	taskname = 'wgain'
	i1 = INDEX(RCS_ID,'.f,v ') + 5
	i2 = INDEX(RCS_ID(i1:),' ') + i1 - 1
	IF(i1.GT.0 .AND. i2.GT.i1) taskname= 'wgain rev_'//RCS_ID(i1:i2)
	IF(INDEX(TASK_ID,'Exp').LE.0) taskname = TASK_ID
C
	verbsty = 0
	bintime = 1D0
	blockio = .FALSE.
	ifcfile = ' '
	gtifile = ' '
	outfile = ' '
	cal1file = ' '
	cal2file = ' '
C
	CALL grunpar(ifcfile,gtifile,cal1file,cal2file,outfile,
     &	    imask,xmax,ymax,bintime,blockio,verbsty)
C
	CALL procifc(ifcfile,gtifile,cal1file,cal2file,outfile,
     &	    imask,xmax,ymax,bintime,blockio,verbsty) 
C
	RETURN
	END

C****************************************************************************** 
C SUBROUTINE:
C      grunpar
C
C DESCRIPTION: 
C      Get parameters from parameter file
C
C AUTHOR:  
C	Janice Tarrant  3/24/92
C	Gerrit Wiersma  96-iv-4
C
C MODIFICATION HISTORY:
C
C NOTES:
C       uses F77/VOS like CALLs to read parameter from .par file
C
C USAGE:
C	CALL grunpar(ifcfile,gtifile,cal1file,cal2file,gainfile,
C    &		bintime,blockio,verbsty)
C
C ARGUMENTS:
C      ifcfile  - input events FITS filename
C      gtifile  - input GTI FITS filename
C      cal1file - calibration filename for WFC1
C      cal2file - calibration filename for WFC2
C      gainfile - output gain filename
C      imask	- array holding selection mask
C      xmax	- maximam x-dimension of imask
C      ymax	- maximam y-dimension of imask
C      bintime  - time in seconds to bin events
C      blockio  - fast block i/o option
C      verbsty  - verbosity volume level
C
C PRIMARY LOCAL VARIABLES:
C	title  - comment string
C	errstr - error message 
C	status - error number
C
C CALLED ROUTINES:
C	subroutine fcerr    - echo error message to terminal
C	subroutine uclgst   - get string parameter
C	subroutine uclgsb   - get boolean parameter
C	subroutine uclgsi   - get integer parameter
C	subroutine uclgsr   - get real parameter
C	subroutine uclgsd   - get double real parameter
C	subroutine wparsavt - get string parameter
C	subroutine wparsavb - get boolean parameter
C	subroutine wparsavi - get integer parameter
C	subroutine wparsavr - get real parameter
C	subroutine wparsavd - get double real parameter
C
C****************************************************************************** 
	SUBROUTINE grunpar(ifcfile,gtifile,cal1file,cal2file,gainfile,
     &		imask,xmax,ymax,bintime,blockio,verbsty)
	    CHARACTER*(*)  ifcfile,gtifile,cal1file,cal2file,gainfile
	    LOGICAL  blockio
	    INTEGER  verbsty
	    INTEGER    xmax, ymax
	    INTEGER*2  imask(0:xmax,0:ymax)
	    DOUBLE PRECISION  bintime
C
	CHARACTER  errstr*80, title*256
	INTEGER  status
C
C  initialize variables
	status = 0
C	
C= get the name of the IFC event FITS file
	CALL uclgst('ifcfile',ifcfile,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get IFCFILE parameter'
	    GOTO 999
	ENDIF
	CALL wparsavt('ifcfile',ifcfile,status)
	IF(status.NE.0) GOTO 990
C
C= get the name of the events GTI FITS file
	CALL uclgst('gtifile',gtifile,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get GTIFILE parameter'
	    GOTO 999
	ENDIF
	CALL wparsavt('gtifile',gtifile,status)
	IF(status.NE.0) GOTO 990
C
C= get the name of the output gain file
	CALL uclgst('gainfile',gainfile,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get GAINFILE parameter'
	    GOTO 999
	ENDIF 
	CALL wparsavt('gainfile',gainfile,status)
	IF(status.NE.0) GOTO 990
C
C= get the name of the WFC1 calibration file
	CALL uclgst('calfile1',cal1file,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get CALFILE1 parameter'
	    GOTO 999
	ENDIF
	CALL wparsavt('calfile1',cal1file,status) 
	IF(status.NE.0) GOTO 990
C
C= get the name of the WFC2 calibration file
	CALL uclgst('calfile2',cal2file,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get CALFILE2 parameter'
	    GOTO 999
	ENDIF
	CALL wparsavt('calfile2',cal2file,status)
	IF(status.NE.0) GOTO 990
C	
C= get the bintime parameter
	CALL uclgsd('bintime',bintime,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get BINTIME parameter'
	    GOTO 999
	ENDIF
	CALL wparsavd('bintime',bintime,status)
	IF(status.NE.0) GOTO 990
C
C= get block io option
	CALL uclgsb('blockio',blockio,status)
	IF(status.NE.0) THEN
	    errstr = 'could not find BLOCKIO factor'
	    GOTO 999
	ENDIF
	CALL wparsavb('blockio',blockio,status)
	IF(status.NE.0) GOTO 990
C
C= get verbosity factor
	CALL uclgsi('verbsty',verbsty,status)
	IF(status.NE.0) THEN
	    errstr = 'could not find VERBSTY factor'
	    GOTO 999
	ENDIF
	CALL wparsavi('verbsty',verbsty,status)
	IF(status.NE.0) GOTO 990
C
C= get the TITLE string
	CALL uclgst('title',title,status)
	IF(status.NE.0) THEN
	    errstr = 'could not get TITLE parameter'
	    GOTO 999
	ENDIF
	CALL wparsavt('title',title,status)
	IF(status.NE.0) GOTO 990
C
C= normal end
	RETURN
C
C= error condition
  990	errstr = 'parameter buffer overflow'
C
  999	CALL fcerr('*****ERROR: '//errstr)
	STOP
C
	END

C****************************************************************************** 
C SUBROUTINE:
C      procifc
C
C DESCRIPTION: 
C      Reads the FITS file and writes the GTI filtered gaib values to
C      a FITS file
C
C AUTHOR:  
C      Gerrit J. Wiersma (SRON) 96-iv-2
C
C MODIFICATION HISTORY:
C
C NOTES:
C
C USAGE:
C      CALL procifc(ifcfile,gtifile,cal1file,cal2file,outfile,
C    &		 imask,xmax,ymax,bintime,blockio,verbsty)
C
C ARGUMENTS:
C      ifcfile  - input events FITS file
C      gtifile  - input GTI FITS file
C      cal1file - calibration file for WFC1
C      cal2file - calibration file for WFC2
C      outfile  - output file name
C      imask	- array holding selection mask
C      xmax	- maximam x-dimension of imask
C      ymax	- maximam y-dimension of imask
C      bintime  - time in seconds to bin eventtimes
C      blockio  - fast block io option
C      verbsty  - verbosity volume level
C
C PRIMARY LOCAL VARIABLES:
C      maxcl    - maximum number of columns
C      filename - name of FITS file
C      ftstat   - FITSIO error number
C      errstr	- error message
C      context  - general message 
C      tcolnum	- column number for time field
C      xcolnum	- column number for x-position field
C      ycolnum	- column number for y-position field
C      zcolnum	- column number for z-position field
C      negflag  - exclude name flag 
C      obstime  - observation start time keyword
C      odate    - observation start date
C      otime    - observation start time
C      toffs    - observation time offset
C      tscale   - observation time unit
C      redfec	- reduction factor
C      regs	- array holding coordinates of selection regions
C      nreg	- number of regions used 
C
C
C CALLED ROUTINES:
C	subroutine evtproc
C	subroutine total_gti
C	subroutine getcalrec  
C	subroutine getadckev
C	subroutine getifcregs
C	subroutine wparfits
C      subroutine fcecho - echo message to terminal
C      subroutine fcerrm - report an error number to terminal
C      subroutine fcgcls - get column list based on parameter
C      subroutine fcpars - parse off filename and extension number
C      function   fcstln - getlength of string
C      subroutine ftadef - define ASCII extension
C      subroutine ftbdef - define binary extension
C      subroutine ftclos - close a FITS file
C      subroutine ftcrhd - create a new header
C      subroutine ftgbnh - get BINARY table header
C      subroutine ftgcno - get column number associated with name
C      subroutine ftgkys - get string keyword value
C      subroutine ftgprh - get primary header keywords from CHU
C      subroutine ftgtbh - get ASCII table header
C      subroutine ftinit - create a new FITS file
C      subroutine ftmahd - absolute move to FITS header
C      subroutine ftmrhd - relative move to FITS header
C      subroutine ftopen - open a FITS file
C      subroutine ftpbnh - put binary table header keywords into CHU
C      subroutine ftpdef - define structure of primary array
C      subroutine ftphis - put history keyword record
C      subroutine ftpprh - put primary header keywords into CHU
C      subroutine ftptbh - put ASCII table header keywords into CHU
C
C******************************************************************************
C 
	SUBROUTINE procifc(ifcfile,gtifile,cal1file,cal2file,outfile,
     &		imask,xmax,ymax,bintime,blockio,verbsty)
	    CHARACTER*(*)  ifcfile,gtifile,cal1file,cal2file,outfile
	    LOGICAL  blockio
	    INTEGER  verbsty
	    INTEGER    xmax, ymax
	    INTEGER*2  imask(0:xmax,0:ymax)
	    DOUBLE PRECISION  bintime
C
	INTEGER maxcl
	PARAMETER (maxcl = 512)
	CHARACTER*160 filename
	CHARACTER*80  context, errstr, tcolnam
	CHARACTER*70  ttype(maxcl), tform(maxcl), tunit(maxcl), extname,
     &	              comment, detcalrev
	LOGICAL  iopen, gopen, oopen, exact, negflag, simple, extend
	INTEGER  ftstat, iunit, gunit, ounit, block, extnum, htype, 
     &           gcount, irows, tfields, varidat,
     &           bitpix, naxis, naxes(99), pcount, mkywd
	INTEGER  tcolnum, cunit, camera, next, redfac, ntotal, nrows, gorow
	DOUBLE PRECISION  toffs, tscale, obstime, mjdmin, mjdmax, totalt
	REAL  adcmean(33),feadcnom,cdadcnom,feadclo,feadchi,cdadclo,cdadchi
C
	REAL  femean,cdmean,feinty,cdinty,cmean
	INTEGER  caccu,i,enomfe,enomcd
	
	INTEGER  regs(4,50),nreg,x,y
	LOGICAL  rawevents
C
	INTEGER  ncdlo, ncdhi, nfelo, nfehi
	DOUBLE PRECISION  chantot,ontime
	COMMON /ifcstat/ ncdlo,ncdhi,nfelo,nfehi,chantot(32),ontime
C
	CHARACTER  ctime*8,cdate*9,fdate*8
	COMMON /wtim/ ctime,cdate,fdate
C
	CHARACTER*40  taskname
	COMMON /task/ taskname
C
C  initialize variables
	ftstat = 0
	iunit = 15 
	gunit = 16
        ounit = 17
	cunit = 18
	iopen = .FALSE.
	gopen = .FALSE.
	oopen = .FALSE.
	negflag = .FALSE.
	exact = .TRUE.
	tscale = 1D0
	mjdmax = 0D0
	mjdmin = 999999D0
	ntotal = 0
	gorow  = 0
C
	nfelo =  8
	nfehi = 20
	ncdlo = 25
	ncdhi = 32
	enomfe =  5.98		!  13.58 chan
	enomcd = 22.60		!  29.55 chan
	DO i = 1,32
	  chantot(i) = 0
	ENDDO
	ontime = 0D0
C
C= get the filename and extension number
	CALL fcpars(ifcfile,filename,extnum,ftstat)
C
C= open the event input FITS file and check contents
C
	CALL FTopen(iunit,filename,0,block,ftstat)
	IF(ftstat.NE.0) THEN
	  errstr = 'unable to open ifcfile'
	  GOTO 990
	ENDIF
	CALL ftgprh(iunit,simple,bitpix,naxis,naxes,pcount,
     &	            gcount,extend,ftstat)
	iopen = .TRUE.
C
	CALL FTgkys(iunit,'CONTENT',context,comment,ftstat)
	rawevents = .TRUE.
	IF(context.EQ.'SCI_IFC') rawevents = .FALSE.
	IF(context.NE.'SCI_NM' .AND. rawevents) THEN
	  errstr = 'illegal input file type'
	  GOTO 990
	ENDIF
C
	camera = 0
	CALL FTgkys(iunit,'INSTRUME',context,comment,ftstat)
	IF(context.EQ.'WFC1') camera = 1
	IF(context.EQ.'WFC2') camera = 2
	IF(camera.LE.0) THEN
	    errstr = 'bad camera specified'
	    GOTO 990
	ENDIF
C
	CALL ftmahd(iunit,2,htype,ftstat)
	IF(ftstat.NE.0) THEN
	    errstr = 'Cannot move to extension header'
	    GOTO 990
	ENDIF	
	obstime = -1D0
	CALL FTgkyd(iunit,'MJD_OBS',obstime,comment,ftstat) 
	IF(ftstat.NE.0) THEN
	   ftstat = 0
	   CALL FTgkyd(iunit,'MJD-OBS',obstime,comment,ftstat) 
	ENDIF
	IF(ftstat.NE.0) THEN
	    errstr = 'No MJD_OBS specified'
	    GOTO 990
	ENDIF
	CALL ftmahd(iunit,1,htype,ftstat)
	IF(ftstat.NE.0) THEN
	    errstr = 'Cannot move to primary header'
	    GOTO 990
	ENDIF	
C
C= read from the calibration FITS file 
	IF(camera.EQ.1) THEN 
	  CALL getcalrev(cunit,cal1file,detcalrev,ftstat)
	  IF(ftstat.NE.0) THEN
	    errstr = 'Cannot open detcalfile wfc1'
	    GOTO 990
	  ENDIF	
	  CALL getadckev(cunit,cal1file,obstime,adcmean,feadcnom,cdadcnom,
     &			feadclo,feadchi,cdadclo,cdadchi,ftstat) 
          IF(cdadcnom.LT.feadcnom) cdadcnom = feadcnom * 2.841
	ELSE 
	  CALL getcalrev(cunit,cal2file,detcalrev,ftstat)
	  IF(ftstat.NE.0) THEN
	    errstr = 'Cannot open detcalfile wfc2'
	    GOTO 990
	  ENDIF	
	  CALL getadckev(cunit,cal2file,obstime,adcmean,feadcnom,cdadcnom,
     &			feadclo,feadchi,cdadclo,cdadchi,ftstat)
          IF(cdadcnom.LT.feadcnom) cdadcnom = feadcnom * 2.958
	ENDIF
	IF(ftstat.NE.0) THEN
	  errstr = 'channel conversion table not found'
	  GOTO 990
	ENDIF
C	
C= init IFC mask array (0 = accept, -1 = reject , +1 = store seperate)
C
	DO y= 0,ymax
	  DO x= 0,xmax
	    IF(x.LT.768 .AND. y.LT.768) THEN
	      IF(rawevents) THEN
	        imask(x,y) = 0
	      ELSE
	        imask(x,y) = 99
	      ENDIF
	    ELSE
	      imask(x,y) = -1
	    ENDIF
	  ENDDO
	ENDDO
C
C- read IFC regions if not filtered eventfile
C
	IF(rawevents) THEN
	  IF(camera.EQ.1)THEN
	    CALL getifcregs(cunit,cal1file,obstime,regs,nreg,ftstat)
	  ELSE
	    CALL getifcregs(cunit,cal2file,obstime,regs,nreg,ftstat)
	  ENDIF
	  IF(ftstat.NE.0) THEN
	    errstr = 'CAL_SOURCES regions not found'
	    GOTO 990
	  ENDIF 
C
	  DO i=1,nreg
	    DO y= regs(2,i),regs(4,i)
	      DO x= regs(1,i),regs(3,i)
	        imask(x,y) = nreg
	      ENDDO
	    ENDDO
	  ENDDO 
	ENDIF
C
C= open the GTI input FITS file
C
	IF(gtifile.NE.' ' .AND. gtifile.NE.'-') THEN
	  CALL FTopen(gunit,gtifile,0,block,ftstat)
	  IF(ftstat.NE.0) THEN
	    errstr = 'unable to open gtifile' 
	    GOTO 990
	  ENDIF
c--	  CALL ftgprh(gunit,simple,bitpix,naxis,naxes,pcount,
c--     &	              gcount,extend,ftstat) 
	  CALL ftmahd(gunit,2,htype,ftstat)
	  IF(ftstat.NE.0) THEN
	    errstr = 'incorrect gtifile' 
	    GOTO 990
	  ENDIF
	  CALL FTgkyj(gunit,'NAXIS2',nrows,comment,ftstat)
	  IF(ftstat.NE.0 .OR. nrows.LE.0) THEN
	    errstr = 'incorrect or empty gtifile' 
	    GOTO 990
	  ENDIF
	  gopen = .TRUE.
	ENDIF
C
C= open the new event output FITS file
	CALL ftinit(ounit,outfile,1,ftstat)
	IF(ftstat .NE. 0) THEN
	  errstr = 'unable to open outfile' 
	  GOTO 990
	ENDIF
	oopen = .TRUE.
C	
C= copy and modify primary header for the new file 
C
	CALL ftmahd(iunit,1,htype,ftstat)
	mkywd = 0
	CALL ftcopy(iunit,ounit,mkywd,ftstat)
	CALL fthdef(ounit,15,ftstat)
C
	CALL ftmkys(ounit,'DATE',fdate,'&',ftstat) 
	IF(ftstat.NE.0) THEN
	  errstr = 'keyword DATE not found'
	  GOTO 990
	ENDIF
C
	CALL FTmkys(ounit,'DATE'    , fdate,'&',ftstat)
	CALL FTmkys(ounit,'CREATOR' , taskname,'&',ftstat)
	CALL FTmkys(ounit,'CONTENT' , 'GAINS','&',ftstat) 
	CALL FTmkyj(ounit,'REVISION', 2,'&',ftstat)  
	CALL FTmkyj(ounit,'WFC_EXTS', 1,'&',ftstat) 
	IF(ftstat.NE.0) THEN
	  errstr = 'standard keywords not found'
	  GOTO 990
	ENDIF
C
	IF(camera.EQ.1) THEN
	  CALL FTmkys(ounit,'CAL_FILE',cal1file,'&',ftstat) 
	  IF(ftstat.NE.0) THEN
	    ftstat = 0	 
	    CALL FTpkls(ounit,'CAL_FILE',cal1file,'Name of det. cal. file',ftstat)  
	  ENDIF
	ELSE
	  CALL FTmkys(ounit,'CAL_FILE',cal2file,'&',ftstat)  
	  IF(ftstat.NE.0) THEN
	    ftstat = 0	 
	    CALL FTpkls(ounit,'CAL_FILE',cal2file,'Name of det. cal. file',ftstat)
	  ENDIF
	ENDIF
	IF(ftstat.NE.0) THEN
	  errstr = 'cannot update CAL_FILE keyword'
	  GOTO 990
	ENDIF
	CALL FTmkys(ounit,'DETREVN',detcalrev,'&',ftstat) 
	IF(ftstat.NE.0) THEN
	  ftstat = 0	 
	  CALL FTpkys(ounit,'DETREVN',detcalrev,'Rev. of det. cal. file',ftstat)  
	ENDIF
	IF(ftstat.NE.0) THEN
	  errstr = 'cannot update DETREVN keyword'
	  GOTO 990
	ENDIF
C
	CALL wparfits(ounit,taskname,ftstat)
	IF(ftstat.NE.0) THEN
	  errstr = 'unable to save run parameters'
	  GOTO 990
	ENDIF
	CALL ftrdef(ounit,ftstat)
	IF(ftstat.NE.0) THEN
	  errstr = 'unable to copy header'
	  GOTO 990
	ENDIF
C
	CALL fonew(ounit,bintime,ftstat)
	IF(ftstat.NE.0) GOTO 999
C
C= create a new extension
C
	next = 1
  500	next = next + 1
C
C= move to the input extension number or stop processing if not available
C
	    CALL ftmahd(iunit,next,htype,ftstat)
	    IF(ftstat.NE.0) GOTO 600
C
C= get the binary extension header 
C
	    CALL ftgbnh(iunit,irows,tfields,ttype,tform,tunit,
     &	                extname,varidat,ftstat)
	    IF(ftstat.NE.0) THEN
		errstr = 'no time column found'
		GOTO 990
	    ENDIF

C= check that column name exists and get its position in the input file
C
	    tcolnam = 'TIME'
	    IF(ttype(1).EQ.tcolnam) THEN
	      CALL ftgcno(iunit,exact,tcolnam,tcolnum,ftstat)
	    ELSE
	      tcolnam = 'TIME_WFC'
	      IF(ttype(1).EQ.tcolnam) THEN
	        CALL ftgcno(iunit,exact,tcolnam,tcolnum,ftstat)
	      ELSE
		errstr = 'no time column found'
		GOTO 990
	      ENDIF
	    ENDIF
	    IF(ftstat.NE.0) THEN
		errstr = 'bad time column'
		GOTO 990
	    ENDIF
C
C  get the time scale  
	    CALL ftgkyd(iunit,'T_SCALE',tscale,comment,ftstat) 
	    IF(ftstat .NE. 0) THEN 
		errstr = 'cannot get time unit'
	  	GOTO 990
	    ENDIF 
	    tscale = tscale / 86400D0
	    
C  get the time offset
	    CALL ftgkyd(iunit,'MJDREF',toffs,comment,ftstat)
	    IF(ftstat.NE.0) THEN
	        toffs = 0.D0 
		ftstat = 0
	    ENDIF
C
C  get the reduction factor 
	    CALL ftgkyj(iunit,'RED_FAC',redfac,comment,ftstat) 
	    IF(ftstat.NE.0) THEN 
		errstr = 'cannot get reduction factor'
	  	GOTO 990
	    ENDIF 
C==
C= process events
C
	    CALL evtproc(iunit,gunit,ounit,blockio,gopen,irows,redfac,
     &		tscale,toffs,camera,bintime,verbsty,imask,xmax,ymax,
     &		adcmean,feadcnom,cdadcnom,feadclo,feadchi,cdadclo,cdadchi,
     &		mjdmin,mjdmax,ntotal,gorow,ftstat) 
	    IF(ftstat.NE.0) GOTO 999
C
C= calculate effective ontime
C
	    IF(gopen) THEN 
	      CALL total_gti(gunit,mjdmin,mjdmax,totalt,verbsty)
	    ELSE
	      totalt = mjdmax - mjdmin
	    ENDIF
	    
	    IF(verbsty.GE.2) WRITE(*,*) (mjdmax-mjdmin)*86400D0, totalt*86400D0	    
	    ontime = ontime + totalt * 86400D0	 
C
C= goto next extension
C
	GOTO 500
C===
C= all extensions processed / modify dates in headers
C
C= calculate mean and intensity of total events
C
  600	ftstat = 0
C
C= correct the number of rows keyword (width is same for output and ifc)
C
	IF(gorow.GT.0) THEN
	  IF(verbsty.GT.0) WRITE(*,*) 'total gain values: ',gorow,' from ',ntotal
	  CALL ftmkyj(ounit,'NAXIS2',gorow,'&',ftstat)
	  CALL ftrdef(ounit, ftstat)
	  IF(ftstat.NE.0) THEN 
	    errstr = 'cannot update keyword'
	    RETURN
	  ENDIF
	ELSE
	  WRITE(*,*) '***WARNING: ZERO GAIN values'
	ENDIF
C
C= store statistics
C
	caccu = 0
	cmean = 0.0
	DO i = nfelo,nfehi
	  caccu = caccu + chantot(i)
	  cmean = cmean + chantot(i) * i
	ENDDO
	femean = cmean / caccu
	feinty = caccu / ontime	
	IF(verbsty.GE.2) WRITE(*,*) 'Fe: ',caccu,cmean,ontime	
	caccu = 0
	cmean = 0.0
	DO i = ncdlo,ncdhi
	  caccu = caccu + chantot(i)
	  cmean = cmean + chantot(i) * i
	ENDDO
	cdmean = cmean / caccu
	cdinty = caccu / ontime	
	IF(verbsty.GE.2) WRITE(*,*) 'Cd: ',caccu,cmean,ontime
C
	CALL FTmkyd(ounit,'ONTIME'  ,ontime,12,'&',ftstat) 
	CALL FTmkye(ounit,'CD_MEAN' ,cdmean,6,'&',ftstat)
	CALL FTmkye(ounit,'FE_MEAN' ,femean,6,'&',ftstat)
	CALL FTmkye(ounit,'CD_INSTY',cdinty,6,'&',ftstat)
	CALL FTmkye(ounit,'FE_INSTY',feinty,6,'&',ftstat)
	IF(ftstat.NE.0) THEN
	  errstr = 'cannot update keywords'
	  GOTO 990
	ENDIF
	IF(verbsty.GE.2) WRITE(*,*) mjdmin,mjdmax,(mjdmax-mjdmin)*86400D0	
C
	ftstat = 0
  	CALL ftmkyd(ounit,'MJD_OBS',mjdmin,12,'Obs. start MJD',ftstat)
	IF(ftstat.NE.0) THEN
	  ftstat = 0
  	  CALL ftpkyd(ounit,'MJD_OBS',mjdmin,12,'Obs. start MJD',ftstat)
	ENDIF
	IF(ftstat.NE.0) THEN 
	  errstr = 'cannot update keyword in header'
	  GOTO 990
	ENDIF 
C
  	CALL ftmkyd(ounit,'MJD_END',mjdmax,12,'Obs. end MJD',ftstat)
	IF(ftstat.NE.0) THEN
	   ftstat = 0
  	  CALL ftpkyd(ounit,'MJD_END',mjdmax,12,'Obs. start MJD',ftstat)
	ENDIF
	IF(ftstat.NE.0) THEN 
	  errstr = 'cannot update keyword in header'
	  GOTO 990
	ENDIF
C
	CALL ftpkyj(ounit,'E_TOTAL',ntotal,'total accepted events',ftstat)
	IF(ftstat.NE.0) THEN 
	  errstr = 'cannot add keyword to header'
	  GOTO 990
	ENDIF
C 
  	CALL upd_mjdobs(ounit,mjdmin,mjdmax,ftstat)
C
C= close files and return on error
C
  990	IF(ftstat.NE.0) CALL fcerr('*****ERROR: '//errstr)
C
  999	IF(ftstat.NE.0) CALL fcerrm(ftstat)
	ftstat = 0
	IF(iopen)  CALL ftclos(iunit,ftstat)
	ftstat = 0
	IF(gopen)  CALL ftclos(gunit,ftstat)
	ftstat = 0
	IF(oopen)  CALL ftclos(ounit,ftstat)
C
 1000	FORMAT(A34,I3)
	RETURN
	END
C******************************************************************************
C SUBROUTINE:
C      fonew
C
C DESCRIPTION:
C      create new extension for output gain file
C
C AUTHOR/DATE:
C      Gerrit J. Wiersma
C
C MODIFICATION HISTORY:
C
C NOTES:
C
C USAGE:
C      CALL fonew(outlu,ftstat) 
C
C ARGUMENTS:
C	outlu	- pointing file LOGICAL unit number
C	bintime - length of the bin in seconds
C	ftstat	- return I/O status on FITS file
C
C PRIMARY LOCAL VARIABLES:
C	tfield	- number of fields (columns) in a record (row)
C
C CALLED ROUTINES:
C      subroutine ftcrhd - create extension header
C      subroutine ftpbnh - initiate binary extension header
C      subroutine fcerr  - oiutput error message
C
C******************************************************************************
	SUBROUTINE fonew(outlu,bintime,ftstat)
		INTEGER  outlu,ftstat 
		DOUBLE PRECISION  bintime
C
	INTEGER  pcount, tfield
	CHARACTER*16 ttype(9),tform(9),tunit(9),extnam
C
C- define columns for FITS extended header
C 
	tfield= 8		! number of fields 
C
	ttype(1)='TIME'   
	tform(1)='1D'	! r8 
	tunit(1)='MJD'
C
	ttype(2)='FE_GAIN'  
	tform(2)='1E'	! r4
	tunit(2)='factor'
C
	ttype(3)='CD_GAIN'  
	tform(3)='1E'	! r4
	tunit(3)='factor'
C
	ttype(4)='FE_MEAN'  
	tform(4)='1E'	! r4 
	tunit(4)='adc_chan'
C
	ttype(5)='CD_MEAN'   
	tform(5)='1E'	! r4
	tunit(5)='adc_chan'
C
	ttype(6)='FE_EFF'  
	tform(6)='1E'	! r4 
	tunit(6)='factor'
C
	ttype(7)='CD_EFF'   
	tform(7)='1E'	! r4
	tunit(7)='factor'
C
	ttype(8)='COUNTS'   
	tform(8)='1J'	! i4
	tunit(8)='number'
C
C= row-length is 8 + 7 * 4 = 36 bytes ( 80 * 36 = 2880 )
C
	ftstat = 0
C
C= init new extension
C 
	extnam = 'GAINS'
	pcount= 0
  	CALL FTcrhd(outlu,ftstat) 
	CALL FTpbnh(outlu,0,tfield,ttype,tform,tunit,extnam,pcount,ftstat) 
	IF(ftstat.NE.0) THEN
	  CALL fcerr('*****ERROR cannot create extension header ')
	  RETURN
	ENDIF
C
C= write some keywords 
C
	CALL FThdef(outlu,15,ftstat) 
	CALL FTpkyd(outlu,'TBIN_GN' ,bintime,12,'gain binning time in s.',ftstat)
	CALL FTpkyd(outlu,'MJDREF'  ,0D0,12,'reference date',ftstat)
	CALL FTpkyd(outlu,'T_SCALE' ,86400D0,12,'time scale',ftstat)
	CALL FTpkyd(outlu,'MJD_OBS' ,0D0,12,'start of obs. (MJD)',ftstat)
	CALL FTpkyd(outlu,'MJD_END' ,0D0,12,'end of obs. (MJD)',ftstat)
	CALL FTpkyd(outlu,'ONTIME'  ,0D0,12,'accumulated seconds',ftstat) 
	CALL FTpkye(outlu,'CD_MEAN' ,0.0,6,'mean channel',ftstat)
	CALL FTpkye(outlu,'FE_MEAN' ,0.0,6,'mean channel',ftstat)
	CALL FTpkye(outlu,'CD_INSTY',0.0,6,'intensity cts/sec.',ftstat)
	CALL FTpkye(outlu,'FE_INSTY',0.0,6,'intensity cts/sec.',ftstat)
	IF(ftstat.NE.0) THEN
	  CALL fcerr('*****ERROR: cannot add keywords')
	  RETURN
	ENDIF
	CALL ftbdef(outlu,tfield,tform,pcount,0,ftstat) 
C
	RETURN
	END

C*****************************************************************************
C SUBROUTINE:
C      evtproc
C
C DESCRIPTION:
C      read IFC events and store in arrays
C
C AUTHOR/DATE:
C      Gerrit J. Wiersma
C
C MODIFICATION HISTORY:
C
C NOTES:
C	- It is assumed that the input file is time ordered
C	- 
C
C USAGE:
C      CALL evtproc(iunit,gunit,ounit,blockio,gopen,irows,redfac,
C		    tscale,toffs,camera,bintime,verbsty,ftstat) 
C
C ARGUMENTS:
C	iunit	- input event file LOGICAL unit number
C	gunit	- input GTI file LOGICAL unit number
C	ounit	- output gain file LOGICAL unit number
C	blockio	- fast blockio reading of file
C	gopen	- use GTI file for time selection
C	irows	- maximum number of rows of input file (extension)
C	toffs	- offset of times in MJD
C	tscale	- scale factor of times units in seconds 
C	verbsty - output verbosity volume
C	bintime	- binning time in seconds for one gain calculation 
C	ntotal  - total number of accepted events
C	gorow	- output row number
C	mjdmin	- time of earliest accepted event
C	mjdmax	- time of latest accepted event
C
C PRIMARY LOCAL VARIABLES:
C	width	- width of table 
C	frow	- beginning row number
C	remain	- number of rows remaining
C	nelem	- number of data elements
C	time 	- observation times
C	time1	- start time of bin
C	time2	- stop time of bin
C	tmid	- middle time of timebin
C
C	esum	- array holding summed energy bins
C
C CALLED ROUTINES:
C	subrouitne anafun - analyse the function
C	
C	subroutine ftgcfe - get real column values
C	subroutine ftgkyj - get INTEGER keyword value
C	subroutine ftmkyj - modify INTEGER keyword value
C
C******************************************************************************
C
	SUBROUTINE evtproc(iunit,gunit,ounit,blockio,gopen,irows,redfac,
     &		  tscale,toffs,camera,bintime,verbsty,imask,xmax,ymax,
     &		  adcmean,feadcnom,cdadcnom,feadclo,feadchi,cdadclo,cdadchi,
     &		  mjdmin,mjdmax,ntotal,gorow,ftstat)
C
	    INTEGER  iunit, gunit, ounit, irows, camera, verbsty
	    INTEGER  redfac, ntotal, gorow
	    DOUBLE PRECISION  bintime, tscale, toffs, mjdmin, mjdmax
	    REAL adcmean(*),feadcnom,cdadcnom,feadclo,feadchi,cdadclo,cdadchi
	    LOGICAL  blockio, gopen
	    INTEGER    xmax, ymax
	    INTEGER*2  imask(0:xmax,0:ymax)
C
	INTEGER  BUFMAX
	PARAMETER (BUFMAX = 1000)
C
	CHARACTER  comment*70
	LOGICAL  anyf, flagvals(BUFMAX)
	INTEGER  width, frow, felem, remain, nelem, ftstat
	INTEGER  i, n, nbin, nadc, nstep
c--	INTEGER  tcolnum,xcolnum,ycolnum,zcolnum
	DOUBLE PRECISION  t,tev,time1,time2,tmid,tstep
	REAL sumfac
	LOGICAL inside_gti, binstart, inside, latest
C
	INTEGER*4  time(BUFMAX) 
	INTEGER*2  eval(BUFMAX), xval(BUFMAX), yval(BUFMAX)
	INTEGER*4  buf4(3,BUFMAX)
	INTEGER*2  buf2(6,BUFMAX) 
	EQUIVALENCE (buf4,buf2)
C
	REAL  esum(32), adc
C
	INTEGER x, y
	REAL  e, zcd, zfe, efffe, effcd
	REAL  cpeak, cmean, cmedn, sigma
	REAL  fwhm, cfwhm, ew50, ew68, ew90
	REAL  cmeancd, cmeanfe, gaincd, gainfe
	INTEGER  ipeak
C
	INTEGER ncdlo, ncdhi, nfelo, nfehi
	DOUBLE PRECISION chantot,ontime
	COMMON /ifcstat/ ncdlo,ncdhi,nfelo,nfehi,chantot(32),ontime
C=
C== end of declarations
C===
	binstart = .FALSE.
	nstep = 1
	tstep = bintime / 86400D0
	sumfac = redfac + 1
C
C- reset sum arrays
C 
	DO n = 1,32
	  esum(n) = 0.
	ENDDO
C
C= get the input row width
C
	ftstat=0
	CALL ftgkyj(iunit,'NAXIS1',width,comment,ftstat)
	IF(ftstat.NE.0) THEN
	  CALL fcerr('*****ERROR: cannot read keyword from event input file')
	  RETURN
	ENDIF
	IF(verbsty.GE.2) WRITE(*,*) 'procevt ',bintime,irows,width
C
C= get (next) buffer of input file records
C
	latest = .FALSE.
	nbin = 0
	frow = 1
	felem = 1
	remain = irows
 10	IF(remain .GE. BUFMAX) THEN
	  nelem = BUFMAX
	ELSE
	  nelem = remain
	ENDIF
C
C= read buffer of BUFMAX events
C
	IF(blockio) THEN 
	  CALL FTgtbb(iunit,frow,1,width*nelem,buf2,ftstat)
	ELSE
	  CALL ftgcfj(iunit,1,frow,felem,nelem,time,flagvals,anyf,ftstat) 
	  CALL ftgcfi(iunit,2,frow,felem,nelem,xval,flagvals,anyf,ftstat) 
	  CALL ftgcfi(iunit,3,frow,felem,nelem,yval,flagvals,anyf,ftstat) 
c--	  CALL ftgcfi(iunit,4,frow,felem,nelem,zval,flagvals,anyf,ftstat) 
	  CALL ftgcfi(iunit,5,frow,felem,nelem,eval,flagvals,anyf,ftstat)
	ENDIF
	IF(ftstat.NE.0) THEN
	  CALL fcerr('*****ERROR: cannot read events from input file')
	  RETURN
	ENDIF
C==
C= process each single event in buffer
C
	DO 50 i = 1, nelem
	  IF(remain.EQ.i) latest = .TRUE.
C
C= get next event time and energy channel
C
	  IF(blockio) THEN
	    t = buf4(1,i)
	    x = buf2(3,i)
	    y = buf2(4,i)
	    e = buf2(6,i) + 1
	  ELSE
	    t = time(i)
	    x = xval(i)
	    y = yval(i)
	    e = eval(i) + 1
	  ENDIF 	  
	  IF(e.GT.32) STOP '*****ERROR in wgain: E .GT. 31'   
C
C- test if event inside time window
C
	  inside = .TRUE.
	  IF(t.LE.0D0) THEN 
	    inside = .FALSE.
	  ELSE
	    tev = ( t * tscale ) + toffs
	    IF(gopen) inside = inside_gti(gunit,tev,verbsty)
	  ENDIF
C
C= skip events until begin timebin with first correct event time
C
	  IF(.NOT.binstart) THEN 
	    IF(inside) THEN
	      time1 = tev
	      time2 = time1 + tstep
	      binstart = .TRUE.
	    ELSE
	      GOTO 50
	    ENDIF
	  ENDIF	  
C
C- test if inside region
C
	  IF(inside) THEN
	    IF(imask(x,y).LE.0) inside = .FALSE.
	  ENDIF
	  IF(inside) THEN
	    IF(tev.GT.mjdmax) mjdmax = tev
	    IF(tev.LT.mjdmin) mjdmin = tev
	    chantot(e) = chantot(e) + sumfac
	    ntotal= ntotal + 1
	  ENDIF
C
C= set end of time bin if last event
C= test whether inside Good Time Interval
C
C if tev G.E. time2 dan is that de eerste voor de volgende bin
C
C time2 = end time of bin
C
	  IF(latest) THEN
	    IF(inside) THEN
	      nbin = nbin + 1
	      esum(e) = esum(e) + sumfac
	    ENDIF
	    tmid = ( time2 + tev - tstep ) * 0.5D0
	  ELSE
	    IF(.NOT.inside) GOTO 50
	    IF(tev.LT.time2) THEN
	      nbin = nbin + 1
	      esum(e) = esum(e) + sumfac
	      GOTO 50
	    ENDIF
	    tmid =  time2 - ( tstep * 0.5D0 )
	  ENDIF
C
C= now process complete energy bin
C
	    IF(verbsty.GE.3) THEN 
	      DO n = 1,32
	        WRITE(*,*) n-1,esum(n)
	      ENDDO
	    ENDIF
	    IF(nbin.GT.1) THEN
C
C= analyze the two different energy peaks its position, mean and fwhm
C
	      CALL anafun(esum,nfelo,nfehi,zfe,ipeak,cpeak,cmean,cmedn,sigma,
     &			fwhm,cfwhm,ew50,ew68,ew90) 
	      nadc = 0
	      cmeanfe = 0.0
	      DO n = 1,32
	        adc = adcmean(n)
	        IF(adc.GT.feadclo .AND. adc.LT.feadchi) THEN
	          nadc = nadc + esum(n)
	          cmeanfe = cmeanfe + adc * esum(n)
	        ENDIF
	      ENDDO
	    
	      IF(nadc.GT.0) THEN 
	        cmeanfe = cmeanfe / nadc
	        gainfe = cmeanfe / feadcnom
	      ELSE
	        IF(verbsty.GT.0) WRITE(*,*) '*WARNING: zero Fe coints'
	        cmeanfe = -1.0
	      ENDIF
	      IF(verbsty.GE.2) WRITE(*,*) nadc,cmeanfe,cmean-1.0,fwhm
C
	      CALL anafun(esum,ncdlo,ncdhi,zcd,ipeak,cpeak,cmean,cmedn,sigma,
     &			fwhm,cfwhm,ew50,ew68,ew90)
	      nadc = 0
	      cmeancd = 0.0
	      DO n = 1,32
	        adc = adcmean(n)
	        IF(adc.GT.cdadclo .AND. adc.LT.cdadchi) THEN
	          nadc = nadc + esum(n)
	          cmeancd = cmeancd + adc * esum(n)
	        ENDIF
	      ENDDO
		
	      IF(nadc.GT.0) THEN 
	        cmeancd = cmeancd / nadc
	        gaincd  = cmeancd / cdadcnom 
	      ELSE
	        IF(verbsty.GT.0) WRITE(*,*) '*WARNING: zero Cd coints'
	        cmeancd = -1.0
	      ENDIF
	      IF(verbsty.GE.2) WRITE(*,*) nadc,cmeancd,cmean-1.0,fwhm
	      IF(verbsty.GE.2) WRITE(*,*) ' '
C
	      efffe = 1.
	      effcd = 1.
C
	      IF(cmeancd.GT.0.0 .AND. cmeanfe.GT.0.0 ) THEN
	        gorow = gorow + 1
	        CALL FTpcld(ounit,1,gorow,1,1,tmid,ftstat) 
	        CALL FTpcle(ounit,2,gorow,1,1,gainfe,ftstat)
	        CALL FTpcle(ounit,3,gorow,1,1,gaincd,ftstat)
	        CALL FTpcle(ounit,4,gorow,1,1,cmeanfe,ftstat)
	        CALL FTpcle(ounit,5,gorow,1,1,cmeancd,ftstat)
	        CALL FTpcle(ounit,6,gorow,1,1,efffe,ftstat)
	        CALL FTpcle(ounit,7,gorow,1,1,effcd,ftstat)
	        CALL FTpclj(ounit,8,gorow,1,1,nbin,ftstat)
	        IF(ftstat.NE.0) THEN
		  CALL fcerr('ERROR: cannot write gain record')
		  RETURN
	        ENDIF
	      ENDIF
C
	    ENDIF
C
   40	    nstep = nstep + 1
	    time2 = time1 + nstep * tstep
	    IF(time2.LE.tev) GOTO 40	      
	    DO n = 1,32
	      IF(verbsty.GE.3) WRITE(*,*) n,esum(n)
	      esum(n) = 0.
	    ENDDO
	    esum(e) = sumfac
	    nbin = 1
C
   50	CONTINUE
C
C= process mext buffer if more events remain
C
	frow = frow + nelem
	remain = remain - nelem
	IF(remain.GT.0) GOTO 10
C
	RETURN
	END
C
C***************************************************************************** 
C SUBROUTINE:
C      xmktie
C
C DESCRIPTION:
C      This subroutine moves the extra keywords,i.e., the
C      keywords which don't contain: XTENSION,BITPIX,
C      NAXIS*,PCOUNT,GCOUNT,TFIELDS,TTYPE*,TFORM*,TUNIT*,
C      END,TBCOL*,and EXTNAME, from the input file to 
C      the output file
C
C AUTHOR/DATE:
C      James Kent Blackburn 6/09/92
C
C MODIFICATION HISTORY:
C
C
C NOTES:
C
C USAGE:
C      call xmktie(iunit,ounit,status)
C
C ARGUMENTS:
C      iunit - input unit number
C      ounit - output unit number
C      fname - input FITS file name
C      status  - error number
C
C PRIMARY LOCAL VARIABLES:
C      i* -index to substring
C      l* - substring presence flag
C      nkeys - number of keywords
C      copyflg - copy keyword flag
C
C CALLED ROUTINES:
C      subroutine ftghsp - get number of keywords in extension
C      subroutine ftgrec - get keyword record
C      subroutine ftprec - put keyword record
C
C***************************************************************************
        subroutine xmktie(iunit,ounit,status)

        integer     iunit,ounit,status
        integer     i1,i2,i3,i4,i5,i6,i7,i8,i9,ia,ib,ic,id
        integer     i,nkeys,nmore
        logical     l1,l2,l3,l4,l5,l6,l7,l8,l9,la,lb,lc,ld
        logical     copyflg
	character*80 record
C
        CALL ftghsp(iunit,nkeys,nmore,status)
        DO 10 i = 1, nkeys
          CALL ftgrec(iunit,i,record,status)
          i1 = index(record(1:8),'XTENSION')
          i2 = index(record(1:6),'BITPIX')
          i3 = index(record(1:5),'NAXIS')
          i4 = index(record(1:6),'PCOUNT')
          i5 = index(record(1:6),'GCOUNT')
          i6 = index(record(1:7),'TFIELDS')
          i7 = index(record(1:5),'TTYPE')
          i8 = index(record(1:5),'TFORM')
          i9 = index(record(1:5),'TUNIT')
          ia = index(record(1:7),'EXTNAME')
          ib = index(record(1:3),'END')
          ic = index(record(1:5),'TBCOL')
	  id = index(record(1:7),'HDUCLAS')
          l1 = i1 .eq. 0
          l2 = i2 .eq. 0
          l3 = i3 .eq. 0
          l4 = i4 .eq. 0
          l5 = i5 .eq. 0
          l6 = i6 .eq. 0
          l7 = i7 .eq. 0
          l8 = i8 .eq. 0
          l9 = i9 .eq. 0
          la = ia .eq. 0
          lb = ib .eq. 0
          lc = ic .eq. 0
	  ld = id .eq. 0
          copyflg = l1 .and. l2 .and. l3 .and. l4 .and. l5
     &        .and. l6 .and. l7 .and. l8 .and. l9 .and. la
     &        .and. lb .and. lc .and. ld
          IF( copyflg ) THEN
            CALL ftprec(ounit,record,status)
          ENDIF
   10   CONTINUE
C
        RETURN
        END

*******************************************************************************