      program oa_das_v1_9

      include 'parameter.const'

      integer stderr, constants
      parameter (stderr=0,constants=15)
      integer TOPEX,ERS1
      integer alt_flag
      integer   idevNull
      parameter (TOPEX=1,TOPEX_ERR=.05,ERS1=15,ERS1_ERR=.15)

c
      real      x(mxobs), y(mxobs), t_since_92(mxobs), h(mxobs),
     >          uopt(mxlim), topt(mxlim), xopt(2,mxlim),
     >          xi(mxint), yi(mxint), ui(mxint), erru(mxint),
     >          e(mxobs),eopt(mxlim)
      real      devNull(4)
      real      km2lon, km2lat, lon_min, lat_min, lon_max, lat_max
      character buff*80
c
      common    /cblock/ c(mxlim)
      external  corr_model
      common    /corr/   imodel, c1, c2, c3, c4, km2lon, km2lat
      
      pi=4*atan(1.0)

      if (iargc() .lt. 14)
     / stop 'USAGE:oa time parameters <BlahBlahIn >BlahBlahOut'
      call getarg(1, buff)
      if ( index(buff,'.') .EQ. 0 ) then
         read(buff,'(i)') i_cen
         t_cen=float(i_cen)
      else
         read(buff,'(f)') t_cen
      endif   
      
      call getarg(2,buff)
      read(buff,'(f)') lat_min
      call getarg(3,buff)
      read(buff,'(f)') lon_min
      call getarg(4,buff)
      read(buff,'(f)') lat_max
      call getarg(5,buff)
      read(buff,'(f)') lon_max
      call getarg(6,buff)
      read(buff,'(f)') del_lat
      call getarg(7,buff)
      read(buff,'(f)') del_lon
      call getarg(8,buff)
      read(buff,'(i)') limit
      call getarg(9,buff)
      read(buff,'(f)') dist_lim
      call getarg(10,buff)
      read(buff,'(f)') time_lim
      call getarg(11,buff)
      read(buff,'(i)') imodel
      call getarg(12,buff)
      read(buff,'(f)') c1
      call getarg(13,buff)
      read(buff,'(f)') c2
      call getarg(14,buff)
      read(buff,'(f)') c3
      if ((imodel.eq.4) .or. (imodel.eq.6)) then
         call getarg(15,buff)
         read(buff,'(f)') c4
      endif

c these are flags in to rescale back to degrees and that the input is
c deviations rather than a mean field.
      irescale=1
      ih=1

      km2lat=110                                     ! 110km/deg of lat
      km2lon=km2lat*cos(.5*(lat_min+lat_max)*pi/180) ! km2lon=km2lat*cos(lat)

c read observations from stdin  data in free-format: t, y, x, h
c will spit to standard out free-formatted output: xi,yi,int_field,int_err
c note only first mxobs will be read.  If there are more than mxobs they
c are read, but not loaded into arrays.

      do 25 i = 1, mxobs
      read (*,*,end=30) t_since_92(i), y(i), x(i), h(i), alt_flag
      if (alt_flag.eq.TOPEX) then
         e(i)=TOPEX_ERR
      else
         e(i)=ERS1_ERR
      endif
      
   25 continue

c   if this is reached more than mxobs exist.  The excess is flush is flushed  
c   from stdin here
   29   read(*,*,end=30) (devNull(k),k=1,4), idevNull
      if (1) goto 29

   30 continue
       most=i-1
c
c  convert locations in degrees into distances in km from lon_min,lat_min
      call scale_space (1,x,y,most,lon_min,lat_min,km2lon,km2lat)


c
c  remove average from the h field
      call remav (h,most,aveh,varh)
c
c set up the interpolation positions
      nx = 1 + abs(lon_max - lon_min)/del_lon
      ny = 1 + abs(lat_max - lat_min)/del_lat
      write(*,*) ny,nx
      nint = nx * ny
           if (nint .gt. mxint) then
                write(stderr,7) nint, mxint
                goto 9999
           endif
      delx = del_lon * km2lon
      dely = del_lat * km2lat
      xmin = 0.
      ymin = 0.
c
      m=0
      do 41 j = 1,ny
         do 40 i = 1,nx
            m = m + 1
            xi(m) = xmin + delx*float(i-1)
            yi(m) = ymin + dely*float(j-1)
   40 continue
   41 continue
