;The procedure speccorrec_ift_wc.pro calculates correction factors for fluxes
;based on signals with spectral degradation for Willow Creek. 
;The h2o and co2 signals from the trailer licors have spectral degradation 
;due to the long delivery tubes.
;
;Methodology is to degrade the temp spectrum based on the shape of the 
;flawed h2o spectrum or based on the theory proposed by Lenschow and 
;Raupach (JGR 1991) for co2. The corrected flux is then given by multiplying
;the flux estimate from the flawed h2o or co2 signals by the ratio of the
;heat flux computed with the normal temp signal and the heat flux computed
;with the degraded temp signal. The degraded temp signal is found by
;computing the inverse Fourier transform of the degraded spectrum. The 
;degraded spectrum is found by multiplying the actual spectrum by a
;transfer function (called phi here). This program is hopefully
;commented enough so that the exact method is clear.
;
;The difference in method between h2o and co2 is that where h2o seems to
;have problems at large scales due to 'sticking' of water on the tube walls,
;the co2 does not and perhaps has only minimal diffusion-related spectral
;degradation at high frequency. The degradation in the c spectrum is
;not visible due to noise at the high frequency end of the spectrum. 
;The h2o spectra, however, are quite clearly degraded well before the 
;noise occurs and provide a way to determine the correction factor as 
;described above and below.
;
;Use is as follows:
;       speccorrec_ift_wc,silent,dq,dc,dtq,dtc,ndt,dwq,dwc,lagc,$
;                         q_cfactor,qflag,$
;                         c_cfactor,cflag
;Where
;        |        silent = 0 to plot spectra used to correct q flux
;        |                 1 otherwise
;        |         dq,dc = detrended and lag adjusted q and c, respectively
; input  |       dwq,dwc = detrended and lag adjusted vertical velocity
;        |                 corresponding to q and c, respectively
;        |       dtq,dtc = detrended and lag adjusted sonic temp
;        |                 corresponding to q and c, respectively.
;        |                 Note: lengths and time should match for
;        |                       dq,dwq,dtq and for dc,dwc,dtc
;        |           ndt = non-detrended sonic temp
;        |          lagc = lag time in seconds for co2
;
;        |     q_cfactor = h2o flux correction factor
;        |                  wq_corr=wq*q_cfactor
;output  |         qflag = quality flag for q_cfactor:
;        |                    = 0 if q_cfactor is good (between 1.4 and 1.0)
;        |                    = bad q_cfactor if q_cfactor is out of range 
;        |                       (then q_cfactor=1)
;        |                    = -10 if data has too many badvals
;        |                    = -20 if data has length problems
;        |                    = -30 if data has too many badvals & length probs
;        |     c_cfactor = co2 flux correction factor
;        |                  wc_corr=wc*c_cfactor
;        |         cflag = quality flag for c_cfactor:
;        |                    = 0 if c_cfactor is good (between 1.1 and 1.0)
;        |                    = bad c_cfactor if c_cfactor is out of range 
;        |                       (then c_cfactor=1)
;        |                    = -10 if data has too many badvals
;        |                    = -20 if data has length problems
;        |                    = -30 if data has too many badvals & length probs
;
;Written by BWBerger 4/20/99 (based heavily on speccorrec_ift.pro for wlef).
;Modified by BWB 4/19/00 to correct phi. See notes written below in
; text for full description.

pro speccorrec_ift_lc,silent,dq,dc,dtq,dtc,ndt,dwq,dwc,lagc,$
                   q_cfactor,qflag,$
                   c_cfactor,cflag

;Frequency of data in Hz.
fs=10

;Set a flag used to be sure the noise 'knee' is found
badkneeflag=0
;Set a flag used to be sure the t spectrum is is greater than the q spectrum
badtspecflag=0

;If operator wants to watch this routine run print a separator line
if silent eq 0 then begin
   subtitle='Licor Spectral Correction'
   print,'_______________________'
endif
printf,8,'_______________________'

;See if there's enough good data to enable finding spectral correction
;get percentage of good values in data and be sure lengths match
lenq=len(dq)
lenc=len(dc)
lentq=len(dtq)
lentc=len(dtc)
lenwq=len(dwq)
lenwc=len(dwc)
qlen_ok=0
clen_ok=0

if (lenq eq lentq) and (lenq eq lenwq) then begin
   qlen_ok=1
endif else begin
   print,'h2o data length problems in speccorrec_ift'
   printf,8,'h2o data length problems in speccorrec_ift'
