************************************
c---------------------------------------------------------------------c
c         N A S A   /  D A O     A N A L Y S I S   G R O U P          c
c    Optimum Interpolation Analysis System   Version 1.2  Fall, 1991  c
c                                                                     c
c.....................................................................c
C                                                                     .
C  Program:    QCSLP                                                  .
C  Language:   Fortran 77 (with Fortran 90 enhancements)              .
C  Computer:   Cray YMP 864                                           .
C                                                                     .
C  Programmer: M. Seablom                                             .
C  Date:       December 16, 1991
C                                                                     .
C  Program QCSLP  computes sea-level pressure and wind deltas, and    .
C  performs version 2.1 quality control.                              .
C                                                                     .
C**********************************************************************
C  Data parallel implementation (CM Fortran) : M. Makivic (NPAC)      *
C  Date:       4/94                                                   *
C**********************************************************************
C**********************************************************************
      SUBROUTINE qcslp(phis,slp,slu,slv,data_mm,nobs,t,ze,zelat,zelon,
     &    maxnze,fgslp,sspct_mm,fcerr)
      include 'params.h'
      include 'processors.h'
      

      
      real data_mm(nobs,NRECORDS)
      integer sspct_mm(nobs,2)

      REAL phis (im, jnp), slp (im, jnp), slu (im, jnp), slv (im, jnp), 
     &zelat (20000), zelon (20000), t (im, 
     &jnp, mlev), fgslp (im, jnp), pm (nobs), work1 (nobs), 
     &work2 (nobs), 
     &fcerr (im, jnp), work3 (nobs), 
     &alpha, ze (20000), phis1 (im,jnp)
      real arrtest(30000)
      INTEGER gc_failures, mgcfpoint (30000), reaccepts
      COMMON / sfcdrag / drag (imr,jmr+1), iwater(13104)
      COMMON / dataptrs / ikx, ilat, ilon, iprs, ityp, islp, ifcerr
      COMMON / cortabs / hcorfe (1000, 3, 20, 2)
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / gridxx / im, jnp, mlev
      COMMON / obs / maxdex, nfgge, lun, hh (13, 13, 20),
     &vcoroe (18, 18, 4), obserr (mxqobs, 4, 19), ww (13, 13, 20),
     &qobs (mxqobs, 5), offlev (mxqobs)
      common / obschar / dataname (5), obsname(mxqobs)
      real obserr1(mxqobs),obserr2(mxqobs),obserr3(mxqobs)
      real obserr4(mxqobs)

      REAL indx, lat, lon, xdat1, xdat2, diff, type
      LOGICAL exit, onestage, obs, 
     &offlev
      CHARACTER * 60 message
      CHARACTER * 8 dataname, obsname
      real,dimension(1:nobs):: clat1,clat2,clond,cdif1
      real,dimension(1:nobs):: distance,weight,temp
      real,dimension(nobs) :: zze,zzelat,zzelon
      integer,dimension(nobs) :: nz,nze
c This may be a bug in seq code, latitude should be real
      real,dimension(1:nobs):: latitude
      real,dimension(1:nobs):: coriolis_sfc_coeff
      real,dimension(1:nobs):: drag_sfc_coeff
c This may be a bug in seq code, longitude should be real
      real,dimension(1:nobs):: longitude
      real,dimension(1:nobs):: datslp
      real,dimension(1:nobs):: kx
      real,dimension(1:nobs):: datfcerr
      real,dimension(1:nobs):: dattypw
      real,dimension(1:nobs):: search_slw,search_lat,search_lon
      real,dimension(1:nobs):: search_typ
      logical,dimension(1:nobs):: warn
      real,dimension(1:nobs):: delta
      integer,dimension(1:nobs):: intkx
      real,dimension(1:nobs):: observation_error 
      real,dimension(1:nobs):: forecast_error
      logical,dimension(1:nobs):: mask1,mask2,mask3

      real,allocatable :: gcfslu(:)
      real,allocatable :: gcfslv(:)
      real :: delta_slwfail, testval
      real,allocatable :: gcflat(:)
      real,allocatable :: gcflon(:)
      real,allocatable :: gcftypw(:)
      real,allocatable :: gcfslp(:)
      integer,allocatable :: gcfikx(:)
      real,allocatable :: gcffcerr(:)
      integer,allocatable :: gcfpoint(:)
      logical :: polar,dateline_case
      real:: latmax,latmin,lonmax,lonmin
