      program fitlinc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                          F I T L I N C
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  scanfile
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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit real*8(a-h,o-z)
      character*72 ifile,ofile
      character*1 ch
      logical lcut,lcol
      real*8 tmin(10000),hmin(10000),d(10000),td(10000)
      real*8 xhr(24000000),rlat(24000000)
      real*4 thr(24000000),h(24000000)
      integer iamp(24000000)
c  collocation dimension: icdim = (maxobs*(maxobs+1))/2+maxobs
      dimension c(5015000)
      dimension nhist(100)
      nmax = 24000000
      maxobs = 10000
      maxhist = 100
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='old')
      open(20,file=ofile,status='unknown')
      open(14,file='fitlinc.ave',status='unknown')
      open(15,file='fitlinc.hst',status='unknown')
c
c  read in all data
c
      write(*,*) '- reading scanner data from: ',ifile
      i = 1
20    if (lcut) then
        read(10,*,end=30) tt,thr(i),xhr(i),h(i),iamp(i)
      else
        read(10,*,end=30) thr(i),rlat(i),xhr(i),h(i),iamp(i)
      endif
      if (h(i).gt.hcutoff) goto 20
      i = i+1
      if (i.gt.nmax) stop 'too many points - increase nmax'
      goto 20
30    n = i-1
      t1 = thr(1)
      write(*,31) n,t1,thr(n)
31    format(' read',i7,' points in interval:',2f9.4)
c
c  now find minimum values 
c
      if (dtmin.eq.0) goto 80
      nint = (thr(n)-t1)/dtmin + 0.5
      if (nint.gt.maxobs) stop 'too many minimum intervals'	
c
      do j = 1, nint
        tmin(j) = -9999
        hmin(j) = 9999
      enddo
c
      do i = 1, n
        j = (thr(i)-t1)/dtmin+1
        if (h(i).lt.hmin(j)) then
          hmin(j) = h(i)
          tmin(j) = thr(i)
        endif
      enddo
c  
      write(*,*) 'mimimum values:'
      write(*,32) (hmin(j),j=1,nint)
32    format(' ',12f6.1)
c
c  now detrend with linear funtion - reject outliers one by one
c  flag rejected points by time 9999
c  ------------------------------------------------------------
c
      iloop = 1
c
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)-t1
        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)-t1
        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+t1,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)
c 
        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
c  
        k = 0
        do 120 i = 1,nsel
        do 120 j = 1,i
          k = k+1
          r = abs(td(i)-td(j))
c
c  second order markov covariance function
c
          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)-t1))
        enddo
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, n
        hc = bias + alfa*(thr(i)-t1)
        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
c
        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,50000).eq.1) 
     .    write(*,50) i,throut,xhrout,hh,iamp(i),hc
          write(20,50) i,throut,xhrout,hh,iamp(i),hc
50        format(i7,2f12.6,f9.2,i4,f7.2)
          rlatsum = rlatsum + throut
          rlonsum = rlonsum + xhrout
        else
          if (mod(i,50000).eq.1) 
     .    write(*,51) throut,rlat(i),xhr(i),hh,iamp(i),hc
          write(20,51) throut,rlat(i),xhr(i),hh,iamp(i),hc
51        format(f9.6,f10.6,f12.6,f8.2,1x,i3,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,f12.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
      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):'
      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
      end
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
c
      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