c
c  Do the OA over a gridpoint-by-gridpoint basis
         do k=1,m
c
c  Select the limit # of closest data points 

            call select(limit,xi(k),yi(k),t_cen,x,y,t_since_92,h,
     x       e,eopt,uopt,xopt,topt,most,nobs,dist_lim,time_lim)
c     
c  Perform the OA 
            call objan(uopt,xopt,topt,eopt,nobs,xi(k),yi(k),
     x          tcen,ui(k),erru(k),ier)
c
c  Add back in the average for the field
            ui(k) = ui(k) + aveh
c
         enddo

c
c
c  convert interp locations from distance in km to lat/lon values
      if (irescale .eq. 1) then
        call scale_space (-1,xi,yi,most,lon_min,lat_min,km2lon,km2lat)
      endif

        do k = 1,m
c        write(*,*) xi(k),yi(k),int(ui(k)),int(100*erru(k))
        write(*,*) int(ui(k)),int(100*erru(k))
        enddo
c 
  500 continue
 9999 continue
    7    format(5x,'number of grid points requested', i5,
     / ' exceeds maxium allotted, ',i5,' exiting')
       end


c
c      
c**************************************************************************
c
      subroutine scale_space (idir,x,y,n,lon_min,lat_min,km2lon,km2lat)
c
c    if idir > 0 then scale_space converts lat and long to km 
c       idir < 0                           km to lat and long 
c
c    all km distances are relative to location (lon_min,lat_min)
c
      dimension x(n),y(n)
      real lon_min, lat_min, km2lon, km2lat
c
      if (idir .gt. 0) then
         do i=1,n
           y(i) = km2lat * (y(i)-lat_min)
           x(i) = km2lon * (x(i)-lon_min)
         enddo
      else 
         do i=1,n
           y(i) = lat_min + y(i) / km2lat
           x(i) = lon_min + x(i) / km2lon
         enddo
      endif
c
      return
      end
c      
c
c**************************************************************************
c
      subroutine scale_time (t_sec_92,t_day_92,n)
c     scale time in seconds since 1992 to days since 1992
c
      dimension t_sec_92(n),t_day_92(n)
c
      sec2day = 60.*60.*24.
      do i=1,n
         t_day_92(i) = t_sec_92(i) / sec2day 
      enddo
c
      return
      end
c
c**************************************************************************
c
      real function corr_model (x,y,t)
c
      common /corr/ imodel, c1, c2, c3, c4, km2lon, km2lat
c
      pi=4.0*atan(1.)

      if (imodel .eq. 1 ) then
         r2 = (x/c1)**2 + (y/c2)**2 + (t/c3)**2 
         corr_model = (1-r2) * exp(-r2/2)
      elseif (imodel .eq. 2) then
         r2 = (x/c1)**2 + (y/c2)**2 + (t/c3)**2 
         corr_model = exp(-r2/2)
      elseif (imodel .eq. 3) then
         r2 = (x/c1)**2 + (y/c2)**2 + (t/c3)**2 
         corr_model = exp(-sqrt(r2))
      elseif (imodel .eq. 4) then
         r2 = ( (x-c1*t)**2 + (y-c2*t)**2 ) / (c3*c3)
         corr_model = (1-r2) * exp(-c4 * r2)
      elseif (imodel .eq. 5) then
         r = sqrt(x*x + y*y)
         corr_model = exp(-r/c1)*cos(2.0*pi*r/c2)*exp(-abs(t)/c3)
      elseif (imodel .eq. 6) then
          x=x-c4*t
          r = sqrt(x*x + y*y)
          corr_model = exp(-r/c1)*cos(2.0*pi*r/c2)*exp(-abs(t)/c3)
      endif
      return
      end
c
c**************************************************************************
c
      subroutine remav(psi,m,aver,sdv)
c  routine to calculate the mean and variance of an array and remove 
c  value 
      dimension psi(m)

      aver=0.
      sdv=0.
      do i=1,m
         aver=aver+psi(i)
         sdv=sdv+psi(i)**2
      enddo
      aver=aver/float(m)
      sdv=(sdv/float(m))-aver**2
      if (m .ne. 1) sdv=(float(m)/float(m-1))*sdv
c
      do i=1,m
         psi(i)=psi(i)-aver
      enddo
c
      return
      end