c   ******** Is this an error in the original OI code? ***********
c   ******** Shouldn't londif be declared as real?     ***********
      real :: londif
      real :: sumofweights
      real :: value_sow
      real :: analyzed_value
      integer,allocatable :: gcfmask1(:),gcfmask2(:),gcfmask3(:)

!hpf$ distribute drag(*,block) onto PROC_ARRAY_1D
!hpf$ distribute iwater(block) onto PROC_ARRAY_1D
!hpf$ distribute data_mm(block,*) onto PROC_ARRAY_1D
!hpf$ distribute sspct_mm(block,*) onto PROC_ARRAY_1D

!hpf$ distribute fcerr(*,block) onto PROC_ARRAY_1D
!hpf$ distribute phis(*,block) onto PROC_ARRAY_1D
!hpf$ distribute fgslp(*,block) onto PROC_ARRAY_1D
!hpf$ distribute slp(*,block) onto PROC_ARRAY_1D
!hpf$ distribute slu(*,block) onto PROC_ARRAY_1D
!hpf$ distribute slv(*,block) onto PROC_ARRAY_1D
!hpf$ distribute t(*,block,*) onto PROC_ARRAY_1D
!hpf$ distribute ze(block) onto PROC_ARRAY_1D
!hpf$ distribute zelat(block) onto PROC_ARRAY_1D
!hpf$ distribute zelon(block) onto PROC_ARRAY_1D


!hpf$ distribute mgcfpoint(block) onto PROC_ARRAY_1D

!hpf$ distribute clat1(block) onto PROC_ARRAY_1D
!hpf$ distribute clat2(block) onto PROC_ARRAY_1D
!hpf$ distribute clond(block) onto PROC_ARRAY_1D
!hpf$ distribute cdif1(block) onto PROC_ARRAY_1D
!hpf$ distribute distance(block) onto PROC_ARRAY_1D
!hpf$ distribute weight(block) onto PROC_ARRAY_1D
!hpf$ distribute temp(block) onto PROC_ARRAY_1D
!hpf$ distribute work1(block) onto PROC_ARRAY_1D
!hpf$ distribute work2(block) onto PROC_ARRAY_1D
!hpf$ distribute work3(block) onto PROC_ARRAY_1D
!hpf$ align with clat1:: zze,zzelat,zzelon,nz,nze
!hpf$ distribute latitude(block) onto PROC_ARRAY_1D
!hpf$ distribute coriolis_sfc_coeff(block) onto PROC_ARRAY_1D
!hpf$ distribute drag_sfc_coeff(block) onto PROC_ARRAY_1D
!hpf$ distribute longitude(block) onto PROC_ARRAY_1D
!hpf$ distribute dattypw(block) onto PROC_ARRAY_1D
!hpf$ distribute search_slw(block) onto PROC_ARRAY_1D
!hpf$ distribute warn(block) onto PROC_ARRAY_1D
!hpf$ distribute delta(block) onto PROC_ARRAY_1D
!hpf$ distribute intkx(block) onto PROC_ARRAY_1D
!hpf$ distribute observation_error(block) onto PROC_ARRAY_1D
!hpf$ distribute mask1(block) onto PROC_ARRAY_1D
!hpf$ distribute mask2(block) onto PROC_ARRAY_1D

!hpf$ distribute gcftypw(block) onto PROC_ARRAY_1D
!hpf$ distribute gcfslp(block) onto PROC_ARRAY_1D
!hpf$ distribute gcffcerr(block) onto PROC_ARRAY_1D
!hpf$ distribute gcfpoint(block) onto PROC_ARRAY_1D


      DATA slp_tolerance / 4. /, slw_tolerance / 2. /, 
     & rad_e / 6.37e6 /,
     & alpha / 5.0625e10 /, fill / 1.e10 /

      message(1:60) = ' '
      work1( 1:nobs ) = 0.
      work2( 1:nobs ) = 0.
      work3( 1:nobs ) = 0.
      pi = 4. * atan(1.)
      drc = pi / 180.
      rdc = 180. / pi
