        program LLF3glas
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                          L L F 3 G L A S
c  
c  program to detrend sea-ice scanner data by linear fit with
c  outlier rejection and optional collocation filtering of minimal values
c  will also make histogram of freeboard data
c    
c  input:
c
c  ICESat file (dech,lat,lon,height)
c  outfile
c  dtmin, rejl, hcutoff, dtave, histsp, lcut, lcol
c  thalf, rmsn    (if lcol true)
c
c  where
c
c  dtmin: time spacing (dechr) for taking min value (0: just cut) 
c  rejl:  limit for rejecting minimum values from linear fit
c  hcutoff: max scan height allowed (remove clouds)
c  dtave:  along-track averaging time for "fitlinc.ave"
c  histsp: histogram spacing 
c  lcut: false: cut output in 0.02 dechr segments for plotting with pse
c        (assuming readscan thr,cross-trk output)
c        true: standard lat,lon output (lgeo=true in readscan)
c  lcol: use collocation after detrending, with corr.length thalf (hr) and 
c        noise rmsn (c0 scaled from residuals). 2nd order Markov cov function.
c        warning: program makes one solution with all data - so keep number of
c        minimum points small (solves no of linear equations = no of min pts)! 
c
c  outputfile contains freeboard heights and amplitude
c  "fitlinc.hst" contains histogram information
c  "fitlinc.ave" contains averaged free-boards (t,lat,lon,h_fb,n_points)
c
c  (c) Rene Forsberg/National Survey and Cadastre (KMS)
c  thule air base 19 may 2002
c  updated Alert 19 may 2004 (reading new format readscan data + collocation)
c  changed by Henriette Skourup (DNSC) July 2005 to read ICESat data files
c  including more than one orbit. Sort on time - for gaps larger than 1h
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       n number of lines in each orbit
c       no number of orbits
c       nint number of minimum points
c       nh number of intervals in histogram 
	implicit real*8(a-h,o-z)
	character*72 ifile,ofile
	character*1 ch
	logical lcut,lcol
	real*8 tmin(1000),hmin(1000),d(1000),td(1000)
	real*8 tminx(1000,100),hminx(1000,100)
	real*8 tminsum(1000),hminsum(1000)
	real*8 hsum50(1000),have50(1000)
	real*8 xhr(25000),rlat(25000)
	real*4 thr(25000),h(25000),gpst(25000)
	real*4 ino(25000)
	real*8 ha(25000)
	integer nlt,nobs,ks,ksend
	integer jint,nint,maxmin,jmin
	integer iamp(25000)
c	collocation dimension: icdim = (maxobs*(maxobs+1))/2+maxobs
	dimension c(501500)
	dimension nhist(100)
	nmax = 25000
	maxobs = 2500
	maxhist = 100
	maxmin = 3
c
	read(*,'(a)') ifile
	read(*,'(a)') ofile
	read(*,*) dtmin,rejl,hcutoff,dtave,histsp,lcut,lcol
	if (lcol) then
	  read(*,*) thalf,rmsn
	endif
c
	write(*,*) 'FITLINC - linear detrending and collocation filtering'
c
	open(10,file=ifile,status='unknown')
	open(20,file=ofile,status='unknown')
        open(14,file='fitlinc.ave',status='unknown')
        open(15,file='fitlinc.hst',status='unknown')

c---------------------------------------------------------------
c	Read in data 
c	More than one orbit if timegap greater than 1h
c---------------------------------------------------------------
	write(*,*) '- reading data from: ',ifile

	tdiff = 0
	ks = 1
	ksend = 0

	nlt=1
	i = 1	
10	read(10,*,end=34) x1,x2,x3,x4,x5
	if (ks.eq.1) then
	 ino(i) = x1
	 rlat(i) = x2
	 xhr(i) = x3
	 gpst(i) = x4
	 h(i) = x5
	 thr(i) = gpst(i)/3600
	 if (i.eq.1) GPSt(0) = GPSt(1)
	 tdiff = GPSt(i) - GPSt(i-1)
	else
	 ino(i) = x1
	 rlat(i) = x2
	 xhr(i) = x3
	 gpst(i) = x4
	 h(i) = x5
	 thr(i) = gpst(i)/3600
	 tdiff = GPSt(i) - GPSt(i-1)
	end if	 
        if (abs(tdiff).gt.3600) goto 35
	i = i+1
	nlt = nlt+1
	goto 10

