SUBROUTINE barnes_cntrl(nx,ny,mxstn,nvar,maxrng,nobs,x,y,               & 2,4
           obs,odiff,obanx,xsta,ysta,isrc,qobs,rely,                    &
           nm_var,npass,kappa,irngsel,range,rngobs,                     &
           wlim,klim,knt,rngwgt,wgtsum,zsum,istat,                      &
           anx,istatus)
!
!  Routine to control Barnes analysis looping.
!  For each pass it calls routines to analyze at grid points,
!  analyze at obs points, computes new obs differences and
!  reports rms statistics.
!
!  Keith Brewster, CAPS
!  March, 1994
!
  IMPLICIT NONE
  INTEGER :: nx,ny,mxstn,nvar,maxrng
  INTEGER :: istat
!
!  Grid arrays
!
  REAL :: x(nx),y(ny)
  REAL :: anx(nx,ny,nvar)
!
!  Observation arrays
!
  INTEGER :: nobs
  REAL :: obs(mxstn,nvar)
  REAL :: odiff(mxstn,nvar)
  REAL :: obanx(mxstn,nvar)
  REAL :: xsta(mxstn),ysta(mxstn)
  REAL :: qobs(mxstn,nvar)
  REAL :: rngobs(mxstn,maxrng)
  INTEGER :: isrc(mxstn)
  INTEGER :: rely(mxstn,nvar)
  CHARACTER (LEN=6) :: nm_var(nvar)
!
!  Barnes control variables
!
  INTEGER :: npass
  INTEGER :: irngsel(nvar)
  REAL :: kappa(npass)
  REAL :: range(nx,ny,maxrng)
  REAL :: wlim
  INTEGER :: klim(npass)
!
!  Barnes scratch space
!
  INTEGER :: knt(nvar)
  REAL :: rngwgt(maxrng)
  REAL :: wgtsum(nvar)
  REAL :: zsum(nvar)
!
  INTEGER :: istatus
!
!  Misc internal variables
!
  INTEGER :: i,j,k,ipass,ista,irng
  INTEGER :: iloc,jloc
  REAL :: delx,dely,rpass
!
!  Initialize rngobs, the appropriate analysis scale range at
!  each observation point.  It is important that this be
!  consistant with the ranges assigned at the grid points.
!  Bilinear interpolation is used rather than an independent
!  calculation.  Obs outside the domain use the nearest grid
!  point in the domain.
!
  delx=x(2)-x(1)
  dely=y(2)-y(1)
  DO ista=1,nobs
    iloc=INT((xsta(ista)-x(1))/delx) + 1
    iloc=MIN(iloc,nx-1)
    iloc=MAX(iloc,1)
    jloc=INT((ysta(ista)-y(1))/dely) + 1
    jloc=MIN(jloc,ny-1)
    jloc=MAX(jloc,1)
    DO irng=1,maxrng
      CALL bilin1(nx,ny,1,xsta(ista),ysta(ista),iloc,jloc,              &
                x,y,range(1,1,irng),rngobs(ista,irng))
    END DO
  END DO
!
!  Establish initial condition.   Analysis is zero and
!  ob differences are equal to obs.
!
  DO k=1,nvar
    DO j=1,ny
      DO i=1,nx
        anx(i,j,k)=0.
      END DO
    END DO
  END DO
!
  DO k=1,nvar
    DO i=1,mxstn
      obanx(i,k)=0.
      odiff(i,k)=obs(i,k)
    END DO
  END DO
!
  DO ipass=1,npass
!
    rpass=kappa(ipass)*kappa(ipass)
!
    CALL bargrid(nx,ny,mxstn,nvar,maxrng,x,y,                           &
                   odiff,xsta,ysta,qobs,isrc,nobs,                      &
                   knt,rngwgt,wgtsum,zsum,                              &
                   rpass,irngsel,range,wlim,klim(ipass),                &
                   anx,istatus)
!
    CALL barobs(mxstn,nvar,maxrng,                                      &
                  odiff,xsta,ysta,qobs,isrc,nobs,                       &
                  knt,rngwgt,wgtsum,zsum,                               &
                  rpass,irngsel,rngobs,wlim,klim(ipass),                &
                  obanx,istatus)
!
    CALL diffnstat(mxstn,nvar,nobs,obs,isrc,obanx,odiff,                &
                       knt,wgtsum,zsum)