endelse
if (lenc eq lentc) and (lenc eq lenwc) then begin
   clen_ok=1
endif else begin
   print,'co2 data length problems in speccorrec_ift'
   printf,8,'co2 data length problems in speccorrec_ift'
endelse
whereq=where(dq eq -999.,countq)
wherec=where(dc eq -999.,countc)
wheretq=where(dtq eq -999.,counttq)
wheretc=where(dtc eq -999.,counttc)
wherewq=where(dwq eq -999.,countwq)
wherewc=where(dwc eq -999.,countwc)
percentbadq=100.0*countq/lenq
percentbadc=100.0*countc/lenc
percentbadtq=100.0*counttq/lentq
percentbadtc=100.0*counttc/lentc
percentbadwq=100.0*countwq/lenwq
percentbadwc=100.0*countwc/lenwc
q_gooddata=0
c_gooddata=0
if (percentbadq lt 20.) and (percentbadtq lt 20.) and (percentbadwq lt 20.) then begin
   q_gooddata=1
endif else begin
   print,'Too many badvals in q data for speccorrec_ift'
   printf,8,'Too many badvals in q data for speccorrec_ift'
endelse
if (percentbadc lt 20.) and (percentbadtc lt 20.) and (percentbadwc lt 20.) then begin
   c_gooddata=1
endif else begin
   print,'Too many badvals in c data for speccorrec_ift'
   printf,8,'Too many badvals in c data for speccorrec_ift'
endelse
;********************************************************

;&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;Get temporary variables that are originals truncated to largest
;factor of 2. This assures quickest fft operation.
;2^14=16384 (2^15=32768 too big, only have 18000 max without lag)
;*****************
dc_temp=dc(0:16384-1)
dq_temp=dq(0:16384-1)
dtc_temp=dtc(0:16384-1)
dtq_temp=dtq(0:16384-1)
dwc_temp=dwc(0:16384-1)
dwq_temp=dwq(0:16384-1)

;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;Find correction factor for h2o (q)
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (qlen_ok eq 1) and (q_gooddata eq 1) then begin
   ;Fill in badvals in data (currently, does this for any gap size and
   ;this may cause problems spectrally if the gaps are large)
   zapbadval,1,transpose(dq_temp),dqc,tableq
   zapbadval,1,transpose(dtq_temp),dtqc,tabletq
   zapbadval,1,transpose(dwq_temp),dwqc,tablewq

   ;Get length of arrays
   nq=len(dqc)

   ;Determine if the length is even or odd. If odd, set to next lowest even.
   ;This is so frequency vector is correct for fft output.
   evq=2*floor(nq/2)
   if nq gt evq then begin
      nq=nq-1
      dqc=dqc(0:nq-1)
      dtqc=dtqc(0:nq-1)
      dwqc=dwqc(0:nq-1)
   endif
   ;Remove any mean that may still exist - variables are detrended coming 
   ;into speccorrec_ift so this mean is probably close to zero
   qmn=mean(dqc)
   tqmn=mean(dtqc)
   wqmn=mean(dwqc)
   q=dqc-qmn
   tq=dtqc-tqmn
   wq=dwqc-wqmn

   ;Do fft and calculate spectrum for q and tq
   aq=fft(q,/double)
   aqhat=double(conj(aq)*aq)
   delfq=(nq-1)/double(fs)
   gqq=2*delfq*aqhat
   ;*******
   atq=fft(tq,/double)
   athat=double(conj(atq)*atq)
   gtt=2*delfq*athat

   ; Want only half of spectrum (remove frequency values larger than Nyquist).
   halfq=nq/2+1
   gqq=gqq(0:halfq-2)
   gtt=gtt(0:halfq-2)
   fq=dindgen(halfq-1)*fs/(nq-1)

   ;**************************************************

   ;Smooth the spectra
   gttsm=marathon(gtt(1:len(gtt)-1),11)
   gqqsm=marathon(gqq(1:len(gqq)-1),11)
   fsm=fq(1:len(fq)-1)

   ;Trim off zeros created by marathon (11 pts are boxcar averaged and
   ;the result is placed at the center of the boxcar, so the 
   ;smoothed series starts at the sixth point of the original series and
   ;ends six points from the end of the original series - the output of
   ;marathon is the same length as the original series because the first
   ;and last five points are set to zero. These are the points trimmed
   ;off here.)
   gttsmc=gttsm(5:len(gttsm)-6)
   gqqsmc=gqqsm(5:len(gqqsm)-6)
   fsmc=fsm(5:len(fsm)-6)

   ;Use the first 15 points to determine the normalizing variance
   ;This is an arbitrary choice in an attempt to use scales where the
   ;spectra should have a similar shape.
   tt=(1/delfq)*total(gttsmc(0:14))
   qq=(1/delfq)*total(gqqsmc(0:14))

   ;Determine what frequency range of q spectrum to fit
   ;Trailer licors have a 'knee' above which is noise that is to be excluded.
   ;In situ licors sometimes have a little high frequency noise, but this
   ;is not excluded at this time (BWB 10/3/99)
   ;Find index in frequency vector for smoothed spectra (fsmc) at knee.
   ;The knee is where noise bends the q spectra up.
   ;This is done using a square wavelet
   edg=edgedetect(gqqsmc,201)
   wav=edg(0,*)
   tm=edg(1,*)

   ;Where the wavelet becomes nearly zero (and constant) is the location
   ;of the knee plus half the width of the wavelet - ie it doesn't become 
   ;constant until the knee sits at the left edge of the wavelet as the
   ;wavelet moves from left (low freq) to right (high freq).
   ;findi finds the index in wav where wav=.001
   tmind=findi(wav,.001,1)
   ;Code added 6/28/99 to deal with tmind=-1 (i.e., could not find knee)
   if tmind eq -1 then begin