C  Define additional pointers to the DATA array.
      CALL getconi('NYMD    ',nymd)
      CALL getconi('NHMS    ',nhms)
      CALL getconh('RHOBAR  ',rhobar)
      CALL getconh('RADEI   ',radei)
      CALL getconh('OMEGA2  ',omega2)
      itypw=ind_type
      data_mm(:,itypw) = 0.
      latitude=data_mm(1:nobs,ind_lat)
      longitude=data_mm(1:nobs,ind_lon)
      datslp=data_mm(1:nobs,ind_slp)
      kx=data_mm(1:nobs,ind_kx)
      write(6,*) 'ready to interpolate'
C  Obtain first-guess pressure values at the locations of the
C  observations.
      CALL terpv(latitude,longitude,work1,nobs,slp)
C  Interpolate forecast errors to observation locations.
      CALL terpv(latitude,longitude,datfcerr,nobs,fcerr)
      data_mm(1:nobs,ind_fcer)=datfcerr
C  Establish loop to compute deltas and perform the gross check of
C  the sea-level pressure data.
      gc_failures = 0
      where (datslp .LT. fill)
        dattypw = 1.
      elsewhere
        dattypw = 0.4
      endwhere
      where (latitude .GT. 90. .OR. latitude .LT. -90. .OR. longitude 
     &        .GT. 180. .OR. longitude .LT. -180.)
        dattypw= 0.4
        warn=.true.
      endwhere
      IF(any(warn)) THEN
         message =  'QCSLP:  Impossible data location; rejected'
     &               // 'point.'
        CALL flash(3,message,1)
      ENDIF
      where (dattypw .EQ. 1.)
        datslp= datslp- work1
      endwhere
      delta = datslp
      data_mm(1:nobs,ind_slp)=datslp
      intkx=int(kx)
      obserr1(1:mxqobs)=obserr(1:mxqobs,1,19)
      obserr2(1:mxqobs)=obserr(1:mxqobs,2,19)
      obserr3(1:mxqobs)=obserr(1:mxqobs,3,19)
      obserr4(1:mxqobs)=obserr(1:mxqobs,4,19)
      observation_error = obserr(intkx,1,19)
      forecast_error = datfcerr
      where (delta ** 2 .GT. (observation_error ** 2 + forecast_error 
     &        ** 2) * slp_tolerance .AND. dattypw .NE. .4)
          dattypw = 0.1
          work3 = 1.
      endwhere
C  Suspect observations automatically fail the gross check.
      where (sspct_mm(:,1).eq.1 .AND. datslp .LT. fill)
          dattypw = 0.3
          work3 = 1.
      endwhere
      CALL wheneq(nobs,work3,1,1.,mgcfpoint,gc_failures)
      IF (gc_failures .GT. 30000) THEN
        message = 'QCSLP:  Excessive gross check failures.'
        CALL flash(5,message,1)
      ENDIF
      WRITE( 6,* ) ' QCSLP:  ',gc_failures,' sea-level pressure',
     &' observations have failed the gross check.'
      reaccepts = 0
C  Analyze data surrounding an observation which failed the gross
C  check. Search for observations within 1000km of the suspect point.

      allocate(gcflat(gc_failures))
      allocate(gcflon(gc_failures))
      allocate(gcftypw(gc_failures))
      allocate(gcfslp(gc_failures))
      allocate(gcfikx(gc_failures))
      allocate(gcffcerr(gc_failures))
      allocate(gcfpoint(gc_failures))
      gcfpoint=mgcfpoint(1:gc_failures)


      gcflat = latitude(gcfpoint)
c      print *, 'gcflat'
c      print *, gcflat
      gcflon = longitude(gcfpoint)