!
    WRITE(6,820) ipass
    820   FORMAT('  Statistics for Barnes pass ',i4,/                   &
                 '  ivar  name    knt      bias        rms')
    DO k=1,nvar
      WRITE(6,822) k,nm_var(k),knt(k),wgtsum(k),zsum(k)
      822     FORMAT(i6,2X,a6,i8,f10.3,f10.3)
    END DO
    IF(istat > 0) THEN
      WRITE(istat,820) ipass
      DO k=1,nvar
        WRITE(istat,822) k,nm_var(k),knt(k),wgtsum(k),zsum(k)
      END DO
    END IF
!
  END DO
  istatus=0
  RETURN
END SUBROUTINE barnes_cntrl
!


SUBROUTINE diffnstat(mxstn,nvar,nobs,obs,isrc,obanx,odiff,              & 1
           knt,dmean,drms)
!
!  Routine to calculate new observation differences
!  and report statistics on the differences.
!
!  Keith Brewster, CAPS
!  March, 1994
!
  IMPLICIT NONE
  INTEGER :: mxstn,nvar
!
  REAL :: obs(mxstn,nvar)
  INTEGER :: isrc(mxstn)
  REAL :: obanx(mxstn,nvar)
  REAL :: odiff(mxstn,nvar)
  INTEGER :: nobs
  INTEGER :: knt(nvar)
  REAL :: dmean(nvar)
  REAL :: drms(nvar)
!
!  Misc internal variables
!
  INTEGER :: i,k
  REAL :: flknt
!
  DO k=1,nvar
    knt(k)=0
    dmean(k)=0.
    drms(k)=0.
  END DO
!
  DO i=1,nobs
    IF(isrc(i) /= 0) THEN
      DO k=1,nvar
        IF(obs(i,k) > -90.) THEN
          odiff(i,k)=obs(i,k)-obanx(i,k)
          knt(k)=knt(k)+1
          dmean(k)=dmean(k)+odiff(i,k)
          drms(k)=drms(k)+odiff(i,k)*odiff(i,k)
        END IF
      END DO
    END IF
  END DO
!
!  Calculate stats from the sums formed above
!
  DO k=1,nvar
    IF(knt(k) > 0) THEN
      flknt=FLOAT(knt(k))
      dmean(k)=dmean(k)/flknt
      drms(k)=SQRT(drms(k)/flknt)
    ELSE
      dmean(k)=0.
      drms(k)=0.
    END IF
  END DO
!
  RETURN
END SUBROUTINE diffnstat


SUBROUTINE barqc(maxsta,nvar,nsrc,nsta,z,zdiff,                         & 1
           obstime,xsta,ysta,isrc,                                      &
           knt,wgtsum,zsum,sqsum,                                       &
           range,wlim,klim,qclim,varnam,snam,                           &
           iqclist,iqcsave,                                             &
           flagz,flagd)
!
!   Routine to objectively quality control data using a
!   Barnes analysis of nearby stations.
!
!   If an observations is found to be bad it is set to missing flag
!   flagz, and the corresponding difference variable is set to flagz.
!
!   Keith Brewster, April, 1991
!   Based on Barnes analysis    Keith Brewster, March, 1989
!
  IMPLICIT NONE
!
!  Observation Arguments
!
  INTEGER :: maxsta,nvar,nsrc
  REAL :: z(maxsta,nvar),     & ! original variables
   zdiff(maxsta,nvar)      ! variable to analyze
  INTEGER :: obstime(maxsta)
  REAL :: xsta(maxsta),ysta(maxsta)  ! x and y locations of stations
  INTEGER :: isrc(maxsta)
  INTEGER :: nsta
!
!  Scratch space
!
  INTEGER :: knt(nvar)
  REAL :: wgtsum(nvar),zsum(nvar),sqsum(nvar)
!
!  File pointers
!
  INTEGER :: iqclist,iqcsave
!
!  Analysis specification arguments
!
!  INTEGER iwgt
  REAL :: range,    & ! e-folding range km**2 of barnes weight
   wlim      ! limit of weight to set max range
  INTEGER :: klim   ! minimum # of stations to influence grid pt
  REAL :: qclim(nsrc,nvar)   ! QC cutoffs
!
!  Diagnostics and other argument variables
!
  CHARACTER (LEN=6) :: varnam(nvar)
  CHARACTER (LEN=5) :: snam(maxsta)
  REAL :: flagz,flagd
!
!  Misc internal variables
!
  REAL :: rlimsq,dxsta,dysta,dist,                                      &
       wgt,zanx,dob,sqanx,chklim
  INTEGER :: jsta,ksta,ivar
  LOGICAL :: listit,saveit