34	continue
	ksend = ks

c       ks number of orbits
c	nobs number of observations in each orbit
35      nobs = i-1

c---------------------------------------------------------------
c	Starting to do the calculation for each interval
c---------------------------------------------------------------
36	if (dtmin.eq.0) goto 200
	nint = (thr(nobs)-thr(1))/dtmin + 0.5
	write(*,*) nint,thr(nobs),thr(1),dtmin,gpst(nobs),gpst(1)
	if (nint.gt.maxobs) stop 'too many minimum intervals'	
c
c---------------------------------------------------------------
c	now find minimum values 
c	nint number of minimum values 
c---------------------------------------------------------------
	do j = 1, nint
	do k = 1, maxmin
		tminx(j,k) = -9999
		hminx(j,k) = 9999
	end do
	end do


	do i = 1,nobs
	jint = (thr(i)-thr(1))/dtmin+1
	if (h(i).lt.hminx(jint,1)) then
		hminx(jint,3) = hminx(jint,2)
		tminx(jint,3) = tminx(jint,2)
		hminx(jint,2) = hminx(jint,1)
		tminx(jint,2) = tminx(jint,1)
		hminx(jint,1) = h(i)
		tminx(jint,1) = thr(i)
	else if (h(i).lt.hminx(jint,2)) then
		hminx(jint,3) = hminx(jint,2)
		tminx(jint,3) = tminx(jint,2)
		hminx(jint,2) = h(i)
		tminx(jint,2) = thr(i)
	else if (h(i).lt.hminx(jint,3)) then
		hminx(jint,3) = h(i)
		tminx(jint,3) = thr(i)
	end if
	end do

	do j = 1, nint
	if (hminx(j,3).eq.9999.or.hminx(j,2).eq.9999.
     .	   or.hminx(j,3).eq.9999) then
	do k = 1, maxmin
		tminx(j,k) = -9999
		hminx(j,k) = 9999
	end do
	end if
	end do
c
	do j = 1, nint
	tminsum = 0
	hminsum = 0
	  do k = 1, maxmin
            tminsum(j) = tminsum(j) + tminx(j,k) 
            hminsum(j) = hminsum(j) + hminx(j,k)
	  end do
	  tmin(j) = tminsum(j)/maxmin 
	  hmin(j) = hminsum(j)/maxmin
	end do

	write(*,31) ks,nobs,thr(1),thr(nobs)
