!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ANXITER                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE anxiter(nx,ny,nz,                                            & 1,12
           nvar,nvarradin,nvarrad,nzua,nzrdr,nzret,                     &
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,                         &
           mxpass,npass,iwstat,xs,ys,zs,icatg,xcor,nam_var,             &
           xsng,ysng,hgtsng,thesng,                                     &
           obsng,odifsng,qobsng,qualsng,isrcsng,icatsng,nobsng,         &
           xua,yua,hgtua,theua,                                         &
           obsua,odifua,qobsua,qualua,isrcua,nlevsua,nobsua,            &
           elvrad,xradc,yradc,                                          &
           distrad,uazmrad,vazmrad,hgtradc,theradc,                     &
           obsrad,odifrad,qobsrad,qualrad,                              &
           irad,isrcrad,nlevrad,ncolrad,                                &
           xretc,yretc,hgtretc,theretc,                                 &
           obsret,odifret,qobsret,qualret,                              &
           iret,isrcret,nlevret,ncolret,                                &
           srcsng,srcua,srcrad,srcret,                                  &
           ianxtyp,iusesng,iuseua,iuserad,iuseret,                      &
           xyrange,kpvrsq,wlim,zrange,zwlim,                            &
           thrng,rngsqi,knt,wgtsum,zsum,                                &
           corsng,corua,corrad,corret,                                  &
           oanxsng,oanxua,oanxrad,oanxret,                              &
           anx,tem1,tem2,tem3,ibegin,iend,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Control ARPS analysis successive corrections iteration.
!  For each pass it calls routines to analyze at grid points,
!  analyze at obs points, computes new obs differences and
!  reports rms statistics.
!
!  On input the analyses at grid points and obs points
!  (anx, oanxsng and oanxua) should be initialized.
!  Either as as zero or with a first guess field and
!  that guess field interpolated to obs sites.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  Version for OLAPS surface Keith Brewster, CAPS
!  March, 1994
!
!  3-D ARPS version Keith Brewster, CAPS
!  July, 1995
!
!  MODIFICATION HISTORY:
!  Update for radar data, misc. improvements
!  Jan, 1996 KB
!
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!  July, 1996 KB
!
!  Changes made to anxiter argument list and calls to Bratseth routines
!  to accomodate changes made to speed up Bratseth routines.
!  Jan 5, 1997 KB
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx         Number of grid points in the x-direction (east/west)
!    ny         Number of grid points in the y-direction (north/south)
!    nz         Number of grid points in the vertical
!    nvar       Number of analysis variables
!    nvarradin  Number of variables in the obsrad array
!    nvarrad    Number of variables in the odifrad array
!    nzua       Maximumnumber of levels in UA observations
!    nzrdr      Maximum number of levels in a radar column
!    nzret      Maximum number of levels in a retrieval column
!    mxsng      Maximum number of single level data
!    mxua       Maximum number of upper air observations
!    mxrad      Maximum number of radars
!    mxcolrad   Maximum number of radar columns
!    mxcolret   Maximum number of retrieval columns
!    mxpass     Maximum number of iterations
!    npass      Number of iterations
!    iwstat     Status indicator for writing statistics
!
!    xs         x location of scalar pts (m)
!    ys         y location of scalar pts (m)
!    zs         z location of scalar pts (m)
!
!    xsng       x location of single-level data
!    ysng       y location of single-level data
!    hgtsng     elevation of single-level data
!    thesng     theta (potential temperature) of single-level data
!
!    obsng      single-level observations
!    odifsng    difference between single-level obs and analysis
!    qobsng     normalized observation error
!    qualsng    single-level data quality indicator
!    isrcsng    index of single-level data source
!    nobsng     number of single-level observations
!
!    xua        x location of upper air data
!    yua        y location of upper air data
!    hgtua      elevation of upper air data
!    theua      theta (potential temperature) of upper air data
!
!    obsua      upper air observations
!    odifua     difference between upper air obs and analysis
!    qobsua     normalized observation error
!    qualua     upper air data quality indicator
!    isrcua     index of upper air data source
!    nlevsua    number of levels of data for each upper air location
!    nobsua     number of upper air observations
!
!    anx        Analyzed variables (or first guess)
!    qback      Background (first guess) error
!
!    nradfil    number of radar files
!    fradname   file name for radar dataset
!    srcrad     name of radar sources
!
!    latrad   latitude of radar  (degrees N)
!    lonrad   longitude of radar (degrees E)
!    elvrad   elevation of feed horn of radar (m MSL)
!    xradc    x location of radar column
!    yradc    y location of radar column
!    irad     radar number
!    nlevrad  number of levels of radar data in each column
!    distrad  distance of radar column from source radar
!    uazmrad  u-component of radar beam for each column
!    vazmrad  v-component of radar beam for each column
!    hgtradc  height (m MSL) of radar observations
!    obsrad   radar observations
!    oanxrad  analysis (first guess) value at radar data location
!    odifrad  difference between radar observation and analysis
!    qobsrad  normalized observation error
!    qualrad  radar data quality indicator
!    ncolrad  number of radar columns read-in
!    istatus  status indicator
!
!    latret   latitude of retrieval radar  (degrees N)
!    lonret   longitude of retrieval radar (degrees E)
!    xretc    x location of retrieval column
!    yretc    y location of retrieval column
!    iret     retrieval number
!    nlevret  number of levels of retrieval data in each column
!    hgtretc  height (m MSL) of retrieval observations
!    obsret   retrieval observations
!    odifret  difference between retrieval observation and analysis
!    qobsret  normalized observation error
!    qualret  retrieval data quality indicator
!    ncolret  number of retr columns read-in
!    istatus  status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Input Sizing Arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  INTEGER :: nvar,nvarradin,nvarrad
  INTEGER :: nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
  INTEGER :: mxpass,npass
!
!-----------------------------------------------------------------------
!
!  Input Grid Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: xs(nx)
  REAL :: ys(ny)
  REAL :: zs(nx,ny,nz)
  INTEGER :: icatg(nx,ny)
  REAL :: xcor(ncat,ncat)
!
!-----------------------------------------------------------------------
!
!  Input Observation Arguments
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=6) :: nam_var(nvar)
  REAL :: xsng(mxsng)
  REAL :: ysng(mxsng)
  REAL :: hgtsng(mxsng)
  REAL :: thesng(mxsng)
  REAL :: obsng(nvar,mxsng)
  REAL :: odifsng(nvar,mxsng)
  REAL :: qobsng(nvar,mxsng)
  INTEGER :: qualsng(nvar,mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: icatsng(mxsng)
  INTEGER :: nobsng
!
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  REAL :: theua(nzua,mxua)
  REAL :: obsua(nvar,nzua,mxua)
  REAL :: odifua(nvar,nzua,mxua)
  REAL :: qobsua(nvar,nzua,mxua)
  INTEGER :: qualua(nvar,nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: elvrad(mxrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: distrad(mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: theradc(nzrdr,mxcolrad)
  REAL :: obsrad(nvarradin,nzrdr,mxcolrad)
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: qobsrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  REAL :: theretc(nzret,mxcolret)
  REAL :: obsret(nvar,nzret,mxcolret)
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: qobsret(nvar,nzret,mxcolret)
  INTEGER :: qualret(nvar,nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Input Analysis Control Variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=8) :: srcsng(nsrcsng)
  CHARACTER (LEN=8) :: srcua (nsrcua )
  CHARACTER (LEN=8) :: srcrad(nsrcrad)
  CHARACTER (LEN=8) :: srcret(nsrcret)

  INTEGER :: ianxtyp(mxpass)
  INTEGER :: iusesng(0:nsrcsng,mxpass)
  INTEGER :: iuseua(0:nsrcua,mxpass)
  INTEGER :: iuserad(0:nsrcrad,mxpass)
  INTEGER :: iuseret(0:nsrcret,mxpass)

  REAL :: xyrange(mxpass)
  REAL :: kpvrsq(nvar)
  REAL :: wlim
  REAL :: zrange(mxpass)
  REAL :: zwlim
  REAL :: thrng(mxpass)
  INTEGER :: iwstat
!
!-----------------------------------------------------------------------
!
!  Scratch Space
!
!-----------------------------------------------------------------------
!
  REAL :: rngsqi(nvar)
  INTEGER :: knt(nvar,nz)
  REAL :: wgtsum(nvar,nz)
  REAL :: zsum(nvar,nz)
!
!-----------------------------------------------------------------------
!
!  Output Variables at Observation Locations
!
!-----------------------------------------------------------------------
!
  REAL :: corsng(nvar,mxsng)
  REAL :: corua(nvar,nzua,mxua)
  REAL :: corrad(nvarrad,nzrdr,mxcolrad)
  REAL :: corret(nvar,nzret,mxcolret)
  REAL :: oanxsng(nvar,mxsng)
  REAL :: oanxua(nvar,nzua,mxua)
  REAL :: oanxrad(nvarrad,nzrdr,mxcolrad)
  REAL :: oanxret(nvar,nzret,mxcolret)
!
!-----------------------------------------------------------------------
!
!  Output Grid
!
!-----------------------------------------------------------------------
!
  REAL :: anx(nx,ny,nz,nvar)
!
!-----------------------------------------------------------------------
!
!  Work arrays
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
  REAL :: tem3(nx,ny,nz)
  INTEGER :: ibegin(ny)
  INTEGER :: iend(ny)
!
!-----------------------------------------------------------------------
!
!  Return status
!
!-----------------------------------------------------------------------
!
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  REAL :: ftabinv,setexp
  INTEGER :: i,j,k,ipass,isrc,sngsw,uasw,radsw,retsw
  REAL :: rpass,zrngsq,thrngsq
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Note in the following the analysis scratch arrays
!  are used to hold the output statistics
!
!-----------------------------------------------------------------------
!
  ipass=0
  WRITE(6,'(a,i4)') 'Computing orig difference stats'
  CALL diffnst3d(nvar,                                                  &
                 nzua,nzret,mxsng,mxua,mxcolret,                        &
                 obsng,odifsng,oanxsng,isrcsng,qualsng,nobsng,          &
                 obsua,odifua,oanxua,isrcua,qualua,nlevsua,nobsua,      &
                 obsret,odifret,oanxret,iret,                           &
                 qualret,nlevret,ncolret,                               &
                 knt,wgtsum,zsum)
!
  WRITE(6,805) ipass
  805   FORMAT('  Statistics for analysis pass ',i4,/                   &
               '  ivar  name     knt    bias        rms')
  DO k=1,nvar
    WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
    812     FORMAT(i6,2X,a6,i6,g12.3,g11.3)
  END DO
  CALL diffnstra(nvar,nvarrad,nvarradin,nzrdr,mxrad,mxcolrad,           &
            elvrad,distrad,hgtradc,                                     &
            uazmrad,vazmrad,qualrad,irad,                               &
            obsrad,odifrad,oanxrad,                                     &
            nlevrad,ncolrad,                                            &
            knt,wgtsum,zsum)
  WRITE(6,815)
  815   FORMAT('  Radar data variables ',i4,/                           &
               '  ivar  name    knt      bias        rms')
  DO k=1,nvar
    WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
  END DO
  WRITE(6,806)
  806   FORMAT(/,' ')

  IF(iwstat > 0) THEN
    WRITE(iwstat,820) ipass
    DO k=1,nvarrad
      WRITE(iwstat,812) k,nam_var(k),                                   &
                        knt(k,1),wgtsum(k,1),zsum(k,1)
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Initialize correlation lookup table.
!
!-----------------------------------------------------------------------
!
  ftabinv=setexp(16.)
!
!-----------------------------------------------------------------------
!
!  Loop through npass analysis iterations
!
!-----------------------------------------------------------------------
!
  DO ipass=1,npass
!
!-----------------------------------------------------------------------
!
!  Set single-level usage switch based in iusesng
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/a,i4)') ' Source usage switches for pass ',ipass
    WRITE(6,'(/a,i4)') '   Single level data'
    sngsw=0
    DO isrc=1,nsrcsng
      IF(iusesng(isrc,ipass) > 0) THEN
        sngsw=1
        WRITE(6,'(a,a)') '      Using     ',srcsng(isrc)
      END IF
    END DO
    IF(sngsw == 0) WRITE(6,'(a)') '      none'
!
!-----------------------------------------------------------------------
!
!  Set multiple-level usage switch based in iuseua
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/a,i4)') '   Multiple level data'
    uasw=0
    DO isrc=1,nsrcua
      IF(iuseua(isrc,ipass) > 0) THEN
        uasw=1
        WRITE(6,'(a,a)') '      Using     ',srcua(isrc)
      END IF
    END DO
    IF(uasw == 0) WRITE(6,'(a)') '      none'
!
!-----------------------------------------------------------------------
!
!  Set radar-data usage switch based in iuserad
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/a,i4)') '   Raw radar data'
    radsw=0
    DO isrc=1,nsrcrad
      IF(iuserad(isrc,ipass) > 0) THEN
        radsw=1
        WRITE(6,'(a,a)') '      Using     ',srcrad(isrc)
      END IF
    END DO
    IF(radsw == 0) WRITE(6,'(a)') '      none'
!
!-----------------------------------------------------------------------
!
!  Set ret usage switch based in iuseret
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/a,i4)') '   Retrieval radar data'
    retsw=0
    DO isrc=1,nsrcret
      IF(iuseret(isrc,ipass) > 0) THEN
        retsw=1
        WRITE(6,'(a,a)') '      Using     ',srcret(isrc)
      END IF
    END DO
    IF(retsw == 0) WRITE(6,'(a)') '      none'

    rpass=xyrange(ipass)*xyrange(ipass)
    zrngsq=zrange(ipass)*zrange(ipass)
    thrngsq=thrng(ipass)*thrng(ipass)
!
    IF(ianxtyp(ipass) == 11) THEN

      WRITE(6,'(/a,a,i4)') 'Barnes on grid using height',               &
                        ' ipass=',ipass
      CALL bargrd3d(nx,ny,nz,                                           &
             nvar,nvarrad,nzua,nzrdr,nzret,                             &
             mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                  &
             nsrcsng,nsrcua,nsrcrad,nsrcret,xs,ys,zs,                   &
             odifsng,xsng,ysng,hgtsng,qualsng,isrcsng,nobsng,           &
             odifua,xua,yua,hgtua,qualua,isrcua,nlevsua,nobsua,         &
             odifrad,xradc,yradc,hgtradc,                               &
             uazmrad,vazmrad,qualrad,                                   &
             irad,isrcrad,nlevrad,ncolrad,                              &
             odifret,xretc,yretc,hgtretc,                               &
             qualret,iret,isrcret,nlevret,ncolret,                      &
             iusesng(0,ipass),iuseua(0,ipass),                          &
             iuserad(0,ipass),iuseret(0,ipass),                         &
             sngsw,uasw,radsw,retsw,kpvrsq,rpass,                       &
             wlim,zrngsq,zwlim,                                         &
             rngsqi,knt,wgtsum,zsum,                                    &
             anx,istatus)
!
      WRITE(6,'(/a,a,i4)') 'Barnes at obs sites using height',          &
                        ' ipass=',ipass
      CALL barobs3d(nvar,nvarrad,nzua,nzrdr,nzret,                      &
             mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                  &
             nsrcsng,nsrcua,nsrcrad,nsrcret,                            &
             odifsng,xsng,ysng,hgtsng,qualsng,isrcsng,nobsng,           &
             odifua,xua,yua,hgtua,qualua,isrcua,nlevsua,nobsua,         &
             odifrad,xradc,yradc,hgtradc,uazmrad,vazmrad,               &
             qualrad,irad,isrcrad,nlevrad,ncolrad,                      &
             odifret,xretc,yretc,hgtretc,                               &
             qualret,iret,isrcret,nlevret,ncolret,                      &
             iusesng(0,ipass),iuseua(0,ipass),                          &
             iuserad(0,ipass),iuseret(0,ipass),                         &
             sngsw,uasw,radsw,retsw,kpvrsq,rpass,                       &
             wlim,zrngsq,zwlim,                                         &
             rngsqi,knt,wgtsum,zsum,                                    &
             oanxsng,oanxua,oanxrad,oanxret,istatus)
