cmm$  Be sure that you have correct vterpv, vterpv_v1,vterpv_v2 etc
c*************************************
CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of DELHUV from delhuv.f
CMAX: ----------------------------------------------------
c......................................................................
c                                                                     .
c  PROGRAM:    DELHUV                                                 .
c  LANGUAGE:   Cray Fortran 77 Version 6                              .
c  COMPUTER:   Cray C90                                               .
c                                                                     .
c  PROGRAMMER: M. SEABLOM                                             .
c  DATE:       November 2, 1993                                       .
c                                                                     .
c  Program DELHUV performs version 2.0 quality control and computes   .
c  deltas for upper air height and wind data.                         .
c                                                                     .
c  Program modifications:                                             .
c    - 11/90 (Seablom)  Corrected error in OBERTRP3                   .
c    -  1/91 (Lamich)   Rejecting observations below ground           .
c                                                                     .
c  Modifications for version 1.1:                                     .
c    -  3/91 (Seablom)  Expanded DATA array to allow >100K obs.       .
c    -  3/91 (Seablom)  Reports tallied from DELHUV rather than RDHUV..
c    -  4/91 (Seablom)  Added date/time stamp to DATA array.          .
c    -  4/91 (Seablom)  Level-dependent tolerance values included.    .
c                       Wind tolerances increased significantly over  .
c                       OI version 1.0.                               .
c                                                                     .
c  Modifications for version 1.2:                                     .
c    - 11/91 (Seablom)  Removal of data ingest.                       .
c                                                                     .
c  Modifications for version 1.3:                                     .
c    -  7/92 (Seablom)  Corrected bug in wind quality control.        .
c    -  7/92 (Seablom)  Code modifications to allow parallelization.  .
c    -  7/92 (Seablom)  Wind standard deviations passed back to driver.
c                                                                     .
c  Modifications for version 1.4:                                     .
c    - 11/93 (Seablom)  Explicitly designated variables for parallel- .
c                       ization.                                      .
c    - 11/93 (Lamich)   Properly passed "suspect" array.              .
c                                                                     .
c  Modifications for version 1.5:                                     .
c    - 11/93 (Seablom)  Implemented "bit vector" syntax for wind      .
c                       quality control following reproducibility     .
c                       problems in parallel environment.             .
c......................................................................
C                                                                     !
C NASA Goddard Space Flight Center       Data Assimilation Office     !
C                                                                     !
C Goddard Earth Observing System                                      !
C Data Analysis System - Version 1.5                                  !
C October, 1993                                                       !
C                                                                     !
c---------------------------------------------------------------------c
      SUBROUTINE delhuv(zmod,umod,vmod,tmod,psmod,z1000,fcerr,u_error,
     &    v_error,data,nobs,ze,zelat,zelon,nze,nymd,nhms,suspect)
      include '/usr/include/cm/CMF_defs.h'
      
c
c     Maximal observation points
      INTEGER MaxObservPoints
c
      PARAMETER (MaxObservPoints = 150000)
      PARAMETER (MaxMlev = 150)
c
c
      REAL zmod (im * jnp, mlev), umod (im * jnp * mlev), u_error (im, j
     &np, mlev), v_error (im, jnp, mlev), psmod (im * jnp), z1000 (im * 
     &jnp), prslev (12), fcerr (im * jnp * mlev), tmod (im * jnp, mlev),
     & vmod (im * jnp * mlev), lat, lon, kx, latmax, latmin, lonmax, lon
     &min, londif, ze (20000), zelat (20000), zelon (20000), srchlat, sr
     &chlon, dist, data (MaxObservPoints * 10), fg1 (MaxObservPoints), 
     &fg2 (MaxObservPoints), psobs (MaxObservPoints), ztol (MaxMlev), 
     &fg3 (MaxObservPoints)
      CHARACTER * 8 dataname, obsname
      CHARACTER * 60 message
      LOGICAL polar, suspect (MaxObservPoints, 2), qobs, offlev, onestag
     &e, dateline_case
      INTEGER gcfpoint (30000)
      INTEGER undergnd_count
      COMMON / dataptrs / ikx, ilat, ilon, iprs, ityp, idel, ife
      COMMON / gridxx / im, jnp, mlev
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / print / msglev
      COMMON / obs / maxdex, nfgge, lun, hh (13, 13, 20),
     &vcoroe (18, 18, 4), obserr (200, 4, 19), ww (13, 13, 20),
     &qobs (200, 5), offlev (200)
      common / obschar / dataname (5), obsname(200)