31      format(///' orbit',i5,' read',i7,' points in interval:',2f9.4)
	write(*,*) 'Number of intervals: ',nint
	write(*,*) 'Mimimum values:'
	write(*,32) (hminx(j,1),j=1,nint)
	write(*,*) ' '
	write(*,32) (hminx(j,2),j=1,nint)
	write(*,*) ' '
	write(*,32) (hminx(j,3),j=1,nint)
	write(*,*) ' '
	write(*,32) (hmin(j),j=1,nint)
32	format(' ',10f7.1)




c---------------------------------------------------------------
c  now detrend with linear funtion - reject outliers one by one
c  flag rejected points by time 9999
c  ------------------------------------------------------------
        iloop = 1
	bias = 0
	alfa = 0

40      nsum = 0
        sumt = 0
        sumt2 = 0
        sumh = 0
        sumth = 0
        resmax = 0
        ressum2 = 0
        do 41 j = 1, nint
          if (tmin(j).eq.9999) goto 41
            tt = tmin(j)-thr(1)
            hh = hmin(j)
            res = abs(hh-(bias+alfa*tt))
          if (iloop.gt.1) then
            ressum2 = ressum2 + res**2
            if (res.gt.resmax) then
              resmax = res
              jmax = j
            endif
          endif
        nsum = nsum + 1
        sumt = sumt + tt
        sumt2 = sumt2 + tt**2
        sumh = sumh + hh
        sumth = sumth + tt*hh
41      continue
c
        if (resmax.ge.rejl) then
          tt = tmin(jmax)-thr(1)
          hh = hmin(jmax)
          nsum = nsum - 1
          sumt = sumt - tt
          sumt2 = sumt2 - tt**2
          sumh = sumh - hh
          sumth = sumth - tt*hh
          tmin(jmax) = 9999
          write(*,42) jmax,tt+thr(1),hh,bias+alfa*tt
42        format(' rej j,t,hh,hlin = ',i3,f9.4,3f9.2)
        elseif (iloop.gt.1) then
          goto 49
        endif
c
        r = nsum*sumt2 - sumt**2
        if (r.eq.0) r = 1.d-99
        bias = (sumt2*sumh - sumt*sumth)/r
        alfa = (-sumt*sumh + nsum*sumth)/r
c
        iloop = iloop+1
        goto 40
c
49      write(*,*) 'number of rejected min values: ',nint-nsum,
     .  ' out of ',nint
        write(*,491) sqrt(ressum2/nsum)
491     format(' r.m.s. residual of linear fit:',f9.2)
        nsel = nsum

c
c  collocation prediction 
c  ----------------------
c  set up collocation equations
c  first move together arrays 
c
80    if (lcol) then 
        c0 = ressum2/nsum
        beta = 0.595*thalf
        cnn = rmsn**2
        write(*,101) thalf,rmsn,nsel
101     format(/' Collocation min value interpolation, thalf:',f6.3,
     .  ' hr, rmsn:',f6.2/' number of equations:',i4)

        j = 1
        do i = 1, nint
          if (tmin(i).lt.9998) then 
	    d(j) = hmin(i)
            td(j) = tmin(i)
            j = j+1
          endif
        enddo

        k = 0
        do 120 i = 1,nsel
        do 120 j = 1,i
	k = k+1
          r = abs(td(i)-td(j))

c  second order markov covariance function
          cov = c0*(1 + r/beta)*exp(-r/beta)
          if (i.eq.j) cov = cov + cnn
          c(k) = cov
120     continue
        irhs = k
        do i = 1,nsel
          c(irhs+i) = d(i) - (bias + alfa*(td(i)-thr(1)))
        end do
c
c  solve equations
c
        call chol(c,nsel,nsing)
        if (nsing.ne.0) write(*,*)
     .  '*** warning: singular equations - nsing = ',nsing
        write(*,*) '- Collocation equations solved -'
      endif
c
c  output data, cut in slices if lcut=true, geo output if lcut=false
c  -----------------------------------------------------------------
c
      hhmin = 9999
      hhmax = -9999
      hhsum = 0
      colsum2 = 0
      if (histsp.le.0) stop '*** histsp must be positive'
      nh = 3.0/histsp+1
      if (nh.gt.maxhist) nh = maxhist
      do i = 1,nh
        nhist(i) = 0
      enddo
c
      tave0 = thr(1)
      tsum = 0
      nsum = 0
      rlatsum = 0
      rlonsum = 0
      havesum = 0 
c  
      do i = 1, nobs
        hc = bias + alfa*(thr(i)-thr(1))
        cbias = 0
        if (lcut) then
          cbias = int((thr(i)-thr(1))*50)/50.d0
          throut = thr(i) - cbias
          xhrout = xhr(i) + cbias*0.1
        else
          throut = thr(i)
        endif
c
        hh = h(i)-hc
c
c  collocation prediction signal
        if (lcol) then
          sum = 0
          do k = 1, nsel
            r = abs(td(k)-throut) 
            cov = c0*(1 + r/beta)*exp(-r/beta)
            sum = sum + cov*c(irhs+k)
          enddo
          hh = hh - sum
          hc = hc + sum
          colsum2 = colsum2 + sum**2
        endif
c
c  histogram
c
        if (hh.lt.hhmin) hhmin = hh
        if (hh.gt.hhmax) hhmax = hh
        hhsum = hhsum + hh
        ii = hh/histsp+2
        if (ii.ge.1.and.ii.le.nh) nhist(ii) = nhist(ii)+1
c
c  output 
c
        if (lcut) then
          if (mod(i,500000).eq.1) 
     .    write(*,50) i,throut,xhrout,hh
c     ,iamp(i),hc
          write(20,50) i,throut,xhrout,hh
c	  ,iamp(i),hc
c50        format(i7,2f11.6,f9.2,i4,f7.2)
50        format(i7,2f11.6,f9.2)
          rlatsum = rlatsum + throut
          rlonsum = rlonsum + xhrout
        else
          if (mod(i,500000).eq.1) 
     .    write(*,51) throut,rlat(i),xhr(i),hh,hc,h(i)
          write(20,51) throut,rlat(i),xhr(i),hh,hc,h(i)
51        format(f10.6,f11.5,f12.5,2f7.2,f8.2)

c          if (mod(i,500000).eq.1) 
c     .    write(*,51) throut,rlat(i),xhr(i),hh,hc,h(i),ha(i)
c          write(20,51) throut,rlat(i),xhr(i),hh,hc,h(i),ha(i)
c	51        format(f10.6,f11.6,f12.6,2f7.2,f8.2,f7.2)
          rlatsum = rlatsum + rlat(i)
          rlonsum = rlonsum + xhr(i)
        endif
        nsum = nsum + 1
        tsum = tsum + throut
        havesum = havesum + hh
c
c  average output
c
        if ((thr(i)-tave0).gt.dtave) then
          if (nsum.gt.0) write(14,52) 
     .    tsum/nsum,rlatsum/nsum,rlonsum/nsum,havesum/nsum,nsum
52        format(f9.6,f10.6,f11.6,f8.2,1x,i5)
          tave0 = thr(i)
          tsum = 0
          nsum = 0
          rlatsum = 0
          rlonsum = 0
          havesum = 0 
        endif
c
      enddo
c
c  statistics and histogram
c  ------------------------
c
      n = nobs
      if (lcol) write(*,54) sqrt(colsum2/n)
54    format(' r.m.s. collocation signal: ',f9.2)
      write(*,55) hhmin,hhmax,hhsum/n
55    format(' min, max and average output values: ',3f9.2)
c
      write(*,*) 
      write(*,*) 'Histogram of freeboards (percent):'
      write(15,58) ks
58    format('Histogram orbit',i5)
      do i = 1,nh
        ch = '*'
        pct = nhist(i)*100.d0/n
        ipct = pct/2.5 + 0.5
        write(*,56) histsp*(i-2),histsp*(i-1),nhist(i),pct,(ch,j=1,ipct)
56      format(f6.2,' -',f5.2,i7,f9.3,' !',41a)
        write(15,57) histsp*(i-2),histsp*(i-1),nhist(i),pct
57      format(f6.2,f6.2,i6,f9.3)
      enddo




c---------------------------------------------------------------
	ino(1) = ino(nobs+1)
	rlat(1) = rlat(nobs+1)
	xhr(1) = xhr(nobs+1)
	gpst(1) = gpst(nobs+1)
	h(1) = h(nobs+1)
	thr(1) = thr(nobs+1)
	if (ks.eq.ksend) goto 200
	ks = ks + 1
	i = 2
	goto 10

200	continue
	write(*,*) ' '
	write(*,*) 'End of File'
        write(*,*) 'Total number of lines in data file: ',nlt-1
	write(*,*) 'Number of orbits in file: ',ks-1
	close(10)
	close(20)
	close(14)
	close(15)
	end
c
c
c
c
c
c
c
	subroutine chol(c,n,nsing)
	implicit double precision(a-h,o-z)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                     c h o l
c
c  subroutine solves positive definite symmetric linear equations
c  using cholesky decomposition. coefficients stored columnwise
c  in c, followed by righthand side. n is number of unknowns.
c  solution is returned as last column in c.
c  'nsing' is number of singularities. it should be zero for a
c  succesful solution.
c
c  rf, nov 85
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      dimension c(501500)
      nsing = 0

      do 50 nr = 1,n+1
        i=nr*(nr-1)/2
        ir=i
        do 40 nc = 1,nr
          if (nc.gt.n) goto 50
          sum=0
          ic=nc*(nc-1)/2
          i=i+1
          nc1=nc-1
          do 30 np = 1,nc1
30        sum = sum - c(ir+np)*c(ic+np)
          ci = c(i)+sum
          if (nr.eq.nc) goto 31
          cc = c(ic+nc)
          if (cc.eq.0) then
            write(*,*) '*** Cholesky reduction zero-division'
            cc = 1.0e9
          endif
          c(i) = ci/cc
          goto 40
31        if (nr.gt.n) goto 40
          if (ci.gt.0) goto 32
            nsing = nsing+1
            c(i) = 1.0e9
            goto 40
32        c(i) = sqrt(ci)
40      continue
50    continue
c
c  back substitution
c
      do 80 nc = n,1,-1
        ir=i
        ic=nc*(nc+1)/2
        cc = c(ic)
        if (cc.eq.0) then
          write(*,*) '*** Cholesky subsitution zero-division'
          cc = 1.0e9
        endif
        c(i) = c(i)/cc
        do 70 np = nc-1,1,-1
          ir=ir-1
          ic=ic-1
          c(ir) = c(ir) - c(i)*c(ic)
70      continue
        i = i-1
80    continue
      return
      end