!
    ELSE IF(ianxtyp(ipass) == 21) THEN
!
      CALL obscor(nvar,nvarrad,nzua,nzrdr,nzret,                        &
             mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                  &
             nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                  &
             qobsng,xsng,ysng,hgtsng,                                   &
             qualsng,isrcsng,icatsng,nobsng,                            &
             qobsua,xua,yua,hgtua,                                      &
             qualua,isrcua,nlevsua,nobsua,                              &
             qobsrad,xradc,yradc,hgtradc,                               &
             uazmrad,vazmrad,                                           &
             qualrad,irad,isrcrad,nlevrad,ncolrad,                      &
             qobsret,xretc,yretc,hgtretc,                               &
             qualret,iret,isrcret,nlevret,ncolret,                      &
             iusesng(0,ipass),iuseua(0,ipass),                          &
             iuserad(0,ipass),iuseret(0,ipass),                         &
             sngsw,uasw,radsw,retsw,kpvrsq,rpass,rngsqi,                &
             wlim,zrngsq,                                               &
             corsng,corua,corrad,corret,istatus)
!
!-----------------------------------------------------------------------
!
!  Calculate covariance sums at observation points for
!
!-----------------------------------------------------------------------
!
      WRITE(6,'(/a,a,i4)') 'Bratseth on grid using height',             &
                        ' ipass=',ipass
      CALL brtgrd3d(nx,ny,nz,nvar,nvarrad,                              &
            nzua,nzrdr,nzret,                                           &
            mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                   &
            nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xs,ys,zs,icatg,xcor,    &
            odifsng,corsng,xsng,ysng,hgtsng,isrcsng,icatsng,nobsng,     &
            odifua,corua,xua,yua,hgtua,isrcua,nlevsua,nobsua,           &
            odifrad,corrad,xradc,yradc,hgtradc,                         &
            uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,               &
            odifret,corret,xretc,yretc,hgtretc,                         &
            iret,isrcret,nlevret,ncolret,                               &
            iusesng(0,ipass),iuseua(0,ipass),                           &
            iuserad(0,ipass),iuseret(0,ipass),                          &
            sngsw,uasw,radsw,retsw,kpvrsq,rpass,                        &
            wlim,zrngsq,rngsqi,                                         &
            anx,tem1,tem2,ibegin,iend,istatus)
!
      WRITE(6,'(a,a,i4)') 'Bratseth at obs sites using height',         &
                        ' ipass=',ipass
      CALL brtobs3d(nvar,nvarrad,nzua,nzrdr,nzret,                      &
            mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                   &
            nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                   &
            odifsng,qobsng,corsng,xsng,ysng,hgtsng,                     &
            isrcsng,icatsng,nobsng,                                     &
            odifua,qobsua,corua,xua,yua,hgtua,isrcua,nlevsua,nobsua,    &
            odifrad,qobsrad,corrad,xradc,yradc,hgtradc,                 &
            uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,               &
            odifret,qobsret,corret,xretc,yretc,hgtretc,                 &
            iret,isrcret,nlevret,ncolret,                               &
            iusesng(0,ipass),iuseua(0,ipass),                           &
            iuserad(0,ipass),iuseret(0,ipass),                          &
            sngsw,uasw,radsw,retsw,kpvrsq,rpass,                        &
            wlim,zrngsq,rngsqi,                                         &
            oanxsng,oanxua,oanxrad,oanxret,istatus)

    ELSE
!
      CALL obscor(nvar,nvarrad,nzua,nzrdr,nzret,                        &
            mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                   &
            nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                   &
            qobsng,xsng,ysng,thesng,                                    &
            qualsng,isrcsng,icatsng,nobsng,                             &
            qobsua,xua,yua,theua,                                       &
            qualua,isrcua,nlevsua,nobsua,                               &
            qobsrad,xradc,yradc,theradc,                                &
            uazmrad,vazmrad,                                            &
            qualrad,irad,isrcrad,nlevrad,ncolrad,                       &
            qobsret,xretc,yretc,theretc,                                &
            qualret,iret,isrcret,nlevret,ncolret,                       &
            iusesng(0,ipass),iuseua(0,ipass),                           &
            iuserad(0,ipass),iuseret(0,ipass),                          &
            sngsw,uasw,radsw,retsw,kpvrsq,rpass,rngsqi,                 &
            wlim,thrngsq,                                               &
            corsng,corua,corrad,corret,istatus)
!
      WRITE(6,'(/a,a,i4)') 'Bratseth on grid using theta',              &
                        ' ipass=',ipass
!
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            tem3(i,j,k)=anx(i,j,k,4)
          END DO
        END DO
      END DO
!
      CALL brtgrd3d(nx,ny,nz,nvar,nvarrad,                              &
            nzua,nzrdr,nzret,                                           &
            mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                   &
            nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,                        &
            xs,ys,tem3,icatg,xcor,                                      &
            odifsng,corsng,xsng,ysng,thesng,                            &
            isrcsng,icatsng,nobsng,                                     &
            odifua,corua,xua,yua,theua,isrcua,nlevsua,nobsua,           &
            odifrad,corrad,xradc,yradc,theradc,                         &
            uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,               &
            odifret,corret,xretc,yretc,theretc,                         &
            iret,isrcret,nlevret,ncolret,                               &
            iusesng(0,ipass),iuseua(0,ipass),                           &
            iuserad(0,ipass),iuseret(0,ipass),                          &
            sngsw,uasw,radsw,retsw,kpvrsq,rpass,                        &
            wlim,thrngsq,rngsqi,                                        &
            anx,tem1,tem2,ibegin,iend,istatus)
!
      WRITE(6,'(a,a,i4)') 'Bratseth at obs sites using theta',          &
                        ' ipass=',ipass
!
!-----------------------------------------------------------------------
!
      CALL brtobs3d(nvar,nvarrad,nzua,nzrdr,nzret,                      &
             mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                  &
             nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                  &
             odifsng,qobsng,corsng,xsng,ysng,thesng,                    &
             isrcsng,icatsng,nobsng,                                    &
             odifua,qobsua,corua,xua,yua,theua,                         &
             isrcua,nlevsua,nobsua,                                     &
             odifrad,qobsrad,corrad,xradc,yradc,theradc,                &
             uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,              &
             odifret,qobsret,corret,xretc,yretc,theretc,                &
             iret,isrcret,nlevret,ncolret,                              &
             iusesng(0,ipass),iuseua(0,ipass),                          &
             iuserad(0,ipass),iuseret(0,ipass),                         &
             sngsw,uasw,radsw,retsw,kpvrsq,rpass,                       &
             wlim,thrngsq,rngsqi,                                       &
             oanxsng,oanxua,oanxrad,oanxret,istatus)
    END IF
!
!-----------------------------------------------------------------------
!
!  Note in the following the analysis scratch arrays
!  are used to hold the output statistics
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(a,i4)') 'Computing new differences, ipass=',ipass
    CALL diffnst3d(nvar,                                                &
                   nzua,nzret,mxsng,mxua,mxcolret,                      &
                   obsng,odifsng,oanxsng,isrcsng,qualsng,nobsng,        &
                   obsua,odifua,oanxua,isrcua,qualua,nlevsua,nobsua,    &
                   obsret,odifret,oanxret,iret,                         &
                   qualret,nlevret,ncolret,                             &
                   knt,wgtsum,zsum)
!
    WRITE(6,820) ipass
    820   FORMAT('  Statistics for analysis pass ',i4,/                 &
                 '  ivar  name    knt      bias        rms')
    DO k=1,nvar
      WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
    END DO
    IF(iwstat > 0) THEN
      WRITE(iwstat,820) ipass
      DO k=1,nvar
        WRITE(iwstat,812) k,nam_var(k),                                 &
                          knt(k,1),wgtsum(k,1),zsum(k,1)
      END DO
      WRITE(iwstat,806)
    END IF
    CALL diffnstra(nvar,nvarrad,nvarradin,nzrdr,mxrad,mxcolrad,         &
              elvrad,distrad,hgtradc,                                   &
              uazmrad,vazmrad,qualrad,irad,                             &
              obsrad,odifrad,oanxrad,                                   &
              nlevrad,ncolrad,                                          &
              knt,wgtsum,zsum)
    WRITE(6,815)
    DO k=1,nvar
      WRITE(6,812) k,nam_var(k),knt(k,1),wgtsum(k,1),zsum(k,1)
    END DO
    WRITE(6,806)
!
  END DO
!
  istatus=0
  RETURN
END SUBROUTINE anxiter
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE DIFFNST3d                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE diffnst3d(nvar,                                              & 2
           nzua,nzret,mxsng,mxua,mxcolret,                              &
           obsng,odifsng,oanxsng,isrcsng,qualsng,nobsng,                &
           obsua,odifua,oanxua,isrcua,qualua,nlevsua,nobsua,            &
           obsret,odifret,oanxret,iret,                                 &
           qualret,nlevret,ncolret,                                     &
           knt,dmean,drms)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate new observation differences
!  and report statistics on the differences.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  March, 1994
!
!  MODIFICATION HISTORY:
!  July, 1995  (Keith Brewster)
!  3D ARPS version
!
!  May, 1996 (KB)
!  Added retrieval data.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nvar
  INTEGER :: nzua,nzret,mxsng,mxua,mxcolret