!
!  Set printing logicals
!
  listit=(iqclist > 0)
  saveit=(iqcsave > 0)
!
!  Based on the minimum weight to consider, set
!  a maximum range
!
  rlimsq=-range*ALOG(wlim)
!
!  print *, ' Analyzing...'
!  print *, '    rlim= ', rlim
!  print *, '    range = ',range
!
!   Uses Barnes weighting function.  see variable wgt.
!
!   Ksta is location that is being examined
!
  DO ksta=nsta,1,-1
    DO ivar=1,nvar
      zsum(ivar)=0.
      sqsum(ivar)=0.
      wgtsum(ivar)=0.
      knt(ivar)=0
    END DO
    DO jsta=1,nsta
      IF(jsta /= ksta) THEN
        dxsta=xsta(ksta)-xsta(jsta)
        dysta=ysta(ksta)-ysta(jsta)
        dist=dxsta*dxsta + dysta*dysta
        IF(dist < rlimsq) THEN
          wgt=EXP(-dist/range)
          DO ivar=1,nvar
            IF(zdiff(jsta,ivar) > flagd) THEN
              knt(ivar)=knt(ivar)+1
              wgtsum(ivar)=wgtsum(ivar)+wgt
              zsum(ivar)=zsum(ivar)+wgt*zdiff(jsta,ivar)
              sqsum(ivar)=sqsum(ivar)+                                  &
                   wgt*zdiff(jsta,ivar)*zdiff(jsta,ivar)
            END IF
          END DO
        END IF
      END IF  ! dont use current station
    END DO
    DO ivar=1,nvar
      IF(zdiff(ksta,ivar) > flagd) THEN
        IF(knt(ivar) > klim) THEN
          sqanx=sqsum(ivar)-                                            &
                (zsum(ivar)*zsum(ivar)/wgtsum(ivar))
          sqanx=AMAX1(sqanx,0.)
          sqanx=3.0*SQRT(sqanx/wgtsum(ivar))
          zanx=zsum(ivar)/wgtsum(ivar)
          dob=zdiff(ksta,ivar)-zanx
          chklim=AMAX1(sqanx,qclim(isrc(ksta),ivar))
          IF(ABS(dob) > chklim) THEN
            PRINT *, ' QC flagging data...sta=',snam(ksta)
            PRINT *, '     var, ob =',                                  &
                        varnam(ivar),z(ksta,ivar)
            PRINT *, '     delta ob =',dob
            PRINT *, '  3.0*stdv,qclimit=',sqanx,                       &
                                 qclim(isrc(ksta),ivar)
            IF(listit) WRITE(iqclist,810)                               &
                snam(ksta),obstime(ksta),ivar,varnam(ivar),             &
                z(ksta,ivar),dob,                                       &
                sqanx,qclim(isrc(ksta),ivar)
            810           FORMAT(2X,a5,i5,2X,i3,1X,a6,f11.2,f11.2,f11.2,f11.2)
            z(ksta,ivar)=flagz
            zdiff(ksta,ivar)=flagd
          END IF
          zsum(ivar)=dob
        ELSE
!          print *, ' Not enuf data to check ivar= ',
!    +                   ivar,' for ',snam(ksta)
          zsum(ivar)=-999.
        END IF
      ELSE
        zsum(ivar)=-999.
      END IF  ! data missing, dont bother checking
    END DO
    IF(saveit) WRITE(iqcsave,820) snam(ksta),obstime(ksta),             &
              (zdiff(ksta,ivar),ivar=1,nvar)
    820   FORMAT(2X,a4,i6,2X,15(f7.2))
  END DO
  RETURN
END SUBROUTINE barqc


SUBROUTINE barobs(mxstn,nvar,maxrng,                                    & 1
           odiff,xsta,ysta,qobs,isrc,nobs,                              &
           knt,rngwgt,wgtsum,zsum,                                      &
           rpass,irngsel,rngobs,wlim,klim,                              &
           obanx,istatus)
!
!  Routine to objectively analyze data at the observation locations
!  Purpose is to create an analysis parallel to the gridded analysis
!  for obtaining a new observation difference for the next pass.
!  Barnes weight function is employed.
!
!  Keith Brewster, March, 1989
!  Streamlined to remove unneeded options, kb, April, 1990
!  Modified for use in OLAPS  Keith Brewster, March, 1994
!
  IMPLICIT NONE