c      print *, 'gcflon'
c      print *, gcflon
      gcftypw=dattypw(gcfpoint)
c      print *, 'gcftypw'
c      print *, gcftypw
      gcfslp=datslp(gcfpoint)
c      print *, 'gcfslp'
c      print *, gcfslp
      gcfikx=int(kx(gcfpoint))
c      print *, 'gcfikx'
c      print *, gcfikx
      gcffcerr=datfcerr(gcfpoint)
c      print *, 'gcffcerr'
c      print *, gcffcerr

      do i = 1,gc_failures   !  ***  Begin SLP buddy check loop ****
 
        polar = .false.
 
        if ( gcflat(i) .gt. 84.9 ) then
            polar = .true.
            latmax = 90.
            latmin = 70.
        else if ( gcflat(i) .lt. -84.9 ) then
            polar = .true.
            latmax = -70.
            latmin = -90.
        else
            latmax = gcflat(i)  + 10.
            latmin = gcflat(i)  - 10.
        end if
 
        if ( .not. polar ) then
 
            londif = atan2( .3139, cos( gcflat(i)*drc ) ) * rdc * .5
c        print *, i,polar,gcflat(i),gcflon(i),drc,rdc,
c     c    cos( gcflat(i)*drc ),atan2( .3139, cos( gcflat(i)*drc ) ),
c     c    londif
            lonmax = gcflon(i) + londif
            lonmin = gcflon(i) - londif
 
!  Take care of points near the international date line.
 
           if ( lonmax .gt. 180. ) then
              lonmax = lonmax - 360.
              dateline_case = .true.
           else if ( lonmin .lt. -180. ) then
              lonmin = lonmin + 360.
              dateline_case = .true.
           else
              dateline_case = .false.
           end if
 
        end if
 
!  Find high quality neighbouring observations.
 
        neighbours = 0
        sumofweights = 0.
        value_sow   = 0.
 
        work1( 1:nobs ) = 0.
        work3( 1:nobs ) = 0.
 
c      print *,i,polar,latmax,latmin,lonmin,lonmax,londif,dateline_case
        mask1 = latitude.le. latmax .and. latitude .ge.
     &         latmin .and. dattypw .eq. 1. 
        mask2=(( .not.(dateline_case).and. longitude.gt.lonmin
     &                                 .and. longitude.lt.lonmax)
     &           .or.  (dateline_case  .and. longitude.gt.lonmin
     &                                 .or.  longitude.lt.lonmax)
     &           .or.   polar )
        mask1=mask1.and.mask2
        where(mask1)
                  clat1 = cos ( latitude * drc )
                  clat2 = cos ( gcflat(i)*drc )
                  clond = cos ((longitude - gcflon(i))*drc)
                  cdif1 = cos ((latitude - gcflat(i))*drc)
                  temp  = 1.-cdif1+ clat1+clat2*(1.-clond)
                  distance =  sqrt ( 2.*temp ) * rad_e
 
                  weight = ( alpha - distance ) / ( alpha + distance )
                  work1 = weight
                  work3 = weight * datslp
         else where
                  work1=0.
                  work3=0.
         end where
 
         sumofweights = sum(work1)
         value_sow    = sum(work3)
 
        if ( sumofweights .gt. 0. ) then

          analyzed_value = value_sow / sumofweights
c          print *, i,gcfslp(i),analyzed_value
c          print *, i,obserr1( gcfikx(i))
c          print *, i,gcffcerr(i),slp_tolerance

 
          if (  (gcfslp(i) - analyzed_value)**2 .lt.
     &          (obserr1( gcfikx(i) )**2 +
     &           gcffcerr(i)**2) * slp_tolerance ) then
 
 
               if ( gcftypw(i) .eq. 0.1 ) then
                 gcftypw(i) = 1.2
                 reaccepts = reaccepts + 1
               else if ( gcftypw(i) .eq. 0.3 ) then
                 gcftypw(i) = 1.3
                 reaccepts = reaccepts + 1
               end if
 
          end if