!
  REAL :: obsng(nvar,mxsng)
  REAL :: odifsng(nvar,mxsng)
  REAL :: oanxsng(nvar,mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: qualsng(nvar,mxsng)
  INTEGER :: nobsng
!
  REAL :: obsua(nvar,nzua,mxua)
  REAL :: odifua(nvar,nzua,mxua)
  REAL :: oanxua(nvar,nzua,mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: qualua(nvar,nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: nobsua
!
  REAL :: obsret(nvar,nzret,mxcolret)
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: oanxret(nvar,nzret,mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: qualret(nvar,nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: ncolret
!
  INTEGER :: knt(nvar)
  REAL :: dmean(nvar)
  REAL :: drms(nvar)
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: ivar,jsta,klev
  REAL :: flknt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO ivar=1,nvar
    knt(ivar)=0
    dmean(ivar)=0.
    drms(ivar)=0.
  END DO
!
!-----------------------------------------------------------------------
!
!  Find differences and form sums at single-level sites
!
!-----------------------------------------------------------------------
!
  DO jsta=1,nobsng
    IF(isrcsng(jsta) /= 0) THEN
      DO ivar=1,nvar
        IF(qualsng(ivar,jsta) > 0) THEN
          odifsng(ivar,jsta)=obsng(ivar,jsta)-oanxsng(ivar,jsta)
          knt(ivar)=knt(ivar)+1
          dmean(ivar)=dmean(ivar)-odifsng(ivar,jsta)
          drms(ivar)=drms(ivar)+                                        &
                     odifsng(ivar,jsta)*odifsng(ivar,jsta)
        END IF
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Find differences at multiple-level sites
!  Add to sums
!
!-----------------------------------------------------------------------
!
  DO jsta=1,nobsua
    IF(isrcua(jsta) /= 0) THEN
      DO klev=1,nlevsua(jsta)
        DO ivar=1,nvar
          IF(qualua(ivar,klev,jsta) > 0) THEN
            odifua(ivar,klev,jsta)=                                     &
                obsua(ivar,klev,jsta)-oanxua(ivar,klev,jsta)
            knt(ivar)=knt(ivar)+1
            dmean(ivar)=dmean(ivar)-odifua(ivar,klev,jsta)
            drms(ivar)=drms(ivar)+                                      &
                odifua(ivar,klev,jsta)*odifua(ivar,klev,jsta)
          END IF
        END DO
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Find differences at retrieval columns
!  Add to sums
!
!-----------------------------------------------------------------------
!
  DO jsta=1,ncolret
    IF(iret(jsta) > 0) THEN
      DO klev=1,nlevret(jsta)
        DO ivar=1,nvar
          IF(qualret(ivar,klev,jsta) > 0) THEN
            odifret(ivar,klev,jsta)=                                    &
                obsret(ivar,klev,jsta)-oanxret(ivar,klev,jsta)
            knt(ivar)=knt(ivar)+1
            dmean(ivar)=dmean(ivar)-odifret(ivar,klev,jsta)
            drms(ivar)=drms(ivar)+                                      &
                odifret(ivar,klev,jsta)*odifret(ivar,klev,jsta)
          END IF
        END DO
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate stats from the sums formed above
!
!-----------------------------------------------------------------------
!
  DO ivar=1,nvar
    IF(knt(ivar) > 0) THEN
      flknt=FLOAT(knt(ivar))
      dmean(ivar)=dmean(ivar)/flknt
      drms(ivar)=SQRT(drms(ivar)/flknt)
    ELSE
      dmean(ivar)=0.
      drms(ivar)=0.
    END IF
  END DO
!
  RETURN
END SUBROUTINE diffnst3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE DIFFNSTRA                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE diffnstra(nvar,nvarrad,nvarradin,nzrdr,mxrad,mxcolrad,       & 2,2
           elvrad,distrad,hgtradc,                                      &
           uazmrad,vazmrad,qualrad,irad,                                &
           obsrad,odifrad,oanxrad,                                      &
           nlevrad,ncolrad,                                             &
           knt,dmean,drms)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate radar observation differences
!  and report statistics on the differences.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  September, 1995
!
!  MODIFICATION HISTORY:
!
!  Jan, 1996 (KB)
!  Added radar data and other improvements.
!
!  August, 2001 (KB)
!  Corrected call to dhdr which now returns dhdr not local
!  elevation angle.  Modified calculation of odifrad to speed
!  convergence to vr.
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nvar,nvarrad,nvarradin,nzrdr,mxrad,mxcolrad
!
  REAL :: elvrad(mxrad)
  REAL :: distrad(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: irad(mxcolrad)
  REAL :: obsrad(nvarradin,nzrdr,mxcolrad)
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: oanxrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: ncolrad
!
  INTEGER :: knt(nvar)
  REAL :: dmean(nvar)
  REAL :: drms(nvar)
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: icol,ivar,kr,jrad
  REAL :: pi,dtr,delz,eleva,range,flknt
  REAL :: vr,vrdiff,avrdif,dhdr,dsdr,dsdrinv
  REAL :: vr_miss
!
!-----------------------------------------------------------------------
!
!  QC and conversion parameters
!
!  vrqcthr   maximum radial velocity diff from background to use
!  vrcnvthr  lower limit of the dot product between the radial
!            direction and the analysis wind component to use data
!
!-----------------------------------------------------------------------
!
  REAL :: vrqcthr,vrcnvthr
  PARAMETER (vrqcthr=20.,   & ! maximum vrdiff from background to use
         vrcnvthr=-0.2) ! lower limit of the dot product to use data
!
!-----------------------------------------------------------------------
!
!  Some initializations
!
!-----------------------------------------------------------------------
!
  vr_miss=-999.
  pi=4.*ATAN(1.0)
  dtr=pi/180.
!
  DO ivar=1,nvar
    knt(ivar)=0
    dmean(ivar)=0.
    drms(ivar)=0.
  END DO
!
!-----------------------------------------------------------------------
!
!  Find radial velocity differences at radar data points
!  Add to sums
!
!-----------------------------------------------------------------------
!
  DO icol=1,ncolrad
    IF(irad(icol) > 0) THEN
      jrad=irad(icol)
!
      DO kr=1,nlevrad(icol)
        IF(qualrad(1,kr,icol) > 0 .AND. qualrad(2,kr,icol) > 0 ) THEN
!
          delz=hgtradc(kr,icol)-elvrad(jrad)
          CALL beamelv(delz,distrad(icol),eleva,range)
          CALL dhdrange(eleva,range,dhdr)
          dsdr=SQRT(MAX(0.,(1.-dhdr*dhdr)))
          IF(dsdr /= 0.) THEN
            dsdrinv=1./dsdr
          ELSE
            dsdrinv=0.
          END IF
!
          vr=(uazmrad(icol)*oanxrad(1,kr,icol) +                        &
              vazmrad(icol)*oanxrad(2,kr,icol)) * dsdr
          vrdiff=obsrad(2,kr,icol)-vr
          avrdif=ABS(vrdiff)

          IF(avrdif < vrqcthr) THEN
!
!-----------------------------------------------------------------------
!
!  Find increment in analysis u and v needed to make vrdiff zero
!
!-----------------------------------------------------------------------
!
            odifrad(1,kr,icol)=dsdrinv*uazmrad(icol)*vrdiff
            odifrad(2,kr,icol)=dsdrinv*vazmrad(icol)*vrdiff
!
            knt(1)=knt(1)+1
            dmean(1)=dmean(1)-odifrad(1,kr,icol)
            drms(1)=drms(1)+                                            &
                odifrad(1,kr,icol)*odifrad(1,kr,icol)

            knt(2)=knt(2)+1
            dmean(2)=dmean(2)-odifrad(2,kr,icol)
            drms(2)=drms(2)+                                            &
                odifrad(2,kr,icol)*odifrad(2,kr,icol)

          ELSE

            odifrad(1,kr,icol)=vr_miss
            qualrad(1,kr,icol)=-199
            odifrad(2,kr,icol)=vr_miss
            qualrad(2,kr,icol)=-199

          END IF
        END IF
        IF(qualrad(3,kr,icol) > 0) THEN
          odifrad(3,kr,icol)=obsrad(1,kr,icol)-oanxrad(3,kr,icol)
!          print *, ' kr,icol,obsrad,oanx= ',kr,icol,
!    +                 obsrad(1,kr,icol),oanxrad(3,kr,icol)
          knt(5)=knt(5)+1
          dmean(5)=dmean(5)-odifrad(3,kr,icol)
          drms(5)=drms(5)+                                              &
              odifrad(3,kr,icol)*odifrad(3,kr,icol)
        END IF
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate stats from the sums formed above
!
!-----------------------------------------------------------------------
!
  DO ivar=1,nvar
    IF(knt(ivar) > 0) THEN
      flknt=FLOAT(knt(ivar))
      dmean(ivar)=dmean(ivar)/flknt
      drms(ivar)=SQRT(drms(ivar)/flknt)
    ELSE
      dmean(ivar)=0.
      drms(ivar)=0.
    END IF
  END DO
!
  RETURN
END SUBROUTINE diffnstra
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BAROBS3d                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE barobs3d(nvar,nvarrad,nzua,nzrdr,nzret,                      & 1
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,                              &
           odifsng,xsng,ysng,hgtsng,qualsng,isrcsng,nobsng,             &
           odifua,xua,yua,hgtua,qualua,isrcua,nlevsua,nobsua,           &
           odifrad,xradc,yradc,hgtradc,uazmrad,vazmrad,                 &
           qualrad,irad,isrcrad,nlevrad,ncolrad,                        &
           odifret,xretc,yretc,hgtretc,                                 &
           qualret,iret,isrcret,nlevret,ncolret,                        &
           iusesng,iuseua,iuserad,iuseret,                              &
           sngsw,uasw,radsw,retsw,kpvrsq,rpass,                         &
           wlim,zrngsq,zwlim,                                           &
           rngsqi,knt,wgtsum,zsum,                                      &
           oanxsng,oanxua,oanxrad,oanxret,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Objectively analyze data at observation pts
!  Barnes weight function is employed.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  July, 1995
!
!  MODIFICATION HISTORY:
!
!  Jan, 1996 (KB)
!  Added radar data and other improvements.
!
!  July, 1996 (KB)
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Sizing arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nvar,nvarrad,nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: odifsng(nvar,mxsng),                                          &
       xsng(mxsng),                                                     &
       ysng(mxsng),                                                     &
       hgtsng(mxsng)
  INTEGER :: qualsng(nvar,mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: nobsng
  REAL :: odifua(nvar,nzua,mxua)             ! variable to analyse
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  INTEGER :: qualua(nvar,nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  INTEGER :: qualret(nvar,nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Scratch space
!
!-----------------------------------------------------------------------
!
  INTEGER :: knt(nvar)
  REAL :: rngsqi(nvar)
  REAL :: wgtsum(nvar)
  REAL :: zsum(nvar)
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: iusesng(0:nsrcsng)
  INTEGER :: iuseua(0:nsrcua)
  INTEGER :: iuserad(0:nsrcrad)
  INTEGER :: iuseret(0:nsrcret)
  INTEGER :: sngsw,uasw,radsw,retsw
  REAL :: kpvrsq(nvar)
  REAL :: rpass
  REAL :: wlim
  REAL :: zrngsq
  REAL :: zwlim
!
!-----------------------------------------------------------------------
!
!  Output arguments
!
!-----------------------------------------------------------------------
!
  REAL :: oanxsng(nvar,mxsng)
  REAL :: oanxua(nvar,nzua,mxua)
  REAL :: oanxrad(nvarrad,nzrdr,mxcolrad)
  REAL :: oanxret(nvar,nzret,mxcolret)
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  REAL :: texp
  REAL :: rlimsq,zlimsq,zsqinv
  REAL :: dist,dzsq,wgt
  INTEGER :: ista,jsta,k,klev,ivar
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rlimsq=-rpass*ALOG(wlim)
  DO ivar=1,nvar
    rngsqi(ivar)=1./(kpvrsq(ivar)*rpass)
  END DO
!
  zlimsq=-zrngsq*ALOG(zwlim)
  zsqinv=1./zrngsq
!
!-----------------------------------------------------------------------
!
!  Loop over single stations.
!
!-----------------------------------------------------------------------
!
  DO ista=1,nobsng
    IF(iusesng(isrcsng(ista)) > 0) THEN
      DO ivar=1,nvar
        knt(ivar)=0
        wgtsum(ivar)=0.
        zsum(ivar)=0.
      END DO
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
      IF(sngsw > 0) THEN
        DO jsta=1,nobsng
          IF(iusesng(isrcsng(jsta)) > 0) THEN
            dist=(xsng(ista)-xsng(jsta))*(xsng(ista)-xsng(jsta))        &
                +(ysng(ista)-ysng(jsta))*(ysng(ista)-ysng(jsta))
            IF(dist < rlimsq) THEN
              dzsq=(hgtsng(ista)-hgtsng(jsta)) *                        &
                   (hgtsng(ista)-hgtsng(jsta))
              IF(dzsq < zlimsq) THEN
                DO ivar=1,nvar
                  IF(qualsng(ivar,jsta) > 0) THEN
                    knt(ivar)=knt(ivar)+1
                    wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                    wgtsum(ivar)=wgtsum(ivar) + wgt
                    zsum(ivar)=zsum(ivar)+                              &
                          wgt*odifsng(ivar,jsta)
                  END IF
                END DO
              END IF
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
      IF(uasw > 0) THEN
        DO jsta=1,nobsua
          IF(iuseua(isrcua(jsta)) > 0) THEN
            dist=(xsng(ista)-xua(jsta))*(xsng(ista)-xua(jsta))          &
                +(ysng(ista)-yua(jsta))*(ysng(ista)-yua(jsta))
            IF(dist < rlimsq) THEN
              DO klev=1,nlevsua(jsta)
                dzsq=(hgtsng(ista)-hgtua(klev,jsta)) *                  &
                     (hgtsng(ista)-hgtua(klev,jsta))
                IF(dzsq < zlimsq) THEN
                  DO ivar=1,nvar
                    IF(qualua(ivar,klev,jsta) > 0) THEN
                      knt(ivar)=knt(ivar)+1
                      wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                      wgtsum(ivar)=wgtsum(ivar) + wgt
                      zsum(ivar)=zsum(ivar)+                            &
                          wgt*odifua(ivar,klev,jsta)
                    END IF
                  END DO
                END IF
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
      IF(radsw > 0) THEN
        DO jsta=1,ncolrad
          IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
            dist=(xsng(ista)-xradc(jsta))*(xsng(ista)-xradc(jsta))      &
                +(ysng(ista)-yradc(jsta))*(ysng(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO klev=1,nlevrad(jsta)
                dzsq=(hgtsng(ista)-hgtradc(klev,jsta)) *                &
                     (hgtsng(ista)-hgtradc(klev,jsta))
                IF(dzsq < zlimsq) THEN
                  IF(qualrad(1,klev,jsta) > 0) THEN
                    knt(1)=knt(1)+1
                    wgt=texp(-dist*rngsqi(1) -dzsq*zsqinv)
                    wgtsum(1)=wgtsum(1)+wgt*                            &
                              ABS(uazmrad(jsta))
                    zsum(1)=zsum(1)+                                    &
                         wgt*odifrad(1,klev,jsta) *                     &
                                ABS(uazmrad(jsta))
                  END IF
                  IF(qualrad(2,klev,jsta) > 0) THEN
                    knt(2)=knt(2)+1
                    wgt=texp(-dist*rngsqi(2) -dzsq*zsqinv)
                    wgtsum(2)=wgtsum(2)+wgt*                            &
                              ABS(vazmrad(jsta))
                    zsum(2)=zsum(2)+                                    &
                         wgt*odifrad(2,klev,jsta) *                     &
                                ABS(vazmrad(jsta))
                  END IF
                  IF(qualrad(3,klev,jsta) > 0) THEN
                    knt(5)=knt(5)+1
                    wgt=texp(-dist*rngsqi(5) -dzsq*zsqinv)
                    wgtsum(5)=wgtsum(5)+wgt
                    zsum(5)=zsum(5)+                                    &
                         wgt*odifrad(3,klev,jsta)
                  END IF
                END IF
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval obs
!
!-----------------------------------------------------------------------
!
      IF(retsw > 0) THEN
        DO jsta=1,ncolret
          IF(iuseret(isrcret(iret(jsta))) > 0) THEN
            dist=(xsng(ista)-xretc(jsta))*(xsng(ista)-xretc(jsta))      &
                +(ysng(ista)-yretc(jsta))*(ysng(ista)-yretc(jsta))
            IF(dist < rlimsq) THEN
              DO klev=1,nlevret(jsta)
                dzsq=(hgtsng(ista)-hgtretc(klev,jsta)) *                &
                     (hgtsng(ista)-hgtretc(klev,jsta))
                IF(dzsq < zlimsq) THEN
                  DO ivar=1,nvar
                    IF(qualret(ivar,klev,jsta) > 0) THEN
                      knt(ivar)=knt(ivar)+1
                      wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                      wgtsum(ivar)=wgtsum(ivar) + wgt
                      zsum(ivar)=zsum(ivar)+                            &
                          wgt*odifret(ivar,klev,jsta)
                    END IF
                  END DO
                END IF
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Divide by sum of weights
!
!-----------------------------------------------------------------------
!
      DO ivar=1,nvar
        IF(wgtsum(ivar) > wlim) oanxsng(ivar,ista)=oanxsng(ivar,ista) + &
                           zsum(ivar)/wgtsum(ivar)
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop over all multi-level obs
!
!-----------------------------------------------------------------------
!
  DO ista=1,nobsua
    IF(iuseua(isrcua(ista)) > 0) THEN
      DO k=1,nlevsua(ista)
        DO ivar=1,nvar
          zsum(ivar)=0.
          wgtsum(ivar)=0.
          knt(ivar)=0
        END DO
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
        IF(sngsw > 0) THEN
          DO jsta=1,nobsng
            IF(iusesng(isrcsng(jsta)) > 0) THEN
              dist=(xua(ista)-xsng(jsta))*(xua(ista)-xsng(jsta))        &
                  +(yua(ista)-ysng(jsta))*(yua(ista)-ysng(jsta))
              IF(dist < rlimsq) THEN
                dzsq=(hgtua(k,ista)-hgtsng(jsta)) *                     &
                     (hgtua(k,ista)-hgtsng(jsta))
                IF(dzsq < zlimsq) THEN
                  DO ivar=1,nvar
                    IF(qualsng(ivar,jsta) > 0) THEN
                      knt(ivar)=knt(ivar)+1
                      wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                      wgtsum(ivar)=wgtsum(ivar) + wgt
                      zsum(ivar)=zsum(ivar)+                            &
                           wgt*odifsng(ivar,jsta)
                    END IF
                  END DO
                END IF
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
        IF(uasw > 0) THEN
          DO jsta=1,nobsua
            IF(iuseua(isrcua(jsta)) > 0) THEN
              dist=(xua(ista)-xua(jsta))*(xua(ista)-xua(jsta))          &
                  +(yua(ista)-yua(jsta))*(yua(ista)-yua(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevsua(jsta)
                  dzsq=(hgtua(k,ista)-hgtua(klev,jsta)) *               &
                       (hgtua(k,ista)-hgtua(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualua(ivar,klev,jsta) > 0) THEN
                        knt(ivar)=knt(ivar)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar)=wgtsum(ivar) + wgt
                        zsum(ivar)=zsum(ivar)+                          &
                            wgt*odifua(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
        IF(radsw > 0) THEN
          DO jsta=1,ncolrad
            IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
              dist=(xua(ista)-xradc(jsta))*(xua(ista)-xradc(jsta))      &
                  +(yua(ista)-yradc(jsta))*(yua(ista)-yradc(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevrad(jsta)
                  dzsq=(hgtua(k,ista)-hgtradc(klev,jsta)) *             &
                       (hgtua(k,ista)-hgtradc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    IF(qualrad(1,klev,jsta) > 0) THEN
                      knt(1)=knt(1)+1
                      wgt=texp(-dist*rngsqi(1) -dzsq*zsqinv)
                      wgtsum(1)=wgtsum(1)+wgt*                          &
                               ABS(uazmrad(jsta))
                      zsum(1)=zsum(1)+                                  &
                           wgt*odifrad(1,klev,jsta) *                   &
                               ABS(uazmrad(jsta))
                    END IF
                    IF(qualrad(2,klev,jsta) > 0) THEN
                      knt(2)=knt(2)+1
                      wgt=texp(-dist*rngsqi(2) -dzsq*zsqinv)
                      wgtsum(2)=wgtsum(2)+wgt*                          &
                               ABS(vazmrad(jsta))
                      zsum(2)=zsum(2)+                                  &
                           wgt*odifrad(2,klev,jsta) *                   &
                                  ABS(vazmrad(jsta))
                    END IF
                    IF(qualrad(3,klev,jsta) > 0) THEN
                      knt(5)=knt(5)+1
                      wgt=texp(-dist*rngsqi(5) -dzsq*zsqinv)
                      wgtsum(5)=wgtsum(5)+wgt
                      zsum(5)=zsum(5)+                                  &
                           wgt*odifrad(3,klev,jsta)
                    END IF
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval obs
!
!-----------------------------------------------------------------------
!
        IF(retsw > 0) THEN
          DO jsta=1,ncolret
            IF(iuseret(isrcret(iret(jsta))) > 0) THEN
              dist=(xua(ista)-xretc(jsta))*(xua(ista)-xretc(jsta))      &
                  +(yua(ista)-yretc(jsta))*(yua(ista)-yretc(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevret(jsta)
                  dzsq=(hgtua(k,ista)-hgtretc(klev,jsta)) *             &
                       (hgtua(k,ista)-hgtretc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualret(ivar,klev,jsta) > 0) THEN
                        knt(ivar)=knt(ivar)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar)=wgtsum(ivar) + wgt
                        zsum(ivar)=zsum(ivar)+                          &
                            wgt*odifret(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Divide by sum of weights
!
!-----------------------------------------------------------------------
!
        DO ivar=1,nvar
          IF(wgtsum(ivar) > wlim) oanxua(ivar,k,ista)=oanxua(ivar,k,ista) + &
              zsum(ivar)/wgtsum(ivar)
        END DO
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop over all radar columns
!
!-----------------------------------------------------------------------
!
  DO ista=1,ncolrad
    IF(iuserad(isrcrad(irad(ista))) > 0) THEN
      DO k=1,nlevrad(ista)
        IF(qualrad(1,k,ista) > 0 .AND. qualrad(2,k,ista) > 0 ) THEN
          DO ivar=1,nvar
            zsum(ivar)=0.
            wgtsum(ivar)=0.
            knt(ivar)=0
          END DO
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
          IF(sngsw > 0) THEN
            DO jsta=1,nobsng
              IF(iusesng(isrcsng(jsta)) > 0) THEN
                dist=(xradc(ista)-xsng(jsta))*(xradc(ista)-xsng(jsta))  &
                    +(yradc(ista)-ysng(jsta))*(yradc(ista)-ysng(jsta))
                IF(dist < rlimsq) THEN
                  dzsq=(hgtradc(k,ista)-hgtsng(jsta)) *                 &
                       (hgtradc(k,ista)-hgtsng(jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualsng(ivar,jsta) > 0) THEN
                        knt(ivar)=knt(ivar)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar)=wgtsum(ivar) + wgt
                        zsum(ivar)=zsum(ivar)+                          &
                             wgt*odifsng(ivar,jsta)
                      END IF
                    END DO
                  END IF
                END IF     ! horizontal dist check
              END IF     ! source check
            END DO
          END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
          IF(uasw > 0) THEN
            DO jsta=1,nobsua
              IF(iuseua(isrcua(jsta)) > 0) THEN
                dist=(xradc(ista)-xua(jsta))*(xradc(ista)-xua(jsta))    &
                    +(yradc(ista)-yua(jsta))*(yradc(ista)-yua(jsta))
                IF(dist < rlimsq) THEN
                  DO klev=1,nlevsua(jsta)
                    dzsq=(hgtradc(k,ista)-hgtua(klev,jsta)) *           &
                         (hgtradc(k,ista)-hgtua(klev,jsta))
                    IF(dzsq < zlimsq) THEN
                      DO ivar=1,nvar
                        IF(qualua(ivar,klev,jsta) > 0) THEN
                          knt(ivar)=knt(ivar)+1
                          wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                          wgtsum(ivar)=wgtsum(ivar) + wgt
                          zsum(ivar)=zsum(ivar)+                        &
                              wgt*odifua(ivar,klev,jsta)
                        END IF
                      END DO
                    END IF
                  END DO
                END IF     ! horizontal dist check
              END IF     ! source check
            END DO
          END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
          IF(radsw > 0) THEN
            DO jsta=1,ncolrad
              IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
                dist=(xradc(ista)-xradc(jsta))*(xradc(ista)-xradc(jsta)) &
                    +(yradc(ista)-yradc(jsta))*(yradc(ista)-yradc(jsta))
                IF(dist < rlimsq) THEN
                  DO klev=1,nlevrad(jsta)
                    dzsq=(hgtradc(k,ista)-hgtradc(klev,jsta)) *         &
                         (hgtradc(k,ista)-hgtradc(klev,jsta))
                    IF(dzsq < zlimsq) THEN
                      IF(qualrad(1,klev,jsta) > 0) THEN
                        knt(1)=knt(1)+1
                        wgt=texp(-dist*rngsqi(1) -dzsq*zsqinv)
                        wgtsum(1)=wgtsum(1)+wgt*                        &
                               ABS(uazmrad(jsta))
                        zsum(1)=zsum(1)+                                &
                            wgt*odifrad(1,klev,jsta) *                  &
                                   ABS(uazmrad(jsta))
                      END IF
                      IF(qualrad(2,klev,jsta) > 0) THEN
                        knt(2)=knt(2)+1
                        wgt=texp(-dist*rngsqi(2) -dzsq*zsqinv)
                        wgtsum(2)=wgtsum(2)+wgt*                        &
                               ABS(vazmrad(jsta))
                        zsum(2)=zsum(2)+                                &
                            wgt*odifrad(2,klev,jsta) *                  &
                                  ABS(vazmrad(jsta))
                      END IF
                      IF(qualrad(3,klev,jsta) > 0) THEN
                        knt(5)=knt(5)+1
                        wgt=texp(-dist*rngsqi(5) -dzsq*zsqinv)
                        wgtsum(5)=wgtsum(5)+wgt
                        zsum(5)=zsum(5)+                                &
                            wgt*odifrad(3,klev,jsta)
                      END IF
                    END IF
                  END DO
                END IF     ! horizontal dist check
              END IF     ! source check
            END DO
          END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar retrieval obs
!
!-----------------------------------------------------------------------
!
          IF(retsw > 0) THEN
            DO jsta=1,ncolret
              IF(iuseret(isrcret(iret(jsta))) > 0) THEN
                dist=(xradc(ista)-xretc(jsta))*(xradc(ista)-xretc(jsta)) &
                    +(yradc(ista)-yretc(jsta))*(yradc(ista)-yretc(jsta))
                IF(dist < rlimsq) THEN
                  DO klev=1,nlevret(jsta)
                    dzsq=(hgtradc(k,ista)-hgtretc(klev,jsta)) *         &
                         (hgtradc(k,ista)-hgtretc(klev,jsta))
                    IF(dzsq < zlimsq) THEN
                      DO ivar=1,nvar
                        IF(qualret(ivar,klev,jsta) > 0) THEN
                          knt(ivar)=knt(ivar)+1
                          wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                          wgtsum(ivar)=wgtsum(ivar) + wgt
                          zsum(ivar)=zsum(ivar)+                        &
                               wgt*odifret(ivar,klev,jsta)
                        END IF
                      END DO
                    END IF
                  END DO
                END IF     ! horizontal dist check
              END IF     ! source check
            END DO
          END IF
!
!-----------------------------------------------------------------------
!
!  Divide by sum of weights
!
!-----------------------------------------------------------------------
!
          IF(wgtsum(1) > wlim) oanxrad(1,k,ista)=oanxrad(1,k,ista) +    &
                zsum(1)/wgtsum(1)
          IF(wgtsum(2) > wlim) oanxrad(2,k,ista)=oanxrad(2,k,ista) +    &
                zsum(2)/wgtsum(2)
          IF(wgtsum(5) > wlim) oanxrad(3,k,ista)=oanxrad(3,k,ista) +    &
                zsum(5)/wgtsum(5)
        END IF
      END DO
    END IF
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop over all radar retrievals
!
!-----------------------------------------------------------------------
!
  DO ista=1,ncolret
    IF(iuseret(isrcret(iret(ista))) > 0) THEN
      DO k=1,nlevret(ista)
        DO ivar=1,nvar
          zsum(ivar)=0.
          wgtsum(ivar)=0.
          knt(ivar)=0
        END DO
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
        IF(sngsw > 0) THEN
          DO jsta=1,nobsng
            IF(iusesng(isrcsng(jsta)) > 0) THEN
              dist=(xretc(ista)-xsng(jsta))*(xretc(ista)-xsng(jsta))    &
                  +(yretc(ista)-ysng(jsta))*(yretc(ista)-ysng(jsta))
              IF(dist < rlimsq) THEN
                dzsq=(hgtretc(k,ista)-hgtsng(jsta)) *                   &
                     (hgtretc(k,ista)-hgtsng(jsta))
                IF(dzsq < zlimsq) THEN
                  DO ivar=1,nvar
                    IF(qualsng(ivar,jsta) > 0) THEN
                      knt(ivar)=knt(ivar)+1
                      wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                      wgtsum(ivar)=wgtsum(ivar) + wgt
                      zsum(ivar)=zsum(ivar)+                            &
                           wgt*odifsng(ivar,jsta)
                    END IF
                  END DO
                END IF
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
        IF(uasw > 0) THEN
          DO jsta=1,nobsua
            IF(iuseua(isrcua(jsta)) > 0) THEN
              dist=(xretc(ista)-xua(jsta))*(xretc(ista)-xua(jsta))      &
                  +(yretc(ista)-yua(jsta))*(yretc(ista)-yua(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevsua(jsta)
                  dzsq=(hgtretc(k,ista)-hgtua(klev,jsta)) *             &
                       (hgtretc(k,ista)-hgtua(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualua(ivar,klev,jsta) > 0) THEN
                        knt(ivar)=knt(ivar)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar)=wgtsum(ivar) + wgt
                        zsum(ivar)=zsum(ivar)+                          &
                            wgt*odifua(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
        IF(radsw > 0) THEN
          DO jsta=1,ncolrad
            IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
              dist=(xretc(ista)-xradc(jsta))*(xretc(ista)-xradc(jsta))  &
                  +(yretc(ista)-yradc(jsta))*(yretc(ista)-yradc(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevrad(jsta)
                  dzsq=(hgtretc(k,ista)-hgtradc(klev,jsta)) *           &
                       (hgtretc(k,ista)-hgtradc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    IF(qualrad(1,klev,jsta) > 0) THEN
                      knt(1)=knt(1)+1
                      wgt=texp(-dist*rngsqi(1) -dzsq*zsqinv)
                      wgtsum(1)=wgtsum(1)+wgt*                          &
                               ABS(uazmrad(jsta))
                      zsum(1)=zsum(1)+                                  &
                           wgt*odifrad(1,klev,jsta) *                   &
                               ABS(uazmrad(jsta))
                    END IF
                    IF(qualrad(2,klev,jsta) > 0) THEN
                      knt(2)=knt(2)+1
                      wgt=texp(-dist*rngsqi(2) -dzsq*zsqinv)
                      wgtsum(2)=wgtsum(2)+wgt*                          &
                               ABS(vazmrad(jsta))
                      zsum(2)=zsum(2)+                                  &
                           wgt*odifrad(2,klev,jsta) *                   &
                                  ABS(vazmrad(jsta))
                    END IF
                    IF(qualrad(3,klev,jsta) > 0) THEN
                      knt(5)=knt(5)+1
                      wgt=texp(-dist*rngsqi(5) -dzsq*zsqinv)
                      wgtsum(5)=wgtsum(5)+wgt
                      zsum(5)=zsum(5)+                                  &
                           wgt*odifrad(3,klev,jsta)
                    END IF
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval obs
!
!-----------------------------------------------------------------------
!
        IF(retsw > 0) THEN
          DO jsta=1,ncolret
            IF(iuseret(isrcret(jsta)) > 0) THEN
              dist=(xretc(ista)-xretc(jsta))*(xretc(ista)-xretc(jsta))  &
                  +(yretc(ista)-yretc(jsta))*(yretc(ista)-yretc(jsta))
              IF(dist < rlimsq) THEN
                DO klev=1,nlevret(jsta)
                  dzsq=(hgtretc(k,ista)-hgtretc(klev,jsta)) *           &
                       (hgtretc(k,ista)-hgtretc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualret(ivar,klev,jsta) > 0) THEN
                        knt(ivar)=knt(ivar)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar)=wgtsum(ivar) + wgt
                        zsum(ivar)=zsum(ivar)+                          &
                            wgt*odifret(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END IF     ! horizontal dist check
            END IF     ! source check
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Divide by sum of weights
!
!-----------------------------------------------------------------------
!
        DO ivar=1,nvar
          IF(wgtsum(ivar) > wlim) oanxret(ivar,k,ista)=oanxret(ivar,k,ista) + &
              zsum(ivar)/wgtsum(ivar)
        END DO
      END DO
    END IF
  END DO
!
  istatus=0
  RETURN
END SUBROUTINE barobs3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BARGRD3d                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE bargrd3d(nx,ny,nz,                                           & 1
           nvar,nvarrad,nzua,nzrdr,nzret,                               &
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,xs,ys,zs,                     &
           odifsng,xsng,ysng,hgtsng,qualsng,isrcsng,nobsng,             &
           odifua,xua,yua,hgtua,qualua,isrcua,nlevsua,nobsua,           &
           odifrad,xradc,yradc,hgtradc,                                 &
           uazmrad,vazmrad,qualrad,irad,isrcrad,nlevrad,ncolrad,        &
           odifret,xretc,yretc,hgtretc,                                 &
           qualret,iret,isrcret,nlevret,ncolret,                        &
           iusesng,iuseua,iuserad,iuseret,                              &
           sngsw,uasw,radsw,retsw,kpvrsq,rpass,                         &
           wlim,zrngsq,zwlim,                                           &
           rngsqi,knt,wgtsum,zsum,                                      &
           anx,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Objectively analyze data on a regular grid specified
!  by x and y.  Barnes weight function is employed.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  September, 1995
!
!  MODIFICATION HISTORY:
!
!  July, 1996 (KB)
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!
!-----------------------------------------------------------------------
!
!  Sizing Arguments
!
  INTEGER :: nx,ny,nz
  INTEGER :: nvar,nvarrad,nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret
!
!-----------------------------------------------------------------------
!
!  Grid arguments
!
!-----------------------------------------------------------------------
!
  REAL :: xs(nx)
  REAL :: ys(ny)
  REAL :: zs(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: odifsng(nvar,mxsng),                                          &
       xsng(mxsng),                                                     &
       ysng(mxsng),                                                     &
       hgtsng(mxsng)
  INTEGER :: qualsng(nvar,mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: nobsng
  REAL :: odifua(nvar,nzua,mxua)             ! variable to analyse
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  INTEGER :: qualua(nvar,nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  INTEGER :: qualret(nvar,nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Scratch space
!
!-----------------------------------------------------------------------
!
  REAL :: rngsqi(nvar)
  INTEGER :: knt(nvar,nz)
  REAL :: wgtsum(nvar,nz)
  REAL :: zsum(nvar,nz)
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: iusesng(0:nsrcsng)
  INTEGER :: iuseua(0:nsrcua)
  INTEGER :: iuserad(0:nsrcrad)
  INTEGER :: iuseret(0:nsrcret)
  INTEGER :: sngsw,uasw,radsw,retsw
  REAL :: kpvrsq(nvar)
  REAL :: rpass,                                                        &
       wlim
  REAL :: zrngsq,                                                       &
       zwlim
!
!-----------------------------------------------------------------------
!
!  Output grid arguments
!
!-----------------------------------------------------------------------
!
  REAL :: anx(nx,ny,nz,nvar)
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  REAL :: texp
  REAL :: rlimsq,zlimsq,zsqinv
  REAL :: dist,dzsq,wgt
  INTEGER :: i,j,k,jsta,klev,ivar
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rlimsq=-rpass*ALOG(wlim)
  DO ivar=1,nvar
    rngsqi(ivar)=1./(kpvrsq(ivar)*rpass)
  END DO
!
  zlimsq=-zrngsq*ALOG(zwlim)
  zsqinv=1./zrngsq
!
!-----------------------------------------------------------------------
!
!  Loop over all grid points
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      DO k=1,nz-1
        DO ivar=1,nvar
          knt(ivar,k)=0
          wgtsum(ivar,k)=0.
          zsum(ivar,k)=0.
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
      IF(sngsw > 0) THEN
        DO jsta=1,nobsng
          IF(iusesng(isrcsng(jsta)) > 0) THEN
            dist=(xs(i)-xsng(jsta))*(xs(i)-xsng(jsta))                  &
                +(ys(j)-ysng(jsta))*(ys(j)-ysng(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nz-1
                dzsq=(zs(i,j,k)-hgtsng(jsta)) *                         &
                     (zs(i,j,k)-hgtsng(jsta))
                IF(dzsq < zlimsq) THEN
                  DO ivar=1,nvar
                    IF(qualsng(ivar,jsta) > 0) THEN
                      knt(ivar,k)=knt(ivar,k)+1
                      wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                      wgtsum(ivar,k)=wgtsum(ivar,k) + wgt
                      zsum(ivar,k)=zsum(ivar,k)+                        &
                            wgt*odifsng(ivar,jsta)
                    END IF
                  END DO
                END IF
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
      IF(uasw > 0) THEN
        DO jsta=1,nobsua
          IF(iuseua(isrcua(jsta)) > 0) THEN
            dist=(xs(i)-xua(jsta))*(xs(i)-xua(jsta))                    &
                +(ys(j)-yua(jsta))*(ys(j)-yua(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nz-1
                DO klev=1,nlevsua(jsta)
                  dzsq=(zs(i,j,k)-hgtua(klev,jsta)) *                   &
                       (zs(i,j,k)-hgtua(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualua(ivar,klev,jsta) > 0) THEN
                        knt(ivar,k)=knt(ivar,k)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar,k)=wgtsum(ivar,k) + wgt
                        zsum(ivar,k)=zsum(ivar,k)+                      &
                               wgt*odifua(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
      IF(radsw > 0) THEN
        DO jsta=1,ncolrad
          IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
            dist=(xs(i)-xradc(jsta))*(xs(i)-xradc(jsta))                &
                +(ys(j)-yradc(jsta))*(ys(j)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nz-1
                DO klev=1,nlevrad(jsta)
                  dzsq=(zs(i,j,k)-hgtradc(klev,jsta)) *                 &
                       (zs(i,j,k)-hgtradc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    IF(qualrad(1,klev,jsta) > 0) THEN
                      knt(1,k)=knt(1,k)+1
                      wgt=texp(-dist*rngsqi(1) -dzsq*zsqinv)
                      wgtsum(1,k)=wgtsum(1,k)+                          &
                                  wgt*uazmrad(jsta)
                      zsum(1,k)=zsum(1,k)+                              &
                          wgt*odifrad(1,klev,jsta) *                    &
                                   uazmrad(jsta)
                    END IF
                    IF(qualrad(2,klev,jsta) > 0) THEN
                      knt(2,k)=knt(2,k)+1
                      wgt=texp(-dist*rngsqi(2) -dzsq*zsqinv)
                      wgtsum(2,k)=wgtsum(2,k)+                          &
                                  wgt*vazmrad(jsta)
                      zsum(2,k)=zsum(2,k)+                              &
                          wgt*odifrad(2,klev,jsta) *                    &
                                   vazmrad(jsta)
                    END IF
                    IF(qualrad(3,klev,jsta) > 0) THEN
                      knt(5,k)=knt(5,k)+1
                      wgt=texp(-dist*rngsqi(5) -dzsq*zsqinv)
                      wgtsum(5,k)=wgtsum(5,k)+wgt
                      zsum(5,k)=zsum(5,k)+                              &
                          wgt*odifrad(3,klev,jsta)
                    END IF
                  END IF
                END DO
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval obs
!
!-----------------------------------------------------------------------
!
      IF(retsw > 0) THEN
        DO jsta=1,ncolret
          IF(iuseret(isrcret(iret(jsta))) > 0) THEN
            dist=(xs(i)-xretc(jsta))*(xs(i)-xretc(jsta))                &
                +(ys(j)-yretc(jsta))*(ys(j)-yretc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nz-1
                DO klev=1,nlevret(jsta)
                  dzsq=(zs(i,j,k)-hgtretc(klev,jsta)) *                 &
                       (zs(i,j,k)-hgtretc(klev,jsta))
                  IF(dzsq < zlimsq) THEN
                    DO ivar=1,nvar
                      IF(qualret(ivar,klev,jsta) > 0) THEN
                        knt(ivar,k)=knt(ivar,k)+1
                        wgt=texp(-dist*rngsqi(ivar) -dzsq*zsqinv)
                        wgtsum(ivar,k)=wgtsum(ivar,k) + wgt
                        zsum(ivar,k)=zsum(ivar,k)+                      &
                               wgt*odifret(ivar,klev,jsta)
                      END IF
                    END DO
                  END IF
                END DO
              END DO
            END IF     ! horizontal dist check
          END IF     ! source check
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Divide by sum of weights
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO ivar=1,nvar
          IF(wgtsum(ivar,k) > wlim) anx(i,j,k,ivar)=anx(i,j,k,ivar) +   &
                            zsum(ivar,k)/wgtsum(ivar,k)
        END DO
      END DO
    END DO
  END DO
!
  istatus=0
  RETURN
END SUBROUTINE bargrd3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE OBSCOR                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE obscor(nvar,nvarrad,nzua,nzrdr,nzret,                        & 2
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                    &
           qobsng,xsng,ysng,hgtsng,                                     &
           qualsng,isrcsng,icatsng,nobsng,                              &
           qobsua,xua,yua,hgtua,                                        &
           qualua,isrcua,nlevsua,nobsua,                                &
           qobsrad,xradc,yradc,hgtradc,                                 &
           uazmrad,vazmrad,                                             &
           qualrad,irad,isrcrad,nlevrad,ncolrad,                        &
           qobsret,xretc,yretc,hgtretc,                                 &
           qualret,iret,isrcret,nlevret,ncolret,                        &
           iusesng,iuseua,iuserad,iuseret,                              &
           sngsw,uasw,radsw,retsw,kpvrsq,rpass,rngsqi,                  &
           wlim,zrngsq,                                                 &
           corsng,corua,corrad,corret,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate covariance sums at observation points for
!  Bratseth analysis method.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  July, 1995
!
!  MODIFICATION HISTORY:
!
!  Jan, 1996 (KB)
!  Added radar data and other improvements.
!
!  July, 1996 (KB)
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!
!  Jan 5, 1997 (KB)
!  New version to reorder summing to speed-up code.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Sizing arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nvar,nvarrad,nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
  REAL :: xcor(ncat,ncat)
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: qobsng(nvar,mxsng),                                           &
       xsng(mxsng),                                                     &
       ysng(mxsng),                                                     &
       hgtsng(mxsng)
  INTEGER :: qualsng(nvar,mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: icatsng(mxsng)
  INTEGER :: nobsng
!
  REAL :: qobsua(nvar,nzua,mxua)
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  INTEGER :: qualua(nvar,nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: qobsrad(nvarrad,nzrdr,mxcolrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  INTEGER :: qualrad(nvarrad,nzrdr,mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: qobsret(nvar,nzret,mxcolret)
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  INTEGER :: qualret(nvar,nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: iusesng(0:nsrcsng)
  INTEGER :: iuseua(0:nsrcua)
  INTEGER :: iuserad(0:nsrcrad)
  INTEGER :: iuseret(0:nsrcret)
  INTEGER :: sngsw,uasw,radsw,retsw
  REAL :: kpvrsq(nvar)
  REAL :: rpass,                                                        &
       wlim
  REAL :: zrngsq
  REAL :: rngsqi(nvar)
!
!-----------------------------------------------------------------------
!
!  Output arguments
!
!-----------------------------------------------------------------------
!
  REAL :: corsng(nvar,mxsng)
  REAL :: corua(nvar,nzua,mxua)
  REAL :: corrad(nvarrad,nzrdr,mxcolrad)
  REAL :: corret(nvar,nzret,mxcolret)
  INTEGER :: istatus
!
  REAL :: ftabinv
  INTEGER :: ntabexp
  PARAMETER (ntabexp=1000)
  REAL :: tabexp(ntabexp)
  COMMON /ftabexp/ ftabinv
  COMMON /tablexp/ tabexp
!
!-----------------------------------------------------------------------
!
!  Misc.local variables
!
!-----------------------------------------------------------------------
!
  REAL :: rsq,rlim,rlimsq,zsqinv
  REAL :: dist,dzsq,rnorm
  INTEGER :: ista,jsta,k,klev,ivar,INDEX
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rlimsq=0.
  DO ivar=1,nvar
    rsq=(kpvrsq(ivar)*rpass)
    rngsqi(ivar)=1./rsq
    rlimsq=AMAX1(rlimsq,rsq)
  END DO
  rlimsq=-rlimsq*ALOG(wlim)
  rlim=0.001*SQRT(rlimsq)
  WRITE(6,'(/a,f10.2,a)')                                               &
       ' OBSCOR: Influence cutoff radius: ',rlim,' km.'

  zsqinv=1./zrngsq
!
!-----------------------------------------------------------------------
!
!  Initialize covariance sums with observation error variance
!
!-----------------------------------------------------------------------
!
! OpenMP:
!$OMP PARALLEL DO PRIVATE(ista,ivar)
  DO ista=1,nobsng
    DO ivar=1,nvar
      corsng(ivar,ista)=qobsng(ivar,ista)
    END DO
  END DO
!$OMP PARALLEL DO PRIVATE(ista,k,ivar)
  DO ista=1,nobsua
    DO k=1,nlevsua(ista)
      DO ivar=1,nvar
        corua(ivar,k,ista)=qobsua(ivar,k,ista)
      END DO
    END DO
  END DO
!$OMP PARALLEL DO PRIVATE(ista,k,ivar)
  DO ista=1,ncolret
    DO k=1,nlevret(ista)
      DO ivar=1,nvar
        corret(ivar,k,ista)=qobsret(ivar,k,ista)
      END DO
    END DO
  END DO
!$OMP PARALLEL DO PRIVATE(ista,k,ivar)
  DO ista=1,ncolrad
    DO k=1,nlevrad(ista)
      DO ivar=1,nvarrad
        corrad(ivar,k,ista)=qobsrad(ivar,k,ista)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
  IF(sngsw > 0) THEN
    DO jsta=1,nobsng
      IF(iusesng(isrcsng(jsta)) > 0) THEN
        DO ivar=1,nvar
          IF(qualsng(ivar,jsta) > 0) THEN
!
!-----------------------------------------------------------------------
!
!  Add to single-level covariance sum
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsng
              dist=(xsng(ista)-xsng(jsta))*(xsng(ista)-xsng(jsta))      &
                  +(ysng(ista)-ysng(jsta))*(ysng(ista)-ysng(jsta))
              IF(dist < rlimsq) THEN
                dzsq=(hgtsng(ista)-hgtsng(jsta)) *                      &
                     (hgtsng(ista)-hgtsng(jsta))
                rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corsng(ivar,ista)=corsng(ivar,ista)+ &
                    tabexp(INDEX)*xcor(icatsng(ista),icatsng(jsta))
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level covariance sum
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsua
              dist=(xua(ista)-xsng(jsta))*(xua(ista)-xsng(jsta))        &
                  +(yua(ista)-ysng(jsta))*(yua(ista)-ysng(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevsua(ista)
                  dzsq=(hgtua(k,ista)-hgtsng(jsta)) *                   &
                       (hgtua(k,ista)-hgtsng(jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corua(ivar,k,ista)=               &
                      corua(ivar,k,ista)+tabexp(INDEX)
                END DO
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval covariance sum
!
!-----------------------------------------------------------------------
!
            DO ista=1,ncolret
              dist=(xretc(ista)-xsng(jsta))*(xretc(ista)-xsng(jsta))    &
                  +(yretc(ista)-ysng(jsta))*(yretc(ista)-ysng(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevret(ista)
                  dzsq=(hgtretc(k,ista)-hgtsng(jsta)) *                 &
                       (hgtretc(k,ista)-hgtsng(jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corret(ivar,k,ista)=              &
                      corret(ivar,k,ista)+tabexp(INDEX)
                END DO
              END IF
            END DO
          END IF
        END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar data covariance sum
!
!-----------------------------------------------------------------------
!
        DO ista=1,ncolrad
          dist=(xradc(ista)-xsng(jsta))*(xradc(ista)-xsng(jsta))        &
              +(yradc(ista)-ysng(jsta))*(yradc(ista)-ysng(jsta))
          IF(dist < rlimsq) THEN
            DO k=1,nlevrad(ista)
              dzsq=(hgtradc(k,ista)-hgtsng(jsta)) *                     &
                   (hgtradc(k,ista)-hgtsng(jsta))
!
              IF(qualsng(1,jsta) > 0) THEN
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corrad(1,k,ista)=                   &
                    corrad(1,k,ista)+tabexp(INDEX)
              END IF
!
              IF(qualsng(2,jsta) > 0) THEN
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corrad(2,k,ista)=                   &
                    corrad(2,k,ista)+tabexp(INDEX)
              END IF
!
              IF(qualsng(5,jsta) > 0) THEN
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corrad(3,k,ista)=                   &
                    corrad(3,k,ista)+tabexp(INDEX)
              END IF
            END DO
          END IF
        END DO
      END IF     ! source check
    END DO
  END IF   ! sngsw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
  IF(uasw > 0) THEN
    DO jsta=1,nobsua
      IF(iuseua(isrcua(jsta)) > 0) THEN
        DO klev=1,nlevsua(jsta)
          DO ivar=1,nvar
            IF(qualua(ivar,klev,jsta) > 0) THEN
!
!-----------------------------------------------------------------------
!
!  Add to single-level covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,nobsng
                dist=(xsng(ista)-xua(jsta))*(xsng(ista)-xua(jsta))      &
                    +(ysng(ista)-yua(jsta))*(ysng(ista)-yua(jsta))
                IF(dist < rlimsq) THEN
                  dzsq=(hgtsng(ista)-hgtua(klev,jsta)) *                &
                       (hgtsng(ista)-hgtua(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corsng(ivar,ista)=                &
                      corsng(ivar,ista)+tabexp(INDEX)
                END IF
              END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,nobsua
                dist=(xua(ista)-xua(jsta))*(xua(ista)-xua(jsta))        &
                    +(yua(ista)-yua(jsta))*(yua(ista)-yua(jsta))
                IF(dist < rlimsq) THEN
                  DO k=1,nlevsua(ista)
                    dzsq=(hgtua(k,ista)-hgtua(klev,jsta)) *             &
                         (hgtua(k,ista)-hgtua(klev,jsta))
                    rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                    INDEX=nint(rnorm+1)
                    IF(INDEX < ntabexp) corua(ivar,k,ista)=             &
                        corua(ivar,k,ista)+tabexp(INDEX)
                  END DO
                END IF
              END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,ncolret
                dist=(xretc(ista)-xua(jsta))*(xretc(ista)-xua(jsta))    &
                    +(yretc(ista)-yua(jsta))*(yretc(ista)-yua(jsta))
                IF(dist < rlimsq) THEN
                  DO k=1,nlevret(ista)
                    dzsq=(hgtretc(k,ista)-hgtua(klev,jsta)) *           &
                         (hgtretc(k,ista)-hgtua(klev,jsta))
                    rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                    INDEX=nint(rnorm+1)
                    IF(INDEX < ntabexp) corret(ivar,k,ista)=            &
                        corret(ivar,k,ista)+tabexp(INDEX)
                  END DO
                END IF
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xua(jsta))*(xradc(ista)-xua(jsta))        &
                +(yradc(ista)-yua(jsta))*(yradc(ista)-yua(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtua(klev,jsta)) *               &
                     (hgtradc(k,ista)-hgtua(klev,jsta))
!
                IF(qualua(1,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(1,k,ista)=                 &
                      corrad(1,k,ista)+tabexp(INDEX)
                END IF
!
                IF(qualua(2,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(2,k,ista)=                 &
                      corrad(2,k,ista)+tabexp(INDEX)
                END IF
!
                IF(qualua(5,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(3,k,ista)=                 &
                      corrad(3,k,ista)+tabexp(INDEX)
                END IF
              END DO
            END IF
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF   ! uasw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval data
!
!-----------------------------------------------------------------------
!
  IF(retsw > 0) THEN
    DO jsta=1,ncolret
      IF(iuseret(isrcret(iret(jsta))) > 0) THEN
        DO klev=1,nlevret(jsta)
          DO ivar=1,nvar
            IF(qualret(ivar,klev,jsta) > 0) THEN
!
!-----------------------------------------------------------------------
!
!  Add to single-level covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,nobsng
                dist=(xsng(ista)-xretc(jsta))*(xsng(ista)-xretc(jsta))  &
                    +(ysng(ista)-yretc(jsta))*(ysng(ista)-yretc(jsta))
                IF(dist < rlimsq) THEN
                  dzsq=(hgtsng(ista)-hgtretc(klev,jsta)) *              &
                       (hgtsng(ista)-hgtretc(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corsng(ivar,ista)=                &
                      corsng(ivar,ista)+tabexp(INDEX)
                END IF
              END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,nobsua
                dist=(xua(ista)-xretc(jsta))*(xua(ista)-xretc(jsta))    &
                    +(yua(ista)-yretc(jsta))*(yua(ista)-yretc(jsta))
                IF(dist < rlimsq) THEN
                  DO k=1,nlevsua(ista)
                    dzsq=(hgtua(k,ista)-hgtretc(klev,jsta)) *           &
                         (hgtua(k,ista)-hgtretc(klev,jsta))
                    rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                    INDEX=nint(rnorm+1)
                    IF(INDEX < ntabexp) corua(ivar,k,ista)=             &
                        corua(ivar,k,ista)+tabexp(INDEX)
                  END DO
                END IF
              END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval covariance sum
!
!-----------------------------------------------------------------------
!
              DO ista=1,ncolret
                dist=(xretc(ista)-xretc(jsta))*                         &
                     (xretc(ista)-xretc(jsta))                          &
                    +(yretc(ista)-yretc(jsta))*                         &
                     (yretc(ista)-yretc(jsta))
                IF(dist < rlimsq) THEN
                  DO k=1,nlevret(ista)
                    dzsq=(hgtretc(k,ista)-hgtretc(klev,jsta)) *         &
                         (hgtretc(k,ista)-hgtretc(klev,jsta))
                    rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                    INDEX=nint(rnorm+1)
                    IF(INDEX < ntabexp) corret(ivar,k,ista)=            &
                        corret(ivar,k,ista)+tabexp(INDEX)
                  END DO
                END IF
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xretc(jsta))*(xradc(ista)-xretc(jsta))    &
                +(yradc(ista)-yretc(jsta))*(yradc(ista)-yretc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtretc(klev,jsta)) *             &
                     (hgtradc(k,ista)-hgtretc(klev,jsta))
!
                IF(qualret(1,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(1,k,ista)=                 &
                      corrad(1,k,ista)+tabexp(INDEX)
                END IF
!
                IF(qualret(2,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(2,k,ista)=                 &
                      corrad(2,k,ista)+tabexp(INDEX)
                END IF
!
                IF(qualret(5,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(3,k,ista)=                 &
                      corrad(3,k,ista)+tabexp(INDEX)
                END IF
              END DO
            END IF
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF   ! retsw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
  IF(radsw > 0) THEN
    DO jsta=1,ncolrad
      IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
        DO klev=1,nlevrad(jsta)
!
!-----------------------------------------------------------------------
!
!  Add to single-level covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsng
            dist=(xsng(ista)-xradc(jsta))*(xsng(ista)-xradc(jsta))      &
                +(ysng(ista)-yradc(jsta))*(ysng(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              dzsq=(hgtsng(ista)-hgtradc(klev,jsta)) *                  &
                   (hgtsng(ista)-hgtradc(klev,jsta))
!
              IF(qualrad(1,klev,jsta) > 0) THEN
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corsng(1,ista)=corsng(1,ista)+      &
                             (tabexp(INDEX)*ABS(uazmrad(jsta)))
              END IF
!
              IF(qualrad(2,klev,jsta) > 0) THEN
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corsng(2,ista)=corsng(2,ista)+      &
                             (tabexp(INDEX)*ABS(vazmrad(jsta)))
              END IF
!
              IF(qualrad(3,klev,jsta) > 0) THEN
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) corsng(5,ista)=corsng(5,ista)+tabexp(INDEX)
              END IF
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsua
            dist=(xua(ista)-xradc(jsta))*(xua(ista)-xradc(jsta))        &
                +(yua(ista)-yradc(jsta))*(yua(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevsua(ista)
                dzsq=(hgtua(k,ista)-hgtradc(klev,jsta)) *               &
                     (hgtua(k,ista)-hgtradc(klev,jsta))
!
                IF(qualrad(1,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corua(1,k,ista)=corua(1,k,ista)+  &
                             (tabexp(INDEX)*ABS(uazmrad(jsta)))
                END IF
!
                IF(qualrad(2,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corua(2,k,ista)=corua(2,k,ista)+  &
                             (tabexp(INDEX)*ABS(vazmrad(jsta)))
                END IF
!
                IF(qualrad(3,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp)                                   &
                      corua(5,k,ista)=corua(5,k,ista)+tabexp(INDEX)
                END IF
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolret
            dist=(xretc(ista)-xradc(jsta))*(xretc(ista)-xradc(jsta))    &
                +(yretc(ista)-yradc(jsta))*(yretc(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevret(ista)
                dzsq=(hgtretc(k,ista)-hgtradc(klev,jsta)) *             &
                     (hgtretc(k,ista)-hgtradc(klev,jsta))
!
                IF(qualrad(1,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corret(1,k,ista)=corret(1,k,ista)+ &
                             (tabexp(INDEX)*ABS(uazmrad(jsta)))
                END IF
!
                IF(qualrad(2,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corret(2,k,ista)=corret(2,k,ista)+ &
                             (tabexp(INDEX)*ABS(vazmrad(jsta)))
                END IF
!
                IF(qualrad(3,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corret(5,k,ista)=                 &
                      corret(5,k,ista)+tabexp(INDEX)
                END IF
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar covariance sum
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xradc(jsta))*(xradc(ista)-xradc(jsta))    &
                +(yradc(ista)-yradc(jsta))*(yradc(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtradc(klev,jsta)) *             &
                     (hgtradc(k,ista)-hgtradc(klev,jsta))
!
                IF(qualrad(1,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(1,k,ista)=corrad(1,k,ista)+ &
                             (tabexp(INDEX)*ABS(uazmrad(jsta)))
                END IF
!
                IF(qualrad(2,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(2,k,ista)=corrad(2,k,ista)+ &
                             (tabexp(INDEX)*ABS(vazmrad(jsta)))
                END IF
!
                IF(qualrad(3,klev,jsta) > 0) THEN
                  rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) corrad(3,k,ista)=                 &
                      corrad(3,k,ista)+tabexp(INDEX)
                END IF
              END DO
            END IF
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF   ! radsw check
!
!-----------------------------------------------------------------------
!
!  Take inverse to speed up later calculations or
!  set covariance inverse to zero as a mask for "bad" obs
!
!-----------------------------------------------------------------------
!
  DO ista=1,nobsng
    DO ivar=1,nvar
      IF(qualsng(ivar,ista) > 0 ) THEN
        corsng(ivar,ista)=1./corsng(ivar,ista)
      ELSE
        corsng(ivar,ista)=0.
      END IF
    END DO
  END DO
!
  DO ista=1,nobsua
    DO k=1,nlevsua(ista)
      DO ivar=1,nvar
        IF(qualua(ivar,k,ista) > 0 ) THEN
          corua(ivar,k,ista)=1./corua(ivar,k,ista)
        ELSE
          corua(ivar,k,ista)=0.
        END IF
      END DO
    END DO
  END DO
!
  DO ista=1,ncolret
    DO k=1,nlevret(ista)
      DO ivar=1,nvar
        IF(qualret(ivar,k,ista) > 0 ) THEN
          corret(ivar,k,ista)=1./corret(ivar,k,ista)
        ELSE
          corret(ivar,k,ista)=0.
        END IF
      END DO
    END DO
  END DO
!
  DO ista=1,ncolrad
    DO k=1,nlevrad(ista)
      DO ivar=1,nvarrad
        IF(qualrad(ivar,k,ista) > 0 ) THEN
          IF(corrad(ivar,k,ista) == 0.) THEN
            PRINT *, ' zero corrad: ',ivar,k,ista
            PRINT *, ' qobsrad: ',qobsrad(ivar,k,ista)
            STOP
          END IF
          corrad(ivar,k,ista)=1./corrad(ivar,k,ista)
        ELSE
          corrad(ivar,k,ista)=0.
        END IF
      END DO
    END DO
  END DO
!
  istatus=0
  RETURN
END SUBROUTINE obscor
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BRTOBS3d                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE brtobs3d(nvar,nvarrad,nzua,nzrdr,nzret,                      & 2
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xcor,                    &
           odifsng,qobsng,corsng,xsng,ysng,hgtsng,                      &
           isrcsng,icatsng,nobsng,                                      &
           odifua,qobsua,corua,xua,yua,hgtua,isrcua,nlevsua,nobsua,     &
           odifrad,qobsrad,corrad,xradc,yradc,hgtradc,                  &
           uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,                &
           odifret,qobsret,corret,xretc,yretc,hgtretc,                  &
           iret,isrcret,nlevret,ncolret,                                &
           iusesng,iuseua,iuserad,iuseret,                              &
           sngsw,uasw,radsw,retsw,kpvrsq,rpass,                         &
           wlim,zrngsq,rngsqi,                                          &
           oanxsng,oanxua,oanxrad,oanxret,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Objectively analyze data at observation pts
!  Bratseth scheme is employed.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  July, 1995
!
!  MODIFICATION HISTORY:
!
!  Jan, 1996 (KB)
!  Added radar data and other improvements.
!  Keith Brewster, CAPS
!
!  July, 1996 (KB)
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!
!  Jan 5, 1997 (KB)
!  New version to reorder summing to speed-up code.
!  Renamed to BRTGRD3D from BRTGRDZ3d and BRTGRDTH3d.
!  When called to analyze using delta-theta, potential
!  temperature data is passed into the height arrays.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Sizing arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nvar,nvarrad,nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
  REAL :: xcor(ncat,ncat)
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: odifsng(nvar,mxsng)
  REAL :: qobsng(nvar,mxsng)
  REAL :: corsng(nvar,mxsng)
  REAL :: xsng(mxsng)
  REAL :: ysng(mxsng)
  REAL :: hgtsng(mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: icatsng(mxsng)
  INTEGER :: nobsng
!
  REAL :: odifua(nvar,nzua,mxua)
  REAL :: qobsua(nvar,nzua,mxua)
  REAL :: corua(nvar,nzua,mxua)
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: qobsrad(nvarrad,nzrdr,mxcolrad)
  REAL :: corrad(nvarrad,nzrdr,mxcolrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: qobsret(nvar,nzret,mxcolret)
  REAL :: corret(nvar,nzret,mxcolret)
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Scratch space
!
!-----------------------------------------------------------------------
!
  REAL :: rngsqi(nvar)
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: iusesng(0:nsrcsng)
  INTEGER :: iuseua(0:nsrcua)
  INTEGER :: iuserad(0:nsrcrad)
  INTEGER :: iuseret(0:nsrcret)
  INTEGER :: sngsw,uasw,radsw,retsw
  REAL :: kpvrsq(nvar)
  REAL :: rpass,                                                        &
       wlim
  REAL :: zrngsq
!
!-----------------------------------------------------------------------
!
!  Output arguments
!
!-----------------------------------------------------------------------
!
  REAL :: oanxsng(nvar,mxsng)
  REAL :: oanxua(nvar,nzua,mxua)
  REAL :: oanxrad(nvarrad,nzrdr,mxcolrad)
  REAL :: oanxret(nvar,nzret,mxcolret)
  INTEGER :: istatus
!
  REAL :: ftabinv
  INTEGER :: ntabexp
  PARAMETER (ntabexp=1000)
  REAL :: tabexp(ntabexp)
  COMMON /ftabexp/ ftabinv
  COMMON /tablexp/ tabexp
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: rsq,rlim,rlimsq,zsqinv
  REAL :: dist,dzsq,rnorm
  INTEGER :: ista,jsta,k,klev,ivar,INDEX
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rlimsq=0.
  DO ivar=1,nvar
    rsq=(kpvrsq(ivar)*rpass)
    rngsqi(ivar)=1./rsq
    rlimsq=AMAX1(rlimsq,rsq)
  END DO
  rlimsq=-rlimsq*ALOG(wlim)
  rlim=0.001*SQRT(rlimsq)
  WRITE(6,'(a,f10.2,a)') ' Influence cutoff radius: ',rlim,' km.'
!
  zsqinv=1./zrngsq
!
!-----------------------------------------------------------------------
!
!  Sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
  IF(sngsw > 0) THEN
    DO jsta=1,nobsng
      IF(iusesng(isrcsng(jsta)) > 0) THEN
        DO ivar=1,nvar
!
!-----------------------------------------------------------------------
!
!  Add to single-level data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsng
            dist=(xsng(ista)-xsng(jsta))*(xsng(ista)-xsng(jsta))        &
                +(ysng(ista)-ysng(jsta))*(ysng(ista)-ysng(jsta))
            IF(dist < rlimsq) THEN
              dzsq=(hgtsng(ista)-hgtsng(jsta)) *                        &
                   (hgtsng(ista)-hgtsng(jsta))
              rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxsng(ivar,ista)=oanxsng(ivar,ista)+ &
                  (tabexp(INDEX)*xcor(icatsng(ista),icatsng(jsta))      &
                   *odifsng(ivar,jsta)*corsng(ivar,jsta))
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsua
            dist=(xua(ista)-xsng(jsta))*(xua(ista)-xsng(jsta))          &
                +(yua(ista)-ysng(jsta))*(yua(ista)-ysng(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevsua(ista)
                dzsq=(hgtua(k,ista)-hgtsng(jsta)) *                     &
                     (hgtua(k,ista)-hgtsng(jsta))
                rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxua(ivar,k,ista)=oanxua(ivar,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifsng(ivar,jsta)*corsng(ivar,jsta))
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolret
            dist=(xretc(ista)-xsng(jsta))*(xretc(ista)-xsng(jsta))      &
                +(yretc(ista)-ysng(jsta))*(yretc(ista)-ysng(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevret(ista)
                dzsq=(hgtretc(k,ista)-hgtsng(jsta)) *                   &
                     (hgtretc(k,ista)-hgtsng(jsta))
                rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxret(ivar,k,ista)=oanxret(ivar,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifsng(ivar,jsta)*corsng(ivar,jsta))
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add contribution due to observation error variance
!
!-----------------------------------------------------------------------
!
          oanxsng(ivar,jsta)=                                           &
                 oanxsng(ivar,jsta)+qobsng(ivar,jsta)*                  &
                 odifsng(ivar,jsta)*corsng(ivar,jsta)
        END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar data analyses
!
!-----------------------------------------------------------------------
!
        DO ista=1,ncolrad
          dist=(xradc(ista)-xsng(jsta))*(xradc(ista)-xsng(jsta))        &
              +(yradc(ista)-ysng(jsta))*(yradc(ista)-ysng(jsta))
          IF(dist < rlimsq) THEN
            DO k=1,nlevrad(ista)
              dzsq=(hgtradc(k,ista)-hgtsng(jsta)) *                     &
                   (hgtradc(k,ista)-hgtsng(jsta))
!
              rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxrad(1,k,ista)=oanxrad(1,k,ista)+  &
                   (tabexp(INDEX) *                                     &
                    odifsng(1,jsta)*corsng(1,jsta))
!
              rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxrad(2,k,ista)=oanxrad(2,k,ista)+  &
                   (tabexp(INDEX) *                                     &
                    odifsng(2,jsta)*corsng(2,jsta))
!
              rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxrad(3,k,ista)=oanxrad(3,k,ista)+  &
                   (tabexp(INDEX) *                                     &
                    odifsng(5,jsta)*corsng(5,jsta))
            END DO
          END IF
        END DO
      END IF     ! source check
    END DO
  END IF   ! sngsw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
  IF(uasw > 0) THEN
    DO jsta=1,nobsua
      IF(iuseua(isrcua(jsta)) > 0) THEN
        DO klev=1,nlevsua(jsta)
          DO ivar=1,nvar
!
!-----------------------------------------------------------------------
!
!  Add to single-level data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsng
              dist=(xsng(ista)-xua(jsta))*(xsng(ista)-xua(jsta))        &
                  +(ysng(ista)-yua(jsta))*(ysng(ista)-yua(jsta))
              IF(dist < rlimsq) THEN
                dzsq=(hgtsng(ista)-hgtua(klev,jsta)) *                  &
                     (hgtsng(ista)-hgtua(klev,jsta))
                rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxsng(ivar,ista)=oanxsng(ivar,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifua(ivar,klev,jsta) *                          &
                       corua(ivar,klev,jsta))
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsua
              dist=(xua(ista)-xua(jsta))*(xua(ista)-xua(jsta))          &
                  +(yua(ista)-yua(jsta))*(yua(ista)-yua(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevsua(ista)
                  dzsq=(hgtua(k,ista)-hgtua(klev,jsta)) *               &
                       (hgtua(k,ista)-hgtua(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) oanxua(ivar,k,ista)=oanxua(ivar,k,ista)+ &
                       (tabexp(INDEX) *                                 &
                        odifua(ivar,klev,jsta) *                        &
                         corua(ivar,klev,jsta))
                END DO
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add term due to observation error
!
!-----------------------------------------------------------------------
!
            oanxua(ivar,klev,jsta)=                                     &
                oanxua(ivar,klev,jsta)+qobsua(ivar,klev,jsta)*          &
                odifua(ivar,klev,jsta)*corua(ivar,klev,jsta)
!
!-----------------------------------------------------------------------
!
!  Add to retrieval data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,ncolret
              dist=(xretc(ista)-xua(jsta))*(xretc(ista)-xua(jsta))      &
                  +(yretc(ista)-yua(jsta))*(yretc(ista)-yua(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevret(ista)
                  dzsq=(hgtretc(k,ista)-hgtua(klev,jsta)) *             &
                       (hgtretc(k,ista)-hgtua(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp)                                   &
                      oanxret(ivar,k,ista)=oanxret(ivar,k,ista)+        &
                       (tabexp(INDEX) *                                 &
                        odifua(ivar,klev,jsta) *                        &
                         corua(ivar,klev,jsta))
                END DO
              END IF
            END DO
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xua(jsta))*(xradc(ista)-xua(jsta))        &
                +(yradc(ista)-yua(jsta))*(yradc(ista)-yua(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtua(klev,jsta)) *               &
                     (hgtradc(k,ista)-hgtua(klev,jsta))
!
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(1,k,ista)=oanxrad(1,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifua(1,klev,jsta) *                             &
                       corua(1,klev,jsta))
!
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(2,k,ista)=oanxrad(2,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifua(2,klev,jsta) *                             &
                       corua(2,klev,jsta))
!
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(3,k,ista)=oanxrad(3,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifua(5,klev,jsta) *                             &
                       corua(5,klev,jsta))
              END DO
            END IF
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF   ! uasw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval data
!
!-----------------------------------------------------------------------
!
  IF(retsw > 0) THEN
    DO jsta=1,ncolret
      IF(iuseret(isrcret(iret(jsta))) > 0) THEN
        DO klev=1,nlevret(jsta)
          DO ivar=1,nvar
!
!-----------------------------------------------------------------------
!
!  Add to single-level data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsng
              dist=(xsng(ista)-xretc(jsta))*(xsng(ista)-xretc(jsta))    &
                  +(ysng(ista)-yretc(jsta))*(ysng(ista)-yretc(jsta))
              IF(dist < rlimsq) THEN
                dzsq=(hgtsng(ista)-hgtretc(klev,jsta)) *                &
                     (hgtsng(ista)-hgtretc(klev,jsta))
                rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxsng(ivar,ista)=oanxsng(ivar,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifret(ivar,klev,jsta) *                         &
                       corret(ivar,klev,jsta))
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,nobsua
              dist=(xua(ista)-xretc(jsta))*(xua(ista)-xretc(jsta))      &
                  +(yua(ista)-yretc(jsta))*(yua(ista)-yretc(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevsua(ista)
                  dzsq=(hgtua(k,ista)-hgtretc(klev,jsta)) *             &
                       (hgtua(k,ista)-hgtretc(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp) oanxua(ivar,k,ista)=oanxua(ivar,k,ista)+ &
                       (tabexp(INDEX) *                                 &
                        odifret(ivar,klev,jsta) *                       &
                         corret(ivar,klev,jsta))
                END DO
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval data analyses
!
!-----------------------------------------------------------------------
!
            DO ista=1,ncolret
              dist=(xretc(ista)-xretc(jsta))*(xretc(ista)-xretc(jsta))  &
                  +(yretc(ista)-yretc(jsta))*(yretc(ista)-yretc(jsta))
              IF(dist < rlimsq) THEN
                DO k=1,nlevret(ista)
                  dzsq=(hgtretc(k,ista)-hgtretc(klev,jsta)) *           &
                       (hgtretc(k,ista)-hgtretc(klev,jsta))
                  rnorm=(dist*rngsqi(ivar) +dzsq*zsqinv)*ftabinv
                  INDEX=nint(rnorm+1)
                  IF(INDEX < ntabexp)                                   &
                      oanxret(ivar,k,ista)=oanxret(ivar,k,ista)+        &
                       (tabexp(INDEX) *                                 &
                        odifret(ivar,klev,jsta) *                       &
                         corret(ivar,klev,jsta))
                END DO
              END IF
            END DO
!
!-----------------------------------------------------------------------
!
!  Add term due to observation error
!
!-----------------------------------------------------------------------
!
            oanxret(ivar,klev,jsta)=                                    &
                oanxret(ivar,klev,jsta)+qobsret(ivar,klev,jsta)*        &
                odifret(ivar,klev,jsta)*corret(ivar,klev,jsta)
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xretc(jsta))*(xradc(ista)-xretc(jsta))    &
                +(yradc(ista)-yretc(jsta))*(yradc(ista)-yretc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtretc(klev,jsta)) *             &
                     (hgtradc(k,ista)-hgtretc(klev,jsta))
!
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(1,k,ista)=oanxrad(1,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifret(1,klev,jsta) *                            &
                       corret(1,klev,jsta))
!
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(2,k,ista)=oanxrad(2,k,ista)+ &
                      (tabexp(INDEX) *                                  &
                       odifret(2,klev,jsta) *                           &
                        corret(2,klev,jsta))
!
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(3,k,ista)=oanxrad(3,k,ista)+ &
                     (tabexp(INDEX) *                                   &
                      odifret(5,klev,jsta) *                            &
                       corret(5,klev,jsta))
              END DO
            END IF
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF   ! retsw check
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
  IF(radsw > 0) THEN
    DO jsta=1,ncolrad
      IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
        DO klev=1,nlevrad(jsta)
!
!-----------------------------------------------------------------------
!
!  Add to single-level data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsng
            dist=(xsng(ista)-xradc(jsta))*(xsng(ista)-xradc(jsta))      &
                +(ysng(ista)-yradc(jsta))*(ysng(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              dzsq=(hgtsng(ista)-hgtradc(klev,jsta)) *                  &
                   (hgtsng(ista)-hgtradc(klev,jsta))
!
              rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxsng(1,ista)=oanxsng(1,ista)+      &
                   (tabexp(INDEX) * ABS(uazmrad(jsta)) *                &
                    odifrad(1,klev,jsta) *                              &
                     corrad(1,klev,jsta))
!
              rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxsng(2,ista)=oanxsng(2,ista)+      &
                    (tabexp(INDEX) * ABS(vazmrad(jsta)) *               &
                     odifrad(2,klev,jsta) *                             &
                      corrad(2,klev,jsta))
!
              rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
              INDEX=nint(rnorm+1)
              IF(INDEX < ntabexp) oanxsng(5,ista)=oanxsng(5,ista)+      &
                   (tabexp(INDEX) *                                     &
                    odifrad(3,klev,jsta) *                              &
                     corrad(3,klev,jsta))
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to multiple-level data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,nobsua
            dist=(xua(ista)-xradc(jsta))*(xua(ista)-xradc(jsta))        &
                +(yua(ista)-yradc(jsta))*(yua(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevsua(ista)
                dzsq=(hgtua(k,ista)-hgtradc(klev,jsta)) *               &
                     (hgtua(k,ista)-hgtradc(klev,jsta))
!
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxua(1,k,ista)=oanxua(1,k,ista)+  &
                     (tabexp(INDEX) * ABS(uazmrad(jsta)) *              &
                      odifrad(1,klev,jsta) *                            &
                       corrad(1,klev,jsta))
!
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxua(2,k,ista)=oanxua(2,k,ista)+  &
                     (tabexp(INDEX) * ABS(vazmrad(jsta)) *              &
                      odifrad(2,klev,jsta) *                            &
                       corrad(2,klev,jsta))
!
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxua(5,k,ista)=oanxua(5,k,ista)+  &
                     (tabexp(INDEX) *                                   &
                      odifrad(3,klev,jsta) *                            &
                       corrad(3,klev,jsta))
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to retrieval data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolret
            dist=(xretc(ista)-xradc(jsta))*(xretc(ista)-xradc(jsta))    &
                +(yretc(ista)-yradc(jsta))*(yretc(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevret(ista)
                dzsq=(hgtretc(k,ista)-hgtradc(klev,jsta)) *             &
                     (hgtretc(k,ista)-hgtradc(klev,jsta))
!
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxret(1,k,ista)=oanxret(2,k,ista)+ &
                      (tabexp(INDEX) * ABS(uazmrad(jsta)) *             &
                         odifrad(1,klev,jsta) *                         &
                          corrad(1,klev,jsta))
!
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxret(2,k,ista)=oanxret(2,k,ista)+ &
                      (tabexp(INDEX) * ABS(vazmrad(jsta)) *             &
                         odifrad(2,klev,jsta) *                         &
                          corrad(2,klev,jsta))
!
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxret(5,k,ista)=oanxret(5,k,ista)+ &
                      (tabexp(INDEX) *                                  &
                         odifrad(3,klev,jsta) *                         &
                          corrad(3,klev,jsta))
              END DO
            END IF
          END DO
!
!-----------------------------------------------------------------------
!
!  Add to radar data analyses
!
!-----------------------------------------------------------------------
!
          DO ista=1,ncolrad
            dist=(xradc(ista)-xradc(jsta))*(xradc(ista)-xradc(jsta))    &
                +(yradc(ista)-yradc(jsta))*(yradc(ista)-yradc(jsta))
            IF(dist < rlimsq) THEN
              DO k=1,nlevrad(ista)
                dzsq=(hgtradc(k,ista)-hgtradc(klev,jsta)) *             &
                     (hgtradc(k,ista)-hgtradc(klev,jsta))
!
                rnorm=(dist*rngsqi(1) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(1,k,ista)=oanxrad(1,k,ista)+ &
                    (tabexp(INDEX) * ABS(uazmrad(jsta)) *               &
                    odifrad(1,klev,jsta) *                              &
                     corrad(1,klev,jsta))
!
                rnorm=(dist*rngsqi(2) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(2,k,ista)=oanxrad(2,k,ista)+ &
                      (tabexp(INDEX) * ABS(vazmrad(jsta)) *             &
                       odifrad(2,klev,jsta) *                           &
                        corrad(2,klev,jsta))
!
                rnorm=(dist*rngsqi(5) +dzsq*zsqinv)*ftabinv
                INDEX=nint(rnorm+1)
                IF(INDEX < ntabexp) oanxrad(3,k,ista)=oanxrad(3,k,ista)+ &
                      (tabexp(INDEX) *                                  &
                       odifrad(3,klev,jsta) *                           &
                        corrad(3,klev,jsta))
              END DO
            END IF
          END DO
          oanxrad(1,klev,jsta)=                                         &
                 oanxrad(1,klev,jsta)+qobsrad(1,klev,jsta)*             &
                 odifrad(1,klev,jsta)*corrad(1,klev,jsta)
          oanxrad(2,klev,jsta)=                                         &
                 oanxrad(2,klev,jsta)+qobsrad(2,klev,jsta)*             &
                 odifrad(2,klev,jsta)*corrad(2,klev,jsta)
          oanxrad(3,klev,jsta)=                                         &
                 oanxrad(3,klev,jsta)+qobsrad(3,klev,jsta)*             &
                 odifrad(3,klev,jsta)*corrad(3,klev,jsta)
        END DO
      END IF     ! source check
    END DO
  END IF   ! radsw check
  istatus=0
  RETURN
END SUBROUTINE brtobs3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BRTGRD3d                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE brtgrd3d(nx,ny,nz,nvar,nvarrad,                              & 2
           nzua,nzrdr,nzret,                                            &
           mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret,                    &
           nsrcsng,nsrcua,nsrcrad,nsrcret,ncat,xs,ys,zs,icatg,xcor,     &
           odifsng,corsng,xsng,ysng,hgtsng,isrcsng,icatsng,nobsng,      &
           odifua,corua,xua,yua,hgtua,isrcua,nlevsua,nobsua,            &
           odifrad,corrad,xradc,yradc,hgtradc,                          &
           uazmrad,vazmrad,irad,isrcrad,nlevrad,ncolrad,                &
           odifret,corret,xretc,yretc,hgtretc,                          &
           iret,isrcret,nlevret,ncolret,                                &
           iusesng,iuseua,iuserad,iuseret,                              &
           sngsw,uasw,radsw,retsw,kpvrsq,rpass,                         &
           wlim,zrngsq,rngsqi,                                          &
           anx,tem1,tem2,ibegin,iend,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Objectively analyze data on a regular grid specified
!  by x and y.
!  Bratseth scheme is employed.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  July, 1995
!
!  MODIFICATION HISTORY:
!
!  Jan, 1996 (KB)
!  Added radar data and other improvements.
!
!  July, 1996 (KB)
!  Added sngsw and uasw to switch on/off use of single-level
!  and sounding data, respectively.
!
!  Jan 5, 1997 (KB)
!  New version to reorder summing to speed-up code.
!  Renamed to BRTGRD3D from BRTGRDZ3d and BRTGRDTH3d.
!  When called to analyze using delta-theta, potential
!  temperature data is passed into the height arrays.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Sizing arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  INTEGER :: nvar,nvarrad,nzua,nzrdr,nzret
  INTEGER :: mxsng,mxua,mxrad,mxcolrad,mxret,mxcolret
  INTEGER :: nsrcsng,nsrcua,nsrcrad,nsrcret,ncat
!
!-----------------------------------------------------------------------
!
!  Grid arguments
!
!-----------------------------------------------------------------------
!
  REAL :: xs(nx)
  REAL :: ys(ny)
  REAL :: zs(nx,ny,nz)
  INTEGER :: icatg(nx,ny)
  REAL :: xcor(ncat,ncat)
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  REAL :: odifsng(nvar,mxsng)             ! variable to analyse
  REAL :: corsng(nvar,mxsng)
  REAL :: xsng(mxsng)
  REAL :: ysng(mxsng)
  REAL :: hgtsng(mxsng)
  INTEGER :: isrcsng(mxsng)
  INTEGER :: icatsng(mxsng)
  INTEGER :: nobsng
!
  REAL :: odifua(nvar,nzua,mxua)
  REAL :: corua(nvar,nzua,mxua)
  REAL :: xua(mxua)
  REAL :: yua(mxua)
  REAL :: hgtua(nzua,mxua)
  INTEGER :: nlevsua(mxua)
  INTEGER :: isrcua(mxua)
  INTEGER :: nobsua
!
  REAL :: odifrad(nvarrad,nzrdr,mxcolrad)
  REAL :: corrad(nvarrad,nzrdr,mxcolrad)
  REAL :: xradc(mxcolrad)
  REAL :: yradc(mxcolrad)
  REAL :: hgtradc(nzrdr,mxcolrad)
  REAL :: uazmrad(mxcolrad)
  REAL :: vazmrad(mxcolrad)
  INTEGER :: nlevrad(mxcolrad)
  INTEGER :: irad(mxcolrad)
  INTEGER :: isrcrad(0:mxrad)
  INTEGER :: ncolrad
!
  REAL :: odifret(nvar,nzret,mxcolret)
  REAL :: corret(nvar,nzret,mxcolret)
  REAL :: xretc(mxcolret)
  REAL :: yretc(mxcolret)
  REAL :: hgtretc(nzret,mxcolret)
  INTEGER :: nlevret(mxcolret)
  INTEGER :: iret(mxcolret)
  INTEGER :: isrcret(0:mxret)
  INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
!  Scratch space
!
!-----------------------------------------------------------------------
!
  REAL :: rngsqi(nvar)
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: iusesng(0:nsrcsng)
  INTEGER :: iuseua(0:nsrcua)
  INTEGER :: iuserad(0:nsrcrad)
  INTEGER :: iuseret(0:nsrcret)
  INTEGER :: sngsw,uasw,radsw,retsw
  REAL :: kpvrsq(nvar)
  REAL :: rpass,                                                        &
       wlim
  REAL :: zrngsq
!
!-----------------------------------------------------------------------
!
!  Output grid arguments
!
!-----------------------------------------------------------------------
!
  REAL :: anx(nx,ny,nz,nvar)
  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
  INTEGER :: ibegin(ny)
  INTEGER :: iend(ny)
  INTEGER :: istatus
!
  REAL :: ftabinv
  INTEGER :: ntabexp
  PARAMETER (ntabexp=1000)
  REAL :: tabexp(ntabexp)
  COMMON /ftabexp/ ftabinv
  COMMON /tablexp/ tabexp
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: rsq,rlim,rlimsq,zsqinv,rnorm
  INTEGER :: i,j,k,jsta,klev,ivar,INDEX,jbegin,jend
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rlimsq=0.
  DO ivar=1,nvar
    rsq=(kpvrsq(ivar)*rpass)
    rngsqi(ivar)=1./rsq
    rlimsq=AMAX1(rlimsq,rsq)
  END DO
  rlimsq=-rlimsq*ALOG(wlim)
  rlim=0.001*SQRT(rlimsq)
  WRITE(6,'(/a,f10.2,a/)') ' Influence cutoff radius: ',rlim,' km.'
!
  zsqinv=1./zrngsq
!
!-----------------------------------------------------------------------
!
!  First sum contributions from single-level obs
!
!-----------------------------------------------------------------------
!
  IF(sngsw > 0) THEN
    DO jsta=1,nobsng
      IF(iusesng(isrcsng(jsta)) > 0) THEN
        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=(xs(i)-xsng(jsta))*(xs(i)-xsng(jsta))           &
                       +(ys(j)-ysng(jsta))*(ys(j)-ysng(jsta))
          END DO
        END DO
!
!-----------------------------------------------------------------------
!
!  Find range in i,j of valid influence
!
!-----------------------------------------------------------------------
!
        jbegin=ny+1
        jend=1
        DO j=1,ny
          DO i=1,nx
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          ibegin(j)=i
          IF(i <= nx) THEN
            jbegin=MIN(jbegin,j)
            jend=MAX(jend,j)
          END IF
          DO i=nx,ibegin(j),-1
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          iend(j)=i
        END DO
!
! OpenMP changed loop order to j,k,i:
!$OMP PARALLEL DO PRIVATE (j,k,i)
        DO j=jbegin,jend
          DO k=1,nz-1
            DO i=ibegin(j),iend(j)
              tem2(i,j,k)=(zs(i,j,k)-hgtsng(jsta)) *                    &
                          (zs(i,j,k)-hgtsng(jsta)) * zsqinv
            END DO
          END DO
        END DO
        DO ivar=1,nvar
! OpenMP changed loop order to j,k,i:
!$OMP PARALLEL DO PRIVATE (j,k,i,rnorm,index)
          DO j=jbegin,jend
            DO k=1,nz-1
              DO i=ibegin(j),iend(j)
                rnorm=(tem1(i,j,1)*rngsqi(ivar) +                       &
                       tem2(i,j,k))*ftabinv
                INDEX=                                                  &
                    MIN(nint(rnorm+1),ntabexp)
                anx(i,j,k,ivar)=anx(i,j,k,ivar) +                       &
                    (tabexp(INDEX)*xcor(icatsng(jsta),icatg(i,j))       &
                    *odifsng(ivar,jsta)*corsng(ivar,jsta))
              END DO
            END DO
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from multiple-level obs
!
!-----------------------------------------------------------------------
!
  IF(uasw > 0) THEN
    DO jsta=1,nobsua
      IF(iuseua(isrcua(jsta)) > 0) THEN
        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=(xs(i)-xua(jsta))*(xs(i)-xua(jsta))             &
                       +(ys(j)-yua(jsta))*(ys(j)-yua(jsta))
          END DO
        END DO
!
!-----------------------------------------------------------------------
!
!  Find range in i,j of valid influence
!
!-----------------------------------------------------------------------
!
        jbegin=ny+1
        jend=1
        DO j=1,ny
          DO i=1,nx
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          ibegin(j)=i
          IF(i <= nx) THEN
            jbegin=MIN(jbegin,j)
            jend=MAX(jend,j)
          END IF
          DO i=nx,ibegin(j),-1
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          iend(j)=i
        END DO
!
        DO klev=1,nlevsua(jsta)
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                tem2(i,j,k)=(zs(i,j,k)-hgtua(klev,jsta)) *              &
                            (zs(i,j,k)-hgtua(klev,jsta)) * zsqinv
              END DO
            END DO
          END DO
          DO ivar=1,nvar
            DO k=1,nz-1
              DO j=jbegin,jend
                DO i=ibegin(j),iend(j)
                  rnorm=(tem1(i,j,1)*rngsqi(ivar) +                     &
                         tem2(i,j,k) ) * ftabinv
                  INDEX=                                                &
                      MIN(nint(rnorm+1),ntabexp)
                  anx(i,j,k,ivar)=anx(i,j,k,ivar) +                     &
                      (tabexp(INDEX) *                                  &
                       odifua(ivar,klev,jsta) *                         &
                       corua(ivar,klev,jsta))
                END DO
              END DO
            END DO
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from radar obs
!
!-----------------------------------------------------------------------
!
  IF(radsw > 0) THEN
    DO jsta=1,ncolrad
      IF(iuserad(isrcrad(irad(jsta))) > 0) THEN
        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=(xs(i)-xradc(jsta))*(xs(i)-xradc(jsta))         &
                       +(ys(j)-yradc(jsta))*(ys(j)-yradc(jsta))
          END DO
        END DO
!
!-----------------------------------------------------------------------
!
!  Find range in i,j of valid influence
!
!-----------------------------------------------------------------------
!
        jbegin=ny+1
        jend=1
        DO j=1,ny
          DO i=1,nx
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          ibegin(j)=i
          IF(i <= nx) THEN
            jbegin=MIN(jbegin,j)
            jend=MAX(jend,j)
          END IF
          DO i=nx,ibegin(j),-1
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          iend(j)=i
        END DO
!
        DO klev=1,nlevrad(jsta)
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                tem2(i,j,k)=(zs(i,j,k)-hgtradc(klev,jsta)) *            &
                            (zs(i,j,k)-hgtradc(klev,jsta)) * zsqinv
              END DO
            END DO
          END DO
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                rnorm=(tem1(i,j,1)*rngsqi(1) +                          &
                       tem2(i,j,k) ) * ftabinv
                INDEX=                                                  &
                    MIN(nint(rnorm+1),ntabexp)
                anx(i,j,k,1)=anx(i,j,k,1) +                             &
                    (tabexp(INDEX) *                                    &
                     ABS(uazmrad(jsta)) *                               &
                     odifrad(1,klev,jsta) *                             &
                     corrad(1,klev,jsta))
              END DO
            END DO
          END DO
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                rnorm=(tem1(i,j,1)*rngsqi(2) +                          &
                       tem2(i,j,k) ) * ftabinv
                INDEX=                                                  &
                    MIN(nint(rnorm+1),ntabexp)
                anx(i,j,k,2)=anx(i,j,k,2) +                             &
                       (tabexp(INDEX) *                                 &
                         ABS(vazmrad(jsta)) *                           &
                         odifrad(2,klev,jsta) *                         &
                         corrad(2,klev,jsta))
              END DO
            END DO
          END DO
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                rnorm=(tem1(i,j,1)*rngsqi(5) +                          &
                       tem2(i,j,k) ) * ftabinv
                INDEX=                                                  &
                    MIN(nint(rnorm+1),ntabexp)
                anx(i,j,k,5)=anx(i,j,k,5) +                             &
                       (tabexp(INDEX) *                                 &
                         odifrad(3,klev,jsta) *                         &
                         corrad(3,klev,jsta))
              END DO
            END DO
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Add contributions from retrieval data
!
!-----------------------------------------------------------------------
!
  IF(retsw > 0) THEN
    DO jsta=1,ncolret
      IF(iuseret(isrcret(iret(jsta))) > 0) THEN
        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=(xs(i)-xretc(jsta))*(xs(i)-xretc(jsta))         &
                       +(ys(j)-yretc(jsta))*(ys(j)-yretc(jsta))
          END DO
        END DO
!
!-----------------------------------------------------------------------
!
!  Find range in i,j of valid influence
!
!-----------------------------------------------------------------------
!
        jbegin=ny+1
        jend=1
        DO j=1,ny
          DO i=1,nx
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          ibegin(j)=i
          IF(i <= nx) THEN
            jbegin=MIN(jbegin,j)
            jend=MAX(jend,j)
          END IF
          DO i=nx,ibegin(j),-1
            IF(tem1(i,j,1) < rlimsq) EXIT
          END DO
          iend(j)=i
        END DO
!
!-----------------------------------------------------------------------
!
!  Apply retrievel corrections
!
!-----------------------------------------------------------------------
!
        DO klev=1,nlevret(jsta)
          DO k=1,nz-1
            DO j=jbegin,jend
              DO i=ibegin(j),iend(j)
                tem2(i,j,k)=(zs(i,j,k)-hgtretc(klev,jsta)) *            &
                       (zs(i,j,k)-hgtretc(klev,jsta)) * zsqinv
              END DO
            END DO
          END DO
          DO ivar=1,nvar
            DO j=jbegin,jend
              DO k=1,nz-1
                DO i=ibegin(j),iend(j)
                  rnorm=(tem1(i,j,1)*rngsqi(ivar) +                     &
                         tem2(i,j,k) ) * ftabinv
                  INDEX=                                                &
                      MIN(nint(rnorm+1),ntabexp)
                  anx(i,j,k,ivar)=anx(i,j,k,ivar) +                     &
                           (tabexp(INDEX) *                             &
                             odifret(ivar,klev,jsta) *                  &
                             corret(ivar,klev,jsta))
                END DO
              END DO
            END DO
          END DO
        END DO
      END IF     ! source check
    END DO
  END IF
!
  istatus=0
  RETURN
END SUBROUTINE brtgrd3d
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE BARQC                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE barqc(maxsta,nvar,nsrc,nprev,nsta,                           & 1
           obs,obstime,xsta,ysta,isrc,qual,                             &
           knt,wgtsum,zsum,sqsum,                                       &
           range,wlim,klim,qclim,varnam,snam,                           &
           iqclist,iqcsave,iqcflag)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Routine to objectively quality control surface data using a
!  Barnes analysis of nearby stations.
!
!  If an observations is found to be bad its quality variable
!  is set to the input flag value, iqcflag
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  September, 1989
!
!  MODIFICATION HISTORY:
!
!  Jul, 1996 (KB)
!  Version for ARPS analysis
!  Keith Brewster, CAPS
!
!  Aug, 1997 (KB)
!  Modifications and clarification of code for QC of analysis
!  variables instead of raw observed variables, before.  Related
!  to changes in data error specification in other routines.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Observation Arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: maxsta,nvar,nsrc
  REAL :: obs(nvar,maxsta)
  INTEGER :: obstime(maxsta)
  REAL :: xsta(maxsta)
  REAL :: ysta(maxsta)
  INTEGER :: isrc(maxsta)
  INTEGER :: nprev
  INTEGER :: nsta
!
!-----------------------------------------------------------------------
!
!  Scratch space
!
!-----------------------------------------------------------------------
!
  INTEGER :: knt(nvar)
  REAL :: wgtsum(nvar)
  REAL :: zsum(nvar)
  REAL :: sqsum(nvar)
!
!-----------------------------------------------------------------------
!
!  File unit numbers
!
!-----------------------------------------------------------------------
!
  INTEGER :: iqclist,iqcsave
!
!-----------------------------------------------------------------------
!
!  Analysis specification arguments
!
!-----------------------------------------------------------------------
!
  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(nvar,nsrc)   ! QC cutoffs
!
!-----------------------------------------------------------------------
!
!  Input/Output quality indicator
!
!-----------------------------------------------------------------------
!
  INTEGER :: qual(nvar,maxsta)
!
!-----------------------------------------------------------------------
!
!  Diagnostics and other argument variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=6) :: varnam(nvar)
  CHARACTER (LEN=5) :: snam(maxsta)
  INTEGER :: iqcflag
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: rlimsq,dist,                                                  &
       wgt,zanx,dob,sqanx,chklim
  INTEGER :: jsta,ksta,ivar
  LOGICAL :: listit,saveit
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  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
!
!-----------------------------------------------------------------------
!
  PRINT *, ' BARQC: nsta= ',nsta
  DO ksta=(nprev+1),(nprev+nsta)
!    print *, ' ksta,isrc= ',ksta,isrc(ksta)
!    print *, '   obs: '
!    print *, (obs(ivar,ksta),ivar=1,nvar)
!    print *, '   qual: '
!    print *, (qual(ivar,ksta),ivar=1,nvar)
    IF(isrc(ksta) > 0) THEN
      DO ivar=1,nvar
        zsum(ivar)=0.
        sqsum(ivar)=0.
        wgtsum(ivar)=0.
        knt(ivar)=0
      END DO
      DO jsta=1,(nprev+nsta)
        IF(jsta /= ksta) THEN
          dist=(xsta(ksta)-xsta(jsta))*(xsta(ksta)-xsta(jsta))          &
              +(ysta(ksta)-ysta(jsta))*(ysta(ksta)-ysta(jsta))
          IF(dist < rlimsq) THEN
            wgt=EXP(-dist/range)
            DO ivar=1,nvar
              IF(qual(ivar,jsta) > 0) THEN
                knt(ivar)=knt(ivar)+1
                wgtsum(ivar)=wgtsum(ivar)+wgt
                zsum(ivar)=zsum(ivar)+wgt*obs(ivar,jsta)
                sqsum(ivar)=sqsum(ivar)+                                &
                     wgt*obs(ivar,jsta)*obs(ivar,jsta)
              END IF
            END DO
          END IF
        END IF  ! dont use current station
      END DO
      DO ivar=1,nvar
        IF(qual(ivar,ksta) > 0) 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=obs(ivar,ksta)-zanx
            chklim=AMAX1(sqanx,qclim(ivar,isrc(ksta)))
            IF(ABS(dob) > chklim) THEN
              PRINT *, ' QC flagging data...sta=',snam(ksta)
              PRINT *, '     var, ob =',                                &
                          varnam(ivar),obs(ivar,ksta)
              PRINT *, '     delta ob =',dob
              PRINT *, '  3.0*stdv,qclimit=',sqanx,                     &
                                   qclim(ivar,isrc(ksta))
              IF(listit) WRITE(iqclist,810)                             &
                  snam(ksta),obstime(ksta),ivar,varnam(ivar),           &
                  obs(ivar,ksta),dob,                                   &
                  sqanx,qclim(ivar,isrc(ksta))
              810           FORMAT(2X,a5,i5,2X,i3,1X,a6,f11.2,f11.2,f11.2,f11.2)
              qual(ivar,ksta)=iqcflag
            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),           &
                (obs(ivar,ksta),ivar=1,nvar)
      820   FORMAT(2X,a5,i5,2X,15(f7.2))
    END IF     ! isrc check
  END DO
  RETURN
END SUBROUTINE barqc

!
!##################################################################
!##################################################################
!######                                                      ######
!######                   FUNCTION TEXP                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION texp(rnorm)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  A faster exponential function, via a lookup table.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:  Keith Brewster, CAPS
!  December, 1995
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  REAL :: texp,setexp
  REAL :: rnorm
!
  REAL :: mxrnrm,ftabinv
  INTEGER :: ntabexp
  PARAMETER (ntabexp=1000)
  REAL :: tabexp(ntabexp)
  COMMON /ftabexp/ ftabinv
  COMMON /tablexp/ tabexp
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,INDEX
  REAL :: ftab
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  INDEX=nint(-(rnorm*ftabinv)) + 1
  INDEX=MAX(INDEX,1)
  INDEX=MIN(INDEX,ntabexp)
  texp=tabexp(INDEX)
  RETURN
!
  ENTRY setexp(mxrnrm)
  ftab=mxrnrm/FLOAT(ntabexp)
  ftabinv=1./ftab
  DO i=1,ntabexp
    tabexp(i)=EXP((-(i-1)*ftab))
  END DO
  setexp=ftabinv
  RETURN
  END FUNCTION texp