!
!  Observation Arguments
!
  INTEGER :: mxstn,nvar,maxrng
  REAL :: odiff(mxstn,nvar),             & ! variable to analyse
   xsta(mxstn),ysta(mxstn),  & ! x and y locations of stations
   qobs(mxstn,nvar)
  INTEGER :: isrc(mxstn)
  INTEGER :: nobs
!
!    Scratch space
!
  INTEGER :: knt(nvar)
  REAL :: rngwgt(maxrng)
  REAL :: wgtsum(nvar),zsum(nvar)
!
!  Input/Output analysis at observation locations.Grid arguments
!
  REAL :: obanx(mxstn,nvar)    ! analysed data --  Input and OUTPUT
  INTEGER :: istatus
!
!  Analysis specification arguments
!
!  integer iwgt
  INTEGER :: irngsel(nvar)
  REAL :: rpass,                                                        &
       rngobs(mxstn,maxrng),    & ! e-folding range km**2 of barnes weight
   wlim      ! limit in km of station influence
  INTEGER :: klim   ! minimum # of stations to influence grid pt
!
!  Misc internal variables
!
  INTEGER :: ista,jsta,ivar,irng
  REAL :: denom,rlimsq,dxsta,dysta,dist,wgt,rmax
!  real wvar
!
!  print *, ' Analyzing...'
!  print *, '    rlim= ', rlim
!  print *, '    range = ',range
!
!   Uses Barnes weighting function.  see variable wgt.
!
  DO ista=1,nobs
    IF(isrc(ista) /= 0) THEN
      rmax=rngobs(ista,1)
      DO irng=1,maxrng
        rmax=AMAX1(rmax,rngobs(ista,irng))
      END DO
      denom=rpass*rmax
      rlimsq=-denom*ALOG(wlim)
      DO ivar=1,nvar
        zsum(ivar)=0.
        wgtsum(ivar)=0.
        knt(ivar)=0
      END DO
      DO jsta=1,nobs
        IF(isrc(jsta) /= 0) THEN
          dxsta=xsta(ista)-xsta(jsta)
          dysta=ysta(ista)-ysta(jsta)
          dist=dxsta*dxsta + dysta*dysta
          IF(dist < rlimsq) THEN
            DO irng=1,maxrng
              rngwgt(irng)=                                             &
                     EXP(-dist/(rpass*rngobs(ista,irng)))
            END DO
            DO ivar=1,nvar
              wgt=rngwgt(irngsel(ivar))
              IF(odiff(jsta,ivar) > -90.) THEN
                knt(ivar)=knt(ivar)+1
!                wvar=wgt/qobs(jsta,ivar)
                wgtsum(ivar)=wgtsum(ivar)+wgt
                zsum(ivar)=zsum(ivar)+wgt*odiff(jsta,ivar)
              END IF
            END DO
          END IF
        END IF
      END DO
      DO ivar=1,nvar
        IF(knt(ivar) > klim)                                            &
            obanx(ista,ivar)=obanx(ista,ivar)+zsum(ivar)/wgtsum(ivar)
      END DO
    END IF  ! valid station?
  END DO
  istatus=0
  RETURN
END SUBROUTINE barobs


SUBROUTINE bargrid(nx,ny,mxstn,nvar,maxrng,x,y,                         & 1
           obs,xsta,ysta,qobs,isrc,nobs,                                &
           knt,rngwgt,wgtsum,zsum,                                      &
           rpass,irngsel,range,wlim,klim,                               &
           anx,istatus)
!
!   Routine to objectively analyze data on a regular grid specified
!   by x and y.  Barnes weight function is employed.
!
!   Keith Brewster, March, 1989
!   Streamlined to remove unneeded options, kb, April, 1990
!   Modified for use in OLAPS  Keith Brewster, March, 1994
!
  IMPLICIT NONE
!
!  Observation Arguments
!
  INTEGER :: mxstn,nvar,maxrng
  REAL :: obs(mxstn,nvar),             & ! variable to analyse
   xsta(mxstn),ysta(mxstn),  & ! x and y locations of stations
   qobs(mxstn,nvar)
  INTEGER :: isrc(mxstn)
  INTEGER :: nobs
!
!  Scratch space
!
  INTEGER :: knt(nvar)
  REAL :: rngwgt(maxrng)
  REAL :: wgtsum(nvar),zsum(nvar)
!
!  Grid arguments
!
  INTEGER :: nx,ny
  REAL :: x(nx),y(ny)        ! x and y locations of grid points
  REAL :: anx(nx,ny,nvar)    ! analysed grid --  Input and OUTPUT
  INTEGER :: istatus