;      print,'Knee cannot be found. No H2O spectral corrections performed. lev=',strtrim(string(lev),2)
;      printf,8,'Knee cannot be found. No H2O spectral corrections
;      performed. lev=',strtrim(string(lev),2)
      print,'Knee cannot be found. No H2O spectral corrections performed.'
      printf,8,'Knee cannot be found. No H2O spectral corrections performed.'
      badkneeflag=1 ;knee couldn't be found
      hiind=len(fsmc)      ;set hiind to end of fsmc and get value of
      hifreq=fsmc(hiind-1) ;fsmc there for plotting purposes if silent=0
   endif else begin ;knee found (hopefully)
      ;Width of wavelet is 201 points so knee is 100 pts to left of where
      ;wav becomes zero and constant.
      hiind=tm(tmind)-100
      if hiind le 2 then begin ;make sure knee isn't right at start of tm
         print,'Difficulty finding knee. No H2O spectral corrections performed'
         printf,8,'Difficulty finding knee. No H2O spectral corrections performed'
         badkneeflag=1 ;knee couldn't be found
         hiind=len(fsmc)      ;again, as for tmind=-1 case, get hifreq for
         hifreq=fsmc(hiind-1) ;plotting purposes
      endif else begin ;knee appears to be good
         ;Finally, we can find the frequency where the knee is.
         hifreq=fsmc(hiind)
         if silent ne 1 then begin
            print,'f at knee=', hifreq
         endif
         printf,8,'f at knee=', hifreq
      endelse
   endelse

   ;**************
   ;Now we're ready to do fits of the t and q spectra
   ;This is done in log coordinates

   ;Get fit coeffiecients by fitting full length of t
;   att=poly_fit(alog10(fsmc),alog10(gttsmc/tt),2)
   att=poly_fit(alog10(fsmc),alog10(gttsmc/tt),3)

   ;Get fit coeffiecients by fitting from lowest f up to knee of q
   ;If for in situ, fit is done up to highest frequency (no knee).
   aqq1flag=0
   aqq=poly_fit(alog10(fsmc(0:hiind-1)),alog10(gqqsmc(0:hiind-1)/qq),2)
   ;The hope is that using the coefficients of the q spectral fit up to
   ;the knee that the fit would predict a decent spectral shape past the
   ;knee. However, it was found that sometimes the curvature of the q 
   ;spectrum resulted in a q fit past the knee that bent upward
   ;inappropriately. This occured when aqq(2) got bigger than 0.5. In
   ;this case, the fit is reduced to first order.
   if aqq(2) gt 0.5 then begin
      aqq1flag=1
      aqq=poly_fit(alog10(fsmc(0:hiind-1)),alog10(gqqsmc(0:hiind-1)/qq),1)
   endif

   ;Create frequency vector from start of fit to highest freq
   fend=fq(6:len(fq)-1)
   ;Obtain log of freq to be consistent with fitting method
   fendl=alog10(fend)
   ;Now get fitted t and q spectra in non-log coordinates (as they are
   ;when they are output by the fft routine).