c          print *, i, gcftypw(i), reaccepts


 
        end if
        dattypw(gcfpoint)=gcftypw
 
2     end do                   ! *** End SLP buddy check loop ********
 
      write(6,'(//," QCSLP: ",i8," pressure reports reaccepted by",
     &      " the buddy check.")' ) reaccepts
      data_mm(:,itypw) = dattypw
      data_mm(1:nobs,ind_slp)=datslp

      data_mm(gcfpoint,ind_type) = gcftypw
 
!  Report quality control results.
 
      call sqcrpt ( data_mm, nobs, 1, nymd, nhms )
      data_mm(:,ind_type)=int(data_mm(:,ind_type))
 
      CALL rptlly(data_mm(:,ind_kx),data_mm(:,ind_type),nobs,1)
 
      deallocate(gcflat)
      deallocate(gcflon)
      deallocate(gcftypw)
      deallocate(gcfslp)
      deallocate(gcfikx)
      deallocate(gcffcerr)
      deallocate(gcfpoint)

C  Perform a similar quality control for the winds.
C  Obtain first-guess sea-level wind values at the locations of the
C  observations.

      open(unit=33,file='data2.out',status='unknown')
      CALL terpv(latitude,longitude,work1,nobs,slu)
      CALL terpv(latitude,longitude,work2,nobs,slv)

      islu=ind_uwnd
      islv=ind_vwnd
      itypw=ind_type
      ifcerr=ind_fcer

      gc_failures = 0

      WHERE(data_mm(:,islu).LT.fill.AND.data_mm(:,islv).LT. fill) 
          data_mm(:,itypw) = data_mm(:,itypw) + 2.
      ELSEWHERE
          data_mm(:,itypw) = data_mm(:,itypw) + 0.04
      ENDWHERE
C  Establish loop to compute deltas and perform the gross check.
      work3 = 0.
      WHERE (latitude .GT. 90. .OR. latitude .LT. -90. .OR. longitude 
     &        .GT. 180. .OR. longitude .LT. -180.) 
          data_mm(:,itypw) = 0.4
          warn=.TRUE.
      ENDWHERE
      IF(any(warn)) THEN
          message = 'QCSLP : Impossible data location; rejected' 
     &        // ' point.'
          CALL flash(3,message,1)
      ENDIF
      delta=0.0
      WHERE(data_mm(:,itypw) .GE. 2.)
          data_mm(:,islu) = data_mm(:,islu) - work1
          data_mm(:,islv) = data_mm(:,islv) - work2
          delta=sqrt(data_mm(:,islu)**2 + data_mm(:,islv)**2)
      ENDWHERE
      observation_error = obserr3(intkx)
C  Calculate the wind forecast error.
      WHERE(abs(latitude) .LT. 1.E-5) 
          latitude = 1.E-5
      ENDWHERE
      print *, rhobar, omega2, drc
      coriolis_sfc_coeff = .1 * (rhobar * omega2 * sin(latitude * drc)
     &      ) / ((rhobar * omega2 * sin(latitude * drc)) ** 2 + 1.E-10)
      drag_sfc_coeff = 1.E-6 / ((rhobar * omega2 * sin(latitude * drc)
     &      ) ** 2 + 1.E-10)
      forecast_error = sqrt(coriolis_sfc_coeff ** 2 + drag_sfc_coeff 
     &      ** 2) * radei * data_mm(:,ifcerr) 
     &      * sqrt(-hcorfe(1,3,1,1))
      WHERE(delta ** 2 .GT. ((observation_error ** 2 + forecast_error 
     &        ** 2) * slw_tolerance) .AND. data_mm(:,itypw) .GE. 2. 
     &        .AND. (sspct_mm(:,2).eq.0))
          data_mm(:,itypw) = data_mm(:,itypw) - 1.99
          work3 = 1.
      ELSEWHERE
          work3 = 0.
      ENDWHERE