!
!  Analysis specification arguments
!
!  integer iwgt
  INTEGER :: irngsel(nvar)
  REAL :: rpass,                                                        &
       range(nx,ny,maxrng),  & ! e-folding range km**2 of barnes weight
   wlim                  ! limit in km of station influence
  INTEGER :: klim   ! minimum # of stations to influence grid pt
!
!  Misc internal variables
!
  REAL :: rlimsq,denom,dxsta,dysta,dist,wgt,rmax
  INTEGER :: i,j,jsta,ivar,irng
!
!   Uses Barnes weighting function.  see variable wgt.
!
  DO j=1,ny
    DO i=1,nx
      rmax=range(i,j,1)
      DO irng=1,maxrng
        rmax=AMAX1(rmax,range(i,j,irng))
      END DO
      denom=rpass*rmax
      rlimsq=-denom*ALOG(wlim)
      DO ivar=1,nvar
        zsum(ivar)=0.
        wgtsum(ivar)=0.
        knt(ivar)=0
      END DO
      DO jsta=1,nobs
        IF(isrc(jsta) /= 0) THEN
          dxsta=x(i)-xsta(jsta)
          dysta=y(j)-ysta(jsta)
          dist=dxsta*dxsta + dysta*dysta
          IF(dist < rlimsq) THEN
            DO irng=1,maxrng
              rngwgt(irng)=EXP(-dist/(rpass*range(i,j,irng)))
            END DO
            DO ivar=1,nvar
              wgt=rngwgt(irngsel(ivar))
              IF(obs(jsta,ivar) > -90.) THEN
                knt(ivar)=knt(ivar)+1
                wgtsum(ivar)=wgtsum(ivar)+wgt
                zsum(ivar)=zsum(ivar)+wgt*obs(jsta,ivar)
              END IF
            END DO
          END IF
        END IF
      END DO
      DO ivar=1,nvar
        IF(knt(ivar) > klim) anx(i,j,ivar)=anx(i,j,ivar)+zsum(ivar)/wgtsum(ivar)
      END DO
    END DO
  END DO
  istatus=0
  RETURN
END SUBROUTINE bargrid


SUBROUTINE bilin1(nx,ny,nvar,xpos,ypos,i,j,x,y,var,varint) 1
!
!  Does bilinear interpolation of variables (nvar variables)
!  to position (xpos,ypos) given a grid of variables "var", their
!  x and y positions of the grid (x,y) and the index of the
!  grid point (i,j) which is to the "lower-left" of xpos and ypos.
!  That means (xpos,ypos) is between (i,j) and (i+1,j+1).
!
!  A saftey feature in case xpos and ypos are outside the grid
!  is built-in.
!
  IMPLICIT NONE
!
!  Arguments, input
!
  INTEGER :: nx,ny,nvar
  REAL :: var(nx,ny,nvar)
  REAL :: x(nx),y(ny),xpos,ypos
  INTEGER :: i,j
!
!  Arguments output
!
  REAL :: varint(nvar)
!
!  Misc Internal Variables
!
  REAL :: xtem,ytem,dxpos,dypos,dx,dy
  REAL :: w11,w12,w21,w22
  INTEGER :: ivar,ip1,jp1
!
  xtem=AMAX1(xpos,x(1))
  xtem=AMIN1(xtem,x(nx))
  ytem=AMAX1(ypos,y(1))
  ytem=AMIN1(ytem,y(ny))
  ip1=i+1
  jp1=j+1
  dx=x(ip1)-x(i)
  dy=y(jp1)-y(j)
  dxpos=xtem-x(i)
  dypos=ytem-y(j)
  IF((dx /= 0.) .AND. (dy /= 0.))THEN
    dx=dxpos/dx
    dy=dypos/dy
    w11=(1.-dx)*(1.-dy)
    w21=dx*(1.-dy)
    w12=dy*(1.-dx)
    w22=dx*dy
  ELSE IF(dx /= 0.) THEN
    dx=dxpos/dx
    w11=1.-dx
    w21=dx
    w12=0.
    w22=0.
  ELSE
    dy=dypos/dy
    w11=1.-dy
    w21=0.
    w12=dy
    w22=0.
  END IF
!
  DO ivar=1,nvar
    varint(ivar)=var(  i,  j,ivar)*w11                                  &
                +var(ip1,  j,ivar)*w21                                  &
                +var(  i,jp1,ivar)*w12                                  &
                +var(ip1,jp1,ivar)*w22
  END DO
  RETURN
END SUBROUTINE bilin1