!hpf$ distribute data(block)
!hpf$ distribute fcerr(block)
!hpf$ distribute fg1(block)
!hpf$ distribute fg2(block)
!hpf$ distribute fg3(block)
!hpf$ distribute gcfpoint(block)
!hpf$ distribute hh(block,block,block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute offlev(block)
!hpf$ distribute prslev(block)
!hpf$ distribute psmod(block)
!hpf$ distribute psobs(block)
!hpf$ distribute qobs(block,block)
!hpf$ distribute suspect(block,block)
!hpf$ distribute tmod(block,block)
!hpf$ distribute umod(block)
!hpf$ distribute u_error(block,block,block)
!hpf$ distribute vcoroe(block,block,block)
!hpf$ distribute vmod(block)
!hpf$ distribute v_error(block,block,block)
!hpf$ distribute ww(block,block,block)
!hpf$ distribute z1000(block)
!hpf$ distribute ze(block)
!hpf$ distribute zelat(block)
!hpf$ distribute zelon(block)
!hpf$ distribute zmod(block,block)
!hpf$ distribute ztol(block)
c   Initialize data.
      LOGICAL lif
      LOGICAL lif100
      REAL  TOLERANCE100(nobs)
!hpf$ distribute TOLERANCE100(block)
      LOGICAL lif101
      LOGICAL lif102
      DATA fill / 1.e10 /, alpha / 9.e10 /, rade / 6.37e6 /
      INTEGER Ti1
      INTEGER Tztol_auto
      pi = 4. * atan(1.)
      drc = pi / 180.
      rdc = 180. / pi
      message(1:60) = ' '
      undergnd_count = 0
      nreacc = 0
c   Initialize tolerance arrays.
      ztol(mlev) = 2.8
      DO 7000 itol = mlev - 1 , 1 , (-1)
        ztol( itol ) = ztol( itol+1 ) + .04
 7000 CONTINUE
c  Compute the geostrophic ind errors.
      CALL werrinit(fcerr,im,jnp,mlev,u_error,v_error)
      WRITE( 6,* ) ' DELHUV: Total number of observations=',nobs
      IF (nobs .EQ. 0) THEN
        RETURN
      ENDIF
      DO ti1 = 1 , nobs
        fg2(ti1) = 0.
      ENDDO
c  Set pointers for winds, observation and forecast errors.
      iz = idel
      iu = iz + nobs
      iv = iu + nobs
      ife = iv + nobs
      ioe = ife + nobs
c  Establish major loop to perform the following:
c   1. Retrieve forecast and obs. errors and store in array "data".
c   2. Compute deltas.
c   3. Pass data through gross check.
c  Obtain first guess height values at observation locations.
      CALL vterpv_v1(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,iprs,maxobservpoints * 10,fg1,1,
     &    maxobservpoints,nobs,zmod)
c  Interpolate the forecast errors to observation locations.
      CALL vterpv_v1(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,iprs,maxobservpoints * 10,data,ife,
     &    maxobservpoints * 10,nobs,fcerr)
c  Interpolate observation error to observation location.
      CALL obertrp0_v1(data,ikx,maxobservpoints * 10,data,iprs,
     &    maxobservpoints * 10,data,ioe,maxobservpoints * 10,obserr,2,
     &    nobs)
c  Interpolate the model surface pressure to the observation
c  locations.
      CALL terpv_v2(data,ilat,maxobservpoints * 10,data,ilon,maxobservpo
     &    ints * 10,psobs,1,maxobservpoints,nobs,psmod)
c  QC gross check / delta computation for height data ----------------
      igcfail = 0
      itest = 0
      message = 'DELHUV: Impossible location; rejected point.'
cmic$ do all vector shared (data, drc, fg1, fg2, fill, ife, ikx, ilat)
cmic$*       shared (ilon, ioe, iprs, ityp, iz, kx, MaxObservPoints)
cmic$*       shared (message)
cmic$*       shared (mlev, nobs, psobs, suspect, undergnd_count)
cmic$*       shared (ztol)
cmic$*       private (delta, fe, i, lat, lon, oe, pointer, prs)
cmic$*       private (tolerance, trwt) autoscope
c  The interpolated surface pressure is used to eliminate observations
c  that are underground as far as the model is concerned.  This check
c  is quite crude in that it will sometimes keep and sometimes remove
c  observations in areas of rapidly varying topography.
c  Now compute the delta if valid data exists. Generated 1000mb height
c  reports are already stored as deltas.
c  Gross check the height data.
c  Tighten tolerance in the tropics.
      DO i = 1 , nobs
        IF (data(iz + i - 1) .LT. fill) THEN
          data(ityp + i - 1) = 1.
        ELSE
          data(ityp + i - 1) = 0.4
        ENDIF
        lat = data(ilat + i - 1)
        lon = data(ilon + i - 1)
        prs = data(iprs + i - 1)
        kx = data(ikx + i - 1)
        IF (lat .GT. 90. .OR. lat .LT. -90. .OR. lon .GT. 180. .OR. lon
     &        .LT. -180.) THEN
          data(ityp + i - 1) = 0.4
        ENDIF
        IF (prs .GT. psobs(i)) THEN
          data(ityp + i - 1) = 0.4
        ENDIF
        IF (data(ityp + i - 1) .GE. 1. .AND. data(ikx + i - 1) .NE. 87.
     &        ) THEN
          data(iz + i - 1) = data(iz + i - 1) - fg1(i)
        ENDIF
        delta = data(iz + i - 1)
        oe = data(ioe + i - 1)
        fe = data(ife + i - 1)
        IF (abs(lat) .LE. 30.0) THEN
          trwt = exp(sin(abs(lat) * drc)) / 1.65
        ELSE
          trwt = 1.
        ENDIF
        pointer = pindx2(prs)
        tolerance = ztol(int(pointer)) ** 2
        IF ((delta ** 2 .GT. (oe ** 2 + fe ** 2) * tolerance * trwt) 
     &        .AND. data(iz + i - 1) .LT. fill) THEN
          data(ityp + i - 1) = 0.1
          fg2(i) = 1.
        ENDIF
        IF (suspect(i,1) .AND. data(iz + i - 1) .LT. fill) THEN
          data(ityp + i - 1) = 0.3
          fg2(i) = 1.
        ENDIF
      ENDDO
celp       call flash( 3,message,1)
c           undergnd_count = undergnd_count + 1
c  Flag suspicious observations and subject them to the gross check.
      IF (undergnd_count .GT. 0) THEN
        WRITE( 6,1001 ) undergnd_count
1001  FORMAT ('DELHUV:',i5,' observations below model surface')
      ENDIF
      CALL wheneq(nobs,fg2,1,1.,gcfpoint,igcfail)
      WRITE( 6,* ) ' DELHUV: ',igcfail,' height observations have',
     &' failed the gross check.'
      IF (igcfail .GT. 30000) THEN
        message = ' DELHUV: Excessive gross check failures.'
        CALL flash(5,message,1)
      ENDIF
c  Find neighbouring observations to those flagged by the gross check.
c  Perform an interpolation of the good quality neighbours (within a
c  1000km radius) to the location of a suspicious report.
cmic$ do all vector shared (alpha, data, drc, gcfpoint, ife, igcfail)
cmic$*       shared (ilat, ilon, ioe, iprs, ityp, iz, MaxObservPoints)
cmic$*       shared (mlev, nneigh, nobs, nreacc, rade, rdc, ztol)
cmic$*       private (cdif1, clat1, clat2, clond, dateline_case)
cmic$*       private (dist, fg1, fg2, gcflat, gpv, i, ifail, isum)
cmic$*       private (j, latmax, latmin, londif, lonmax, lonmin)
cmic$*       private (pointer, polar, srchiz, srchlat, srchlon)
cmic$*       private (srchprs, srchtyp, swts, swtsv, temp, tolerance)
cmic$*       private (trwt) autoscope
c  Define the search region.
c       print*,' box (latx,latn,lonx,lonn):', latmax,latmin,lonmax,
c    &           lonmin
c  Find the neighbours.
      nneigh = 0
      DO i = 1 , igcfail
        ifail = gcfpoint(i)
        polar = .FALSE.
        IF (data(ifail + ilat - 1) .GT. 84.9) THEN
          latmax = 90.
          latmin = 70.
          polar = .TRUE.
        ELSEIF (data(ifail + ilat - 1) .LT. -84.9) THEN
          latmax = -70.
          latmin = -90.
          polar = .TRUE.
        ELSE
          latmax = data(ifail + ilat - 1) + 10.
          latmin = data(ifail + ilat - 1) - 10.
        ENDIF
        IF (.NOT.polar) THEN
          gcflat = data(ifail + ilat - 1)
          londif = atan2(.3139E0,cos(gcflat * drc)) * rdc / 2.
          lonmax = data(ifail + ilon - 1) + londif
          lonmin = data(ifail + ilon - 1) - londif
c  Case where region overlaps the international date line.
          IF (lonmax .GT. 180.) THEN
            lonmax = lonmax - 360.
            dateline_case = .TRUE.
          ELSEIF (lonmin .LT. -180.) THEN
            lonmin = lonmin + 360.
            dateline_case = .TRUE.
          ELSE
            dateline_case = .FALSE.
          ENDIF
        ENDIF
        swts = 0.
        swtsv = 0.
        fg1(:nobs) = 0.
        DO ti1 = 1 , nobs
          fg2(ti1) = 0.
        ENDDO
        DO j = 1 , nobs
          srchlat = data(ilat + j - 1)
          srchlon = data(ilon + j - 1)
          srchiz = data(iz + j - 1)
          srchtyp = data(ityp + j - 1)
          srchprs = data(iprs + j - 1)
c   Compute distance between neighbour and suspect point.
          IF (srchlat .LE. latmax .AND. srchlat .GE. latmin .AND. 
     &          srchtyp .EQ. 1. .AND. (abs(srchprs - data(ifail + iprs -
     &           1)) .LT. 15.) .AND. (dateline_case .AND. srchlon 
     &          .GT. lonmin .OR. srchlon .LT. lonmax .OR. polar)) THEN
            clat1 = cos(srchlat * drc)
            clat2 = cos(data(ifail + ilat - 1) * drc)
            clond = cos((srchlon - data(ifail + ilon - 1)) * drc)
            cdif1 = cos((srchlat - data(ifail + ilat - 1)) * drc)
            temp = 1. - cdif1 + clat1 * clat2 * (1. - clond)
            dist = sqrt(2.0 * temp) * rade
c               print*,' neighbour:', srchlat, srchlon, srchiz
c  Compute the weight (alpha is the spatial convergence parameter).
c               weight = ( alpha - dist) / ( alpha + dist)
c               swts = swts + weight
c               swtsv = swtsv + weight*srchiz
            fg1(j) = (alpha - dist) / (alpha + dist)
            fg2(j) = ((alpha - dist) / (alpha + dist)) * srchiz
          ENDIF
        ENDDO
        DO isum = 1 , nobs
          swtsv = swtsv + fg2(isum)
        ENDDO
        DO isum = 1 , nobs
          swts = swts + fg1(isum)
        ENDDO
        lif = swts .GT. 0.
        IF (swts .GT. 0.) THEN
          gpv = swtsv / swts
          IF (abs(data(ifail + ilat - 1)) .LE. 30.) THEN
            trwt = exp(sin(abs(data(ifail + ilat - 1)) * drc)) / 1.65
          ELSE
            trwt = 1.
          ENDIF
        ENDIF
        IF (lif) THEN
          pointer = pindx2(data(ifail + iprs - 1))
          tolerance = ztol(int(pointer)) * trwt
c               data(ifail+ityp-1) = 1.2
          lif100 = ((data(ifail + iz - 1) - gpv) ** 2 .LT. (data(ifail 
     &        + ioe - 1) ** 2 + data(ifail + ife - 1) ** 2) * tolerance)
     &        .AND. data(ifail + ioe - 1) .LE. 900.
          IF (((data(ifail + iz - 1) - gpv) ** 2 .LT. (data(ifail + ioe
     &          - 1) ** 2 + data(ifail + ife - 1) ** 2) * tolerance) 
     &          .AND. data(ifail + ioe - 1) .LE. 900.) THEN
            IF (data(ifail + ityp - 1) .EQ. 0.1) THEN
              data(ifail + ityp - 1) = 1.2
            ENDIF
          ENDIF
        ENDIF
        IF (lif .AND. lif100 .AND. data(ifail + ityp - 1) .EQ. 0.3) THEN
          data(ifail + ityp - 1) = 1.3
        ENDIF
        IF (lif .AND. lif100) THEN
          nreacc = nreacc + 1
        ENDIF
      ENDDO
      WRITE( 6,'(1x,"DELHUV:",i8," height reports reaccepted by",       
     &" the buddy check.")' ) nreacc
      CALL qcreport(data,nobs,1,mlev,nymd,nhms)
      data(ityp:ityp + nobs - 1) = int(data(ityp:ityp + nobs - 1))
      CALL rptlly_v1(data,ikx,maxobservpoints * 10,data,ityp,maxobservpo
     &    ints * 10,nobs,4)
c  QC gross check / delta computation for wind data ----------------
      igcfail = 0
      DO ti1 = 1 , maxobservpoints
        fg3(ti1) = 0.
      ENDDO
c  Compute the forecast wind speed errors.
      DO 7016 i = 1 , im
        DO 7014 j = 1 , jnp
          DO 7012 k = 1 , mlev
            u_error(i,j,k) = sqrt(u_error(i,j,k) ** 2 + v_error(i,j,k) 
     &          ** 2)
 7012     CONTINUE
 7014   CONTINUE
 7016 CONTINUE
c  Interpolate first guess winds to locations of observations.
      CALL vterpv_v1(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,iprs,maxobservpoints * 10,fg1,1,
     &    maxobservpoints,nobs,umod)
      CALL vterpv_v2(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,iprs,maxobservpoints * 10,fg2,nobs,
     &    vmod)
      nreacc = 0
c  Interpolate the forecast errors to observation locations.
      CALL vterpv_v1(data,ilat,maxobservpoints * 10,data,ilon,
     &    maxobservpoints * 10,data,iprs,maxobservpoints * 10,data,ife,
     &    maxobservpoints * 10,nobs,u_error)
c  Interpolate observation error to observation location.
      CALL obertrp0_v1(data,ikx,maxobservpoints * 10,data,iprs,
     &    maxobservpoints * 10,data,ioe,maxobservpoints * 10,obserr,3,
     &    nobs)
cmic$ do all vector shared(data,fg1,fg2,fill,gcfpoint,ife,igcfail)
cmic$*       shared (ikx, ilat, ilon, ioe, iprs, ityp, iu, iv)
cmic$*       private ( lat,lon, kx)
cmic$*       shared (MaxObservPoints, nobs, psobs, suspect, fg3)
cmic$*       private (delta, fe, i, oe, prs, tolerance) autoscope
c  Now compute the delta if valid data exists.
c  Gross check the wind data.
      tolerance100(1:nobs) = 2.
      DO i = 1 , nobs
        IF (data(iu + i - 1) .LT. fill .AND. data(iv + i - 1) .LT. fill
     &        .AND. data(iprs + i - 1) .LE. psobs(i)) THEN
          data(ityp + i - 1) = data(ityp + i - 1) + 2.
        ELSE
          data(ityp + i - 1) = data(ityp + i - 1) + 0.04
        ENDIF
        lat = data(ilat + i - 1)
        lon = data(ilon + i - 1)
        prs = data(iprs + i - 1)
        kx = data(ikx + i - 1)
        IF (data(ityp + i - 1) .GE. 2.) THEN
          data(iu + i - 1) = data(iu + i - 1) - fg1(i)
          data(iv + i - 1) = data(iv + i - 1) - fg2(i)
        ENDIF
        delta = sqrt(data(iu + i - 1) ** 2 + data(iv + i - 1) ** 2)
        oe = data(ioe + i - 1)
        fe = data(ife + i - 1)
        IF (prs .LE. 500. .AND. prs .GT. 70.) THEN
          tolerance100(i) = 2.5
        ENDIF
        IF ((delta ** 2 .GT. (oe ** 2 + fe ** 2) * tolerance100(i)) 
     &        .AND. data(ityp + i - 1) .GE. 2. .AND. (.NOT.suspect(i,2))
     &        ) THEN
          data(ityp + i - 1) = data(ityp + i - 1) - 1.99
          fg3(i) = 1.
        ENDIF
        IF (suspect(i,2) .AND. data(ityp + i - 1) .GE. 2.) THEN
          data(ityp + i - 1) = data(ityp + i - 1) - 1.97
          fg3(i) = 1.
        ENDIF
      ENDDO
      TOLERANCE = TOLERANCE100(nobs)
      CALL wheneq(nobs,fg3,1,1.,gcfpoint,igcfail)
      WRITE( 6,* ) ' DELHUV: ',igcfail,' wind observations have',
     &' failed the gross check.'
c  Find neighbouring observations to those flagged by the gross check.
c  Perform an interpolation of the good quality neighbours (within a
c  1000km radius) to the location of a suspicious report.
cmic$ do all vector shared (alpha, data, drc, gcfpoint, ife, igcfail)
cmic$*       shared (ilat, ilon, ioe, iprs, ityp, iu, iv, MaxObservPoint
cmic$*       shared (nobs, nreacc, rade, rdc, srchiz, srchws)
cmic$*       private (cdif1, clat1, clat2, clond, dateline_case)
cmic$*       private (delta, dist, fg1, fg2, gcflat, gpv, i, ifail)
cmic$*       private (isum, j, latmax, latmin, londif, lonmax)
cmic$*       private (lonmin, polar, prs, srchlat, srchlon, srchprs)
cmic$*       private (srchtyp, swts, swtsv, temp, testval, tolerance)
cmic$*       autoscope
c  Define the search region.
c  Find the neighbours.
      DO i = 1 , igcfail
        ifail = gcfpoint(i)
        polar = .FALSE.
        IF (data(ifail + ilat - 1) .GT. 84.9) THEN
          latmax = 90.
          latmin = 70.
          polar = .TRUE.
        ELSEIF (data(ifail + ilat - 1) .LT. -84.9) THEN
          latmax = -70.
          latmin = -90.
          polar = .TRUE.
        ELSE
          latmax = data(ifail + ilat - 1) + 10.
          latmin = data(ifail + ilat - 1) - 10.
        ENDIF
        IF (.NOT.polar) THEN
          gcflat = data(ifail + ilat - 1)
          londif = atan2(.31390,cos(gcflat * drc)) * rdc / 2.
          lonmax = data(ifail + ilon - 1) + londif
          lonmin = data(ifail + ilon - 1) - londif
c  Case where box overlaps the international date line.
          IF (lonmax .GT. 180.) THEN
            lonmax = lonmax - 360.
            dateline_case = .TRUE.
          ELSEIF (lonmin .LT. -180.) THEN
            lonmin = lonmin + 360.
            dateline_case = .TRUE.
          ELSE
            dateline_case = .FALSE.
          ENDIF
        ENDIF
        fg1(:nobs) = 0.
        DO ti1 = 1 , nobs
          fg2(ti1) = 0.
        ENDDO
        swts = 0.
        swtsv = 0.
        DO j = 1 , nobs
          srchws = sqrt(data(iu + j - 1) ** 2 + data(iv + j - 1) ** 2)
        ENDDO
        DO j = 1 , nobs
          srchlat = data(ilat + j - 1)
          srchlon = data(ilon + j - 1)
          srchtyp = data(ityp + j - 1)
          srchprs = data(iprs + j - 1)
c   Compute distance between neighbour and suspect point.
          IF (srchlat .LE. latmax .AND. srchlat .GE. latmin .AND. 
     &          srchtyp .GE. 2. .AND. (abs(srchprs - data(ifail + iprs -
     &           1)) .LT. 15.) .AND. (dateline_case .AND. srchlon 
     &          .GT. lonmin .OR. srchlon .LT. lonmax .OR. polar)) THEN
            clat1 = cos(srchlat * drc)
            clat2 = cos(data(ifail + ilat - 1) * drc)
            clond = cos((srchlon - data(ifail + ilon - 1)) * drc)
            cdif1 = cos((srchlat - data(ifail + ilat - 1)) * drc)
            temp = 1. - cdif1 + clat1 * clat2 * (1. - clond)
            dist = sqrt(2.0 * temp) * rade
c               print*,' neighbour:', srchlat, srchlon, srchws
c  Compute the weight (alpha is the spatial convergence parameter).
            fg1(j) = (alpha - dist) / (alpha + dist)
            fg2(j) = ((alpha - dist) / (alpha + dist)) * srchiz
c               weight = ( alpha - dist) / ( alpha + dist)
c               swts = swts + weight
c               swtsv = swtsv + weight*srchws
          ENDIF
        ENDDO
        DO isum = 1 , nobs
          swtsv = swtsv + fg2(isum)
        ENDDO
        DO isum = 1 , nobs
          swts = swts + fg1(isum)
        ENDDO
        lif101 = swts .GT. 0.
        IF (swts .GT. 0.) THEN
          gpv = swtsv / swts
          prs = data(ifail + iprs - 1)
          IF (prs .LE. 70.) THEN
            tolerance = 2.
          ENDIF
        ENDIF
        IF (lif101 .AND. prs .GT. 500.) THEN
          tolerance = 2.
        ENDIF
        IF (lif101 .AND. prs .LE. 500. .AND. prs .GT. 70.) THEN
          tolerance = 2.5
        ENDIF
        IF (lif101) THEN
          delta = sqrt(data(ifail + iu - 1) ** 2 + data(ifail + iv - 1) 
     &        ** 2)
          lif102 = (delta - gpv) ** 2 .LT. (data(ifail + ioe - 1) ** 2 +
     &        data(ifail + ife - 1) ** 2) * tolerance .AND. data(ifail 
     &        + ioe - 1) .LE. 900.
          IF ((delta - gpv) ** 2 .LT. (data(ifail + ioe - 1) ** 2 + data
     &          (ifail + ife - 1) ** 2) * tolerance .AND. data(ifail 
     &          + ioe - 1) .LE. 900.) THEN
            testval = data(ifail + ityp - 1)
            IF ((int(testval * 100.) - int(testval * 10.) * 10.) .EQ. 3.
     &            ) THEN
              data(ifail + ityp - 1) = data(ifail + ityp - 1) + 2.
            ELSE
              data(ifail + ityp - 1) = data(ifail + ityp - 1) + 2.01
            ENDIF
          ENDIF
        ENDIF
        IF (lif101 .AND. lif102) THEN
          nreacc = nreacc + 1
        ENDIF
      ENDDO
      PRINT *, nreacc,' wind reports reaccepted by buddy check.'
      CALL qcreport(data,nobs,2,mlev,nymd,nhms)
      data(ityp:ityp + nobs - 1) = int(data(ityp:ityp + nobs - 1))
      WRITE( 91 ) nobs,nymd,nhms,'ZUV'
      DO i = 1 , nobs
        WRITE( 91 ) data(ilat + i - 1),data(ilon + i - 1),
     &      data(ityp + i - 1),data(iprs + i - 1),data(ikx + i - 1),
     &      data(iu + i - 1),data(iv + i - 1),data(iz + i - 1)
      ENDDO
      CALL rptlly_v1(data,ikx,maxobservpoints * 10,data,ityp,maxobservpo
     &    ints * 10,nobs,5)
c  Call offlvl to alter the type array to allow/disallow off-level
c  data.
      CALL offlvl(data,nobs)
      RETURN
      END