c
c**************************************************************************
c
      subroutine select (limit, xi, yi, tcen, x, y, t, psi,
     1         e,ei, parsi, parx, tobs, n, nobs, dist, tim)
c ***
c *** this routine has been heavily modified so that observations 
c *** are "quasi-independent"- less correlated with each other than with 
c *** the interpolation point.  This was implemented in the paper
c *** "Synoptic Analysis and Dynamical Adjustments of GOES 3 and Seasat
c *** Altimeter Eddy Fields in the Northwest Atlantic.
c *** De Mey and Menard 1989 JGR Vol 94 C5 6221-6230.
c ***
c *** Not sure if this is what they had in mind, but the inversion
c *** routine is much happier.
c ***/
c
c  routine to select the limit nearby data values to interpolation 
c      point xi,yi,tcen
c  limit is the maximum number of points to use
c  dist is the spatial radius of influential points in km
c  tim is the temporal radius of influential points in days
c
      include 'parameter.const'
c
      dimension x(mxobs),y(mxobs),t(mxobs),psi(mxobs),parsi(mxlim),
     >          tobs(mxlim), parx(2,mxlim),e(mxobs),ei(mxlim)
      dimension index(mxobs),cor(mxobs)
      common    /cblock/ c(mxlim)
      common    /corr/   imodel, c1, c2, c3, c4, km2lon, km2lat
      external  corr_model

c /* good will contain the indicies of observations to be used.
c  * this index will refer to the original arrays of mxobs length.
c  */ 
      integer good(mxlim)
c
      nobs=0
      do 50 i=1,n
         delx = xi-x(i)
         dely = yi-y(i)
         delt = tcen-t(i)
         r=sqrt(delx**2+dely**2)
         if (abs(delt) .gt. tim) goto 50
         if (r .gt. dist) goto 50
           nobs        = nobs+1
           index(nobs) = i
           cor(nobs)   = corr_model(delx,dely,delt)
   50 continue

      if (nobs .eq. 0) goto 75
c *
c * For all the observations-
c * 1) sort them by correlation with the gridpoint
c * 2) loop through the reordered observations, start with the most
c *    correlated, excluding observations that are more correlated with 
c *    the current observation, than the current observation is with the
c *    gridpoint.  
c * 3) stop when limit observations are found or out of observations within
c *    the search space
c *
c *  this should ensure that influential points are quasi-independent-
c *  less correlated with each other than to the interpolation point.
c *   this idea is from DeMey and Menard 1989 JGR 6225.
c */