;   gttfit=10^(fendl^2*att(2)+fendl*att(1)+att(0))
   gttfit=10^(fendl^3*att(3)+fendl^2*att(2)+fendl*att(1)+att(0))
   if aqq1flag eq 0 then begin
      gqqfit=10^(fendl^2*aqq(2)+fendl*aqq(1)+aqq(0));
   endif else begin
      gqqfit=10^(fendl*aqq(1)+aqq(0));
   endelse

   ;While moving to lower frequency starting at the knee for trailer q or the
   ;highest frequency for in situ q, find the first index
   ;where gqqfit is larger than gttfit. This should be at the frequency where 
   ;the transfer function starts differing from a value of 1. This may not
   ;always happen, but it has been found that when the fluxes are large
   ;the spectral behavior is usually good and this method works well.
   ;Otherwise, the transfer function may be close to 1, and if not, the
   ;flux is so low that the correction contributes little to total daily flux.
   phiflag=0
   if hiind ge 6 then begin ;findi.pro can be a little finicky getting tmind 
                           ;and thus hiind; this assures hiind is big enough
      for i=hiind-6,0,-1 do begin 
         if (gqqfit(i) ge gttfit(i)) and (phiflag eq 0) then begin
            phistart=i
            phiflag=1
         endif

         if i le 0 and  phiflag eq 0 then begin   ;was i=1        8/27/99
            phistart=0                            ;was phistart=1    "
         endif
      endfor
   endif else begin
      phistart=0
   endelse

   lofreq=fq(6+phistart)                       ;was 5+phistart    "
   if silent ne 1 then begin
      print,'f at phistart=',lofreq
   endif
   printf,8,'f at phistart=',lofreq

   ;Create transfer function phi
   ;Note: It may not be clear why the fits are done in log coordinates
   ;and then converted to non-log coordinates again. This is done because 
   ;the shapes of the spectra are clearer (more linear?) in log
   ;coordinates. Once the fits are returned to non-log coordinates the
   ;ratio of the spectra allows the t spectra to be degraded to what the
   ;q spectra would hopefully look like without the noise. 
   ;Basically: gqqfit/gttfit=gtt_degraded/gtt for each frequency.
   ;(An analogous transfer function phi can be defined in log coordinates 
   ;that degrades the t spectra based on the difference between the
   ;spectral fits - ie I think this is something akin to 
   ;log(A/B)=log(A)-log(B)).
   ;Anyway, phi is assumed to be equal to 1 up to the point where the q
   ;spectra start to have bad shape (phistart) and is equal to the ratio
   ;mentioned above from phistart to the highest frequency.
   ;phi1 goes from f=0 up to the Nyquist frequency.
   qphi1=[[fltarr(1,6+phistart)+1],[transpose(gqqfit(phistart:len(gqqfit)-1)/gttfit(phistart:len(gttfit)-1))]]

   ;To do the inverse fft, need the mirror of phi1 to get full length phi.
   qphi2=fltarr(1,len(qphi1))
   for i=0,len(qphi1)-1 do begin
      qphi2(i)=qphi1(len(qphi1)-i-1)
   endfor

   qphi=[[qphi1],[qphi2]]

   ;note: qphi is really phi squared. It's just not called that here.
   ;      Taking the sqrt of qphi below is correct. Things have been
   ;      changed, however, for co2 (see that section below).  BWB 4/19/00
   ;*******************

   ;Get degraded t time series by ifft. The q in the name tdgq denotes that
   ;this is the t time series degraded to help fix q.
   tdgq=double(fft(atq*sqrt(qphi),/inverse,/double))
   ;total degraded t time series (not used at this time)
   tdgqt=tdgq+tqmn

   ;Get statistics for ploting
   baremusic,tdgq,5,1,1,gttdg,fttdg,ttdg
   ;ttdg=(gttdg(0)/(delfq*2))+(1/delfq)*total(gttdg(1:5+phistart))
   gttdgsm=marathon(gttdg(1:len(gttdg)-1),11)
   gttdgsmc=gttdgsm(5:len(gttdgsm)-6)
   ttdg=(1/delfq)*total(gttdgsmc(0:14))

   if silent eq 0 then begin
      subtitle='Spectral Correction'
      xmin=min(fsmc)
      xmax=max(fend)
      ymin=min([gttsmc/tt,gqqsmc/qq,gttfit,gqqfit])
      ymax=max([gttsmc/tt,gqqsmc/qq,gttfit,gqqfit])
      plot,[xmin,xmax],[ymin,ymax],/nodata,/xlog,/ylog,title=subtitle
      oplot,fsmc,gttsmc/tt
      oplot,fsmc,gttdgsmc/ttdg,linestyle=1
      oplot,fsmc,gqqsmc/qq,color=7600000
      oplot,fend,gttfit,color=7999999
      oplot,fend,gqqfit,color=7999999,linestyle=2
      oplot,[hifreq,hifreq],[1e10,1e-12],linestyle=0
      oplot,[lofreq,lofreq],[1e10,1e-12],linestyle=1
      if badtspecflag eq 1 then begin
         xyouts,0.2,0.15,'T spec < Q spec; no correction done.',/normal
      endif
      if badkneeflag eq 1 then begin
         xyouts,0.2,0.1,'Knee difficulty; no correction done.',/normal
      endif
   endif

   ;Obtain correction factor which is the ratio of the heat flux from the 
   ;non-degraded t to the heat flux using the degraded t.
   wtq=total(wq*tq)/nq
   wtdgq=total(wq*tdgq)/nq
   q_cfactor=wtq/wtdgq
   qflag=0
   ;If correction factor is outside of acceptable range set correction
   ;factor to 1 and flag = bad correction factor, otherwise flag = 0.
   if q_cfactor gt 1.40 or q_cfactor lt 1.0 then begin
      qflag=q_cfactor
      q_cfactor=1
   endif