C  Suspect observations automatically fail the gross check.
      WHERE(sspct_mm(:,2).eq.1 .AND. data_mm(:,islu) .LT. fill)
          data_mm(:,itypw) = data_mm(:,itypw) - 1.97
          work3 = 1.
      ENDWHERE

      CALL wheneq(nobs,work3,1,1.,mgcfpoint,gc_failures)
      IF (gc_failures .GT. 30000) THEN
        message = 'QCSLP:  Excessive gross check failures.'
        CALL flash(5,message,1)
      ENDIF
      WRITE( 6,* ) ' QCSLP:  ',gc_failures,' sea-level wind',' observati
     &ons have failed the gross check.'
      reaccepts = 0
C  Analyze data surrounding an observation which failed the gross
C  check. Search for observations within 1000km of the suspect point.
C  Find high quality neighbouring observations.


      allocate(gcflat(gc_failures))
      allocate(gcflon(gc_failures))
      allocate(gcftypw(gc_failures))
      allocate(gcfslu(gc_failures))
      allocate(gcfslv(gc_failures))
      allocate(gcfikx(gc_failures))
      allocate(gcffcerr(gc_failures))
      allocate(gcfpoint(gc_failures))
      gcfpoint=mgcfpoint(1:gc_failures)

      latitude=data_mm(:,ind_lat)
      longitude=data_mm(:,ind_lon)
      search_lat=latitude
      search_lon=longitude
      dattypw = data_mm(:, itypw )
      gcflat = latitude(gcfpoint)
      gcflon = longitude(gcfpoint)
      gcftypw=dattypw(gcfpoint)
      gcfikx=int(kx(gcfpoint))
      gcffcerr=datfcerr(gcfpoint)
      gcfslu=data_mm(gcfpoint,islu)
      gcfslv=data_mm(gcfpoint,islv)

      search_slw = sqrt(data_mm(:, islu)**2 + data_mm(:,islv)**2 )

      do 3 i = 1,gc_failures   !  ***  Begin SLW buddy check loop ****
 
        polar = .false.
        ifail = gcfpoint( i )

 
        if ( gcflat(i) .gt. 84.9 ) then
            polar = .true.
            latmax = 90.
            latmin = 70.
        else if ( gcflat(i) .lt. -84.9 ) then
            polar = .true.
            latmax = -70.
            latmin = -90.
        else
            latmax = gcflat(i)  + 10.
            latmin = gcflat(i)  - 10.
        end if
 
        if ( .not. polar ) then
 
            londif = atan2( .3139, cos( gcflat(i)*drc ) ) * rdc * .5
            lonmax = gcflon(i) + londif
            lonmin = gcflon(i) - londif
 
!  Take care of points near the international date line.
 
           if ( lonmax .gt. 180. ) then
              lonmax = lonmax - 360.
              dateline_case = .true.
           else if ( lonmin .lt. -180. ) then
              lonmin = lonmin + 360.
              dateline_case = .true.
           else
              dateline_case = .false.
           end if

 
        end if
 
!  Find high quality neighbouring observations.
 
        neighbours = 0
        sumofweights = 0.
        value_sow   = 0.
 
        work1( 1:nobs ) = 0.
        work3( 1:nobs ) = 0.

        search_typ=data_mm(:,itypw)
        mask1 = latitude .le. latmax .and. latitude .ge.
     &         latmin .and. search_typ .ge. 2.
        mask2 = (( .not.(dateline_case).and. longitude.gt.lonmin
     &                                 .and. longitude.lt.lonmax)
     &           .or.  (dateline_case  .and. longitude.gt.lonmin
     &                                 .or.  longitude.lt.lonmax)
     &           .or.   polar )
        mask1=mask1.and.mask2

        where(mask1)
                  clat1 = cos ( latitude * drc )
                  clat2 = cos ( gcflat(i)*drc )
                  clond = cos ((longitude - gcflon(i))*drc)
                  cdif1 = cos ((latitude - gcflat(i))*drc)
                  temp  = 1.-cdif1+ clat1+clat2*(1.-clond)
                  distance =  sqrt ( 2.*temp ) * rad_e

                  weight = ( alpha - distance ) / ( alpha + distance )
                  work1 = weight
                  work3 = weight * search_slw
         else where
                  work1=0.
                  work3=0.
                  clat1=0.
                  clat2=0.
                  clond=0.
                  cdif1=0.
                  temp=0.
                  distance=0.
                  weight=0.

         end where