c
c * sort does a sort in place.  It also sorts the index to the 
c * observations that are within the search region (distance limit and 
c * time limit).  abs_cor sorts by ordering abs(cor).  This takes into 
c * account the possiblity of strong negative correlation (negative 
c * lobes.
c */
         call abs_sort(cor,index,nobs)

         j=0  !j count number of usable observations

         do 80 i=1,nobs
            if (cor(i).ne.-999) then  !test if pt has been excluded
               j=j+1
               c(j)=cor(i)
               good(j)=index(i)
               if (j.eq.limit) goto 81
               do 90 k=(i+1),nobs
                  l=index(i)
                  m=index(k)
                  delx=x(l)-x(m)
                  dely=y(l)-y(m)
                  delt=t(l)-t(m)
                  if (abs(corr_model(delx,dely,delt))
     /                              .gt. abs(cor(i))) then
                     cor(k)=-999 !exclude the observation
                  endif
  90            continue
            endif               
  80      continue

  81      continue
         nobs=j
         do 100 i=1,nobs
            parx(1,i)=x(good(i))
            parx(2,i)=y(good(i))
            tobs(i)=t(good(i))
            parsi(i)=psi(good(i))
            ei(i)=e(good(i))
 100     continue

  75      continue
       return
       end
c     
c**************************************************************************
c
c**************************************************************************
c
      subroutine sort(cor,index,n)
c     a shell sort routine to sort index and cor down
c     according to the values of cor
      dimension cor(n),index(n)
      igap=n
    5 if (igap .le. 1) return
      igap=igap/2
      imax=n-igap
   10 iex=0
      do 20 i=1,imax
         iplusg=i+igap
         if (cor(i) .ge. cor(iplusg)) goto 20
         save=cor(i)
         cor(i)=cor(iplusg)
         cor(iplusg)=save
         isave=index(i)
         index(i)=index(iplusg)
         index(iplusg)=isave
         iex=1
   20 continue
      if (iex .ne. 0) goto 10
      goto 5
      end
c
c*********************************************************************
c
c*******************************************************************
      subroutine abs_sort(cor,index,n)
c     a shell sort routine to sort index and cor down
c     according to the values of abs(cor)
      dimension cor(n),index(n)
      igap=n
    5 if (igap .le. 1) return
      igap=igap/2
      imax=n-igap
   10 iex=0
      do 20 i=1,imax
         iplusg=i+igap
         if (abs(cor(i)) .ge. abs(cor(iplusg))) goto 20
         save=cor(i)
         cor(i)=cor(iplusg)
         cor(iplusg)=save
         isave=index(i)
         index(i)=index(iplusg)
         index(iplusg)=isave
         iex=1
   20 continue
      if (iex .ne. 0) goto 10
      goto 5
      end
c
c**************************************************************************
c
      subroutine objan(psi,l,t,eopt,n,x,y,tcen,b,w,ier)
c
c     the space-time objective analysis routine
c     version for 1 interpolation point
c     uses 2 space and 1 time dimension
c     note     delta t = tcen - t(j)
c     l is the array of observation positions, in km
c     t is the time of observation, in days
c     x is the array of interpolation positions, in km
c     tcen is the central interpolation time
c     psi is the array of observation values
c     b is the interpolated value
c     w is the interpolation error limit
c     n is the number of observation points
c     ier is an error flag, zero for no error
c                           -1     no data (warning)
c                           -3     poor matrix inversion (warning)
      include 'parameter.const'
c
      dimension psi(mxlim),t(mxlim),l(2,mxlim),eopt(mxlim)
      real*8    a(mxlim,mxlim)
      real l
      common   /cblock/ c(mxlim)
      external corr_model
      b=0.
      w=1.0
      ier=-1
      if (n .le. 0) goto 500
c  calculate the inverted autocorrelation matrix for the observations
      call seta(a,l,t,n,ier,eopt)

      call getave(a,psi,n,ave)
c  calculate the matrix c
c  -- already calculated in this version, common block cblock
c
c  calculate the rms interpolation error,w
c  and calculate the interpolated value b
      w=0.
      b=0.0
      do 150 i=1,n
         h=0.0
         dumc=c(i)
         do 140 j=1,n
            p=dumc*c(j)*sngl(a(i,j))
            w=w+p
            p2=sngl(a(i,j))*psi(j)
            h=h+p2
  140    continue
         dumy=dumc*h
         b=b+dumy
  150 continue
      b=b+ave
      w=abs(1.-w)
  500 continue
      return
      end
c
c********************************************************************
c
      subroutine seta(a,parx,t,nobs,ier,e)
c
c     this routine calculates the autocorrelation matrix for the
c     observations given the positions, parx and times, t
c     it returns the inverted matrix
c
      include 'parameter.const'
c
      dimension parx(2,mxlim),t(mxlim),e(mxlim)
      real*8 a(mxlim,mxlim),det
      integer ip(2*mxlim)
      external corr_model
      na = mxlim

      do 20 i=1,nobs
         do 10 j=i,nobs
            delt=t(i)-t(j)
            delx=parx(1,i)-parx(1,j)
            dely=parx(2,i)-parx(2,j)
            a(i,j)=dble(corr_model(delx,dely,delt))
            a(j,i)=a(i,j)
   10    continue
       a(i,i)=a(i,i)+dble(e(i))
   20 continue
c     invert the nobs*nobs matrix a
      call invmtx2(a,na,nobs)
      return
      end
c
c**************************************************************************
c
      subroutine getave(a,psi,n,ave)
c     calculate and remove the weighted mean
c
      include 'parameter.const'
c
      dimension psi(n)
      dimension c(mxlim),d(mxlim)
      real*8    a(mxlim,mxlim)
      do 10 i=1,n
         c(i)=0
         d(i)=0
         do 10 k=1,n
            c(i)=c(i)+a(i,k)*psi(k)
            d(i)=d(i)+a(i,k)
   10 enddo
      sum1=0
      sum2=0
      do 20 i=1,n
         sum1=sum1+c(i)
         sum2=sum2+d(i)
   20 enddo
      ave=sum1/sum2
      do 30 i=1,n
           psi(i)=psi(i)-ave
   30 enddo
      return
      end