endif else begin
   if (q_gooddata eq 0) and (qlen_ok eq 1) then begin
      q_cfactor=1.
      qflag=-10
   endif else begin
      if (q_gooddata eq 1) and (qlen_ok eq 0) then begin
         q_cfactor=1.
         qflag=-20
      endif else begin
         if (q_gooddata eq 0) and (qlen_ok eq 0) then begin
            q_cfactor=1.
            qflag=-30
         endif
      endelse
   endelse
endelse
;If T spectrum is less than Q spectrum, speccorrec is not possible.
if badtspecflag eq 1 then begin
   q_cfactor=1.
   qflag=-60
endif
;If knee can't be found (for trailer licors only), speccorrec cannot be done
if badkneeflag eq 1 then begin
   q_cfactor=1.
   qflag=-50
endif

;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;Now find correction factor for co2 (c)
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (clen_ok eq 1) and (c_gooddata eq 1) then begin
   ;Fill in badvals in data (currently, does this for any gap size and
   ;this may cause problems spectrally if the gaps are large)
   zapbadval,1,transpose(dc_temp),dcc,tablec
   zapbadval,1,transpose(dtc_temp),dtcc,tabletc
   zapbadval,1,transpose(dwc_temp),dwcc,tablewc
   ;Get length of arrays
   nc=len(dcc)

   ;Determine if the length is even or odd. If odd, set to next lowest even.
   ;This is so frequency vector is correct for fft output.
   evc=2*floor(nc/2)
   if nc gt evc then begin
      nc=nc-1
      dcc=dcc(0:nc-1)
      dtcc=dtcc(0:nc-1)
      dwcc=dwcc(0:nc-1)
   endif
   ;Remove any mean that may still exist - variables are detrended coming 
   ;into speccorrec_ift so this mean is probably close to zero
   cmn=mean(dcc)
   tcmn=mean(dtcc)
   wcmn=mean(dwcc)
   c=dcc-cmn
   tc=dtcc-tcmn
   wc=dwcc-wcmn

   ;Do fft and calculate spectrum for tc
   ;Note: Only the Fourier components and frequency are needed here.
   atc=fft(tc,/double)

   ;Get frequency vector up to Nyquist
   halfc=nc/2+1
   fc=dindgen(halfc-1)*fs/(nc-1)

   ;From Lenshow and Raupach phi(f)=exp(-79.7*a*f^2*L*Re^(-1/8)*U^(-2))
   ;where a=tube radius, f=frequency, L=length of tubing,
   ;U=mean flow speed, and Re=2*a*U/nu (Reynolds number) with
   ;nu=molecular viscosity. More simply, phi=exp(el*f^2) where
   ;el=-79.7*a*L*Re^(-1/8)*U^(-2).
   a=0.000794 ;tube inner radius in meters (=1/32") - according to Bruce Cook
   a=0.0019812 ; tube inner radius in meters--according to Bruce and Ankur (email), 6/18/2003

   L=2.5 ;tube length in meters - according to Bruce Cook for Wcreek
   L=2.7 ; maesured by Cook, 2003, July
   ;U=L/lag with lag in seconds, L in meters
   u=L/(lagc+0.0001); to prevent lagc=0, add 0.0001
   print,'u',u
   printf,8,'u',u
   ;nu is determined with a second order fit based on temperature using 
   ;data I found in the text "Fundamentals of Heat and Mass Transfer", 
   ;Incropera and De Witt, 1990, Wiley.

   temp=meanbadval(ndt)
   if silent ne 1 then begin
      print,'mean temp',temp
   endif
   printf,8,'mean temp',temp
   nu=(1.18e-4*temp^2+8.844e-2*temp+13.427)*1e-6
   if silent ne 1 then begin
      print,'nu',nu
   endif
   printf,8,'nu',nu

   ;Reynolds number
   re=2*a*u/nu
   if silent ne 1 then begin
      print,'re',re
   endif
   printf,8,'re',re

   ;Lenschow and Raupach formulation (note: el is -K in methodology paper)
   if re gt 2300. then begin
      if silent ne 1 then begin
         print,'Turbulent tube flow'
      endif
      printf,8,'Turbulent tube flow'
      el=-79.7*a*L*re^(-1/8)*u^(-2) ;turbulent case
   endif else begin
      if silent ne 1 then begin
         print,'Laminar tube flow'
      endif
      printf,8,'Laminar tube flow'
      el=-0.41*a*L*re*nu/(0.16e-4*u^2) ;laminar case el=exp(-.41*Re*Sc*a*L/U^2)
   endelse
   cphi1=exp(el*transpose(fc)^2)

   if silent ne 1 then begin
      print,'el',el
   endif
   printf,8,'el',el

   ;To do the inverse fft, need the mirror of phi1 to get full length phi.
   cphi2=fltarr(1,len(cphi1))
   for i=0,len(cphi1)-1 do begin
      cphi2(i)=cphi1(len(cphi1)-i-1)
   endfor

   cphi=[[cphi1],[cphi2]]
   ;note: This is the correct phi from Lenschow and Raupach. However,
   ;      it originally was thought by myself to be the spectral ratio when
   ;      in fact it is the concentration ratio. phi squared is the correct
   ;      spectral ratio. Therefore, when the degraded temp is calculated
   ;      below, cphi multiplies the fft of the temp directly, not the
   ;      square root of phi as was written before. BWB 4/19/00

   ;Get degraded t time series by ifft.
   tdgc=double(fft(atc*cphi,/inverse,/double)) ;note: used to be sqrt(cphi)
                                               ;that was incorrect

   ;total degraded t (not used for anything at this time)
   tdgct=tdgc+tcmn

   ;Obtain correction factor which is the ratio of the heat flux from the 
   ;non-degraded t to the heat flux using the degraded t.
   wtc=total(wc*tc)/nc
   wtdgc=total(wc*tdgc)/nc
   c_cfactor=wtc/wtdgc
   cflag=0
   ;If correction factor is outside of acceptable range set correction
   ;factor to 1 and flag = bad correction factor, otherwise flag = 0.
   if c_cfactor gt 1.3 or c_cfactor lt 1.0 then begin
      cflag=c_cfactor
      c_cfactor=1
   endif

endif else begin
   if (c_gooddata eq 0) and (clen_ok eq 1) then begin
      c_cfactor=1.
      cflag=-10
   endif else begin
      if (c_gooddata eq 1) and (clen_ok eq 0) then begin
         c_cfactor=1.
         cflag=-20
      endif else begin
         if (c_gooddata eq 0) and (clen_ok eq 0) then begin
            c_cfactor=1.
            cflag=-30
         endif
      endelse
   endelse

endelse

if silent ne 1 then begin
   print,'c_cfactor=',strtrim(string(c_cfactor),2),', cflag=',strtrim(string(cflag),2)
   print,'q_cfactor=',strtrim(string(q_cfactor),2),', qflag=',strtrim(string(qflag),2)
   print,'_______________________'
endif
printf,8,'c_cfactor=',strtrim(string(c_cfactor),2),', cflag=',strtrim(string(cflag),2)
printf,8,'q_cfactor=',strtrim(string(q_cfactor),2),', qflag=',strtrim(string(qflag),2)
printf,8,'_______________________'

return
end