c        if(i.eq.88) then
c              write(33,*) work3
c          print *, latmax,latmin,lonmax,lonmin,dateline_case,polar
c        end if

        sumofweights = sum(work1)
        value_sow    = sum(work3)

 
        if ( sumofweights .gt. 0. ) then
 
          analyzed_value = value_sow / sumofweights
          delta_slwfail = sqrt( data_mm(ifail,islu)**2 +
     &                       data_mm(ifail,islv)**2 )
 
          if (  (delta_slwfail - analyzed_value)**2 .lt.
     &          (obserr( int(data_mm( ifail,ind_kx )),3,19 )**2 +
     &           data_mm(ifail,ifcerr)**2) * slw_tolerance ) then
 
               testval = data_mm( ifail,itypw )
               if ((int(testval*100.)-int(testval*10.)*10.).eq.3.)then
                 data_mm( ifail,itypw ) = data_mm( ifail,itypw ) + 2.
                 reaccepts = reaccepts + 1
               else
                 data_mm( ifail,itypw ) = data_mm(ifail,itypw) + 2.01
                 reaccepts = reaccepts + 1
               end if
 
          end if
 
        end if
 
3     end do                   ! *** End SLW buddy check loop ********
 
      print*,reaccepts,' wind reports reaccepted by buddy check.'

      call sqcrpt ( data_mm, nobs, 2, nymd, nhms )
      data_mm(:,itypw) = int(data_mm(:,itypw))

      work1 = 0.
      work2 = 0.
C  Create 1000mb heights based upon lowest level first guess temperature
C  and surface pressure.
      r_d = 287.
      g = 9.81
C  Interpolate temperatures to locations of observations.
      CALL terpv(latitude,longitude,work1,nobs,t(:,:,mlev))
      CALL terpv(latitude,longitude,work2,nobs,fgslp)

      WHERE (data_mm(:,ind_kx).GE.3. .AND.data_mm(:,ind_kx).LE.6.
     &        .AND. data_mm(:,ind_slp) .LT. 1.E10)
          nz =  1
          zze = (r_d * work1 / g) * ((1000.
     &        / work2) * (data_mm(:,ind_slp)
     &        / work2))
          zzelat = latitude
          zzelon = longitude
          warn=.true.
      ELSEWHERE
          nz=0
          warn=.false.
      ENDWHERE
c      forall(i=1:nobs) nze(i)=sum(nz(1:i),dim=1)
c      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000)
c     &             ze(nze(i))=zze(i)
c      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000)
c     &             zelat(nze(i))=zzelat(i)
c      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000)
c     &             zelon(nze(i))=zzelon(i)
c      maxnze=maxval(nze)
c      IF (maxnze .GE. 20000) THEN
c          message = ' Excessive GEN 1000MB heights. First 20000 used.'
c          CALL flash(3,message,1)
c          exit = .TRUE.
c      ENDIF
cC  Write out to DEL dataset.
c      WRITE( 91 ) nobs,nymd,nhms,'SLP'
c      DO i = 1 , nobs
c        WRITE( 91 ) data_mm(i,ind_lat),data_mm(i,ind_lon),
c     &  data_mm(i,ind_type),data_mm(i,ind_slp),data_mm(i,ind_kx),
c     &  data_mm(i,islu),data_mm(i,islv),data_mm(i,ind_slp)
c      ENDDO
cC  Tally all observations.
      CALL rptlly(data_mm(:,ind_kx),data_mm(:,ind_type),nobs,1)
      CALL rptlly(data_mm(:,ind_kx),data_mm(:,ind_type),nobs,2)

      deallocate(gcflat)
      deallocate(gcflon)
      deallocate(gcftypw)
      deallocate(gcfslu)
      deallocate(gcfslv)
      deallocate(gcfikx)
      deallocate(gcffcerr)
      deallocate(gcfpoint)

      RETURN
      END
