C************************************
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'

      
      real data_mm(nobs,NRECORDS)
!hpf$ distribute data_mm(block,block)
      integer sspct_mm(nobs,2)
!hpf$ distribute sspct_mm(block,block)

      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)
      INTEGER gc_failures, mgcfpoint (30000), reaccepts
      COMMON / sfcdrag / drag (imr,jmr+1), iwaterflag (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
!hpf$ distribute drag(block,block)
!hpf$ distribute fcerr(block,block)
!hpf$ distribute fgslp(block,block)
!hpf$ distribute mgcfpoint(block)
!hpf$ distribute hcorfe(block,block,block,block)
!hpf$ distribute hh(block,block,block)
!hpf$ distribute iwaterflag(block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute obserr1(block)
!hpf$ distribute obserr2(block)
!hpf$ distribute obserr3(block)
!hpf$ distribute obserr4(block)
!hpf$ distribute offlev(block)
!hpf$ distribute phis(block,block)
!hpf$ distribute phis1(block,block)
!hpf$ distribute pm(block)
!hpf$ distribute qobs(block,block)
!hpf$ distribute slp(block,block)
!hpf$ distribute slu(block,block)
!hpf$ distribute slv(block,block)
!hpf$ distribute t(block,block,block)
!hpf$ distribute vcoroe(block,block,block)
!hpf$ distribute work1(block)
!hpf$ distribute work2(block)
!hpf$ distribute work3(block)
!hpf$ distribute ww(block,block,block)
!hpf$ distribute ze(block)
!hpf$ distribute zelat(block)
!hpf$ distribute zelon(block)
      real,dimension(nobs) :: zze,zzelat,zzelon
      integer,dimension(nobs) :: nz,nze
!hpf$ distribute nze(block)
!hpf$ distribute nz(block)
!hpf$ distribute zze(block)
!hpf$ distribute zzelat(block)
!hpf$ distribute zzelon(block)
      real,dimension(1:nobs):: latitude
!hpf$ distribute latitude(block)
      real,dimension(1:nobs):: coriolis_sfc_coeff
!hpf$ distribute coriolis_sfc_coeff(block)
      real,dimension(1:nobs):: drag_sfc_coeff
!hpf$ distribute drag_sfc_coeff(block)
      real,dimension(1:nobs):: longitude
!hpf$ distribute longitude(block)
      real,dimension(1:nobs):: datslp
!hpf$ distribute datslp(block)
      real,dimension(1:nobs):: kx
!hpf$ distribute kx(block)
      real,dimension(1:nobs):: datfcerr
!hpf$ distribute datfcerr(block)
      real,dimension(1:nobs):: dattypw
!hpf$ distribute dattypw(block)
      real,dimension(1:nobs):: datslu
!hpf$ distribute datslu(block)
      real,dimension(1:nobs):: datslv
!hpf$ distribute datslv(block)
      real,dimension(1:nobs):: search_slw
!hpf$ distribute search_slw(block)
      real,dimension(1:nobs):: search_typ
!hpf$ distribute search_typ(block)
      logical,dimension(1:nobs):: warn
!hpf$ distribute warn(block)
      real,dimension(1:nobs):: delta
!hpf$ distribute delta(block)
      integer,dimension(1:nobs):: intkx
!hpf$ distribute intkx(block)
      real,dimension(1:nobs):: observation_error 
!hpf$ distribute observation_error(block)
      real,dimension(1:nobs):: forecast_error
!hpf$ distribute forecast_error(block)
      real,dimension(1:nobs):: clat1
!hpf$ distribute clat1(block)
      real,dimension(1:nobs):: clat2
!hpf$ distribute clat2(block)
      real,dimension(1:nobs):: clond
!hpf$ distribute clond(block)
      real,dimension(1:nobs):: cdif1
!hpf$ distribute cdif1(block)
      real,dimension(1:nobs):: temp
!hpf$ distribute temp(block)
      real,dimension(1:nobs):: distance
!hpf$ distribute distance(block)
      real,dimension(1:nobs):: weight
!hpf$ distribute weight(block)
      logical,dimension(1:nobs):: mask1,mask2,mask3
!hpf$ distribute mask1(block)
!hpf$ distribute mask2(block)
!hpf$ distribute mask3(block)

      real,allocatable :: gcfslu(:)
!hpf$ distribute gcfslu(block)
      real,allocatable :: gcfslv(:)
!hpf$ distribute gcfslv(block)
      real,allocatable :: testval(:)
!hpf$ distribute testval(block)
      real,allocatable :: delta_slwfail(:)
!hpf$ distribute delta_slwfail(block)
      real,allocatable :: gcflat(:)
!hpf$ distribute gcflat(block)
      real,allocatable :: gcflon(:)
!hpf$ distribute gcflon(block)
      real,allocatable :: gcftypw(:)
!hpf$ distribute gcftypw(block)
      real,allocatable :: gcfslp(:)
!hpf$ distribute gcfslp(block)
      integer,allocatable :: gcfikx(:)
!hpf$ distribute gcfikx(block)
      real,allocatable :: gcffcerr(:)
!hpf$ distribute gcffcerr(block)
      integer,allocatable :: gcfpoint(:)
!hpf$ distribute gcfpoint(block)
      integer,allocatable :: polar(:)
!hpf$ distribute polar(block)
      real,allocatable :: latmax(:)
!hpf$ distribute latmax(block)
      real,allocatable :: latmin(:)
!hpf$ distribute latmin(block)
      real,allocatable :: lonmax(:)
!hpf$ distribute lonmax(block)
      real,allocatable :: lonmin(:)
!hpf$ distribute lonmin(block)
      real,allocatable :: londif(:)
!hpf$ distribute londif(block)
      real,allocatable :: sumofweights(:)
!hpf$ distribute sumofweights(block)
      real,allocatable :: value_sow(:)
!hpf$ distribute value_sow(block)
      real,allocatable :: analyzed_value(:)
!hpf$ distribute analyzed_value(block)
      integer,allocatable :: dateline_case(:)
!hpf$ distribute dateline_case(block)
      integer,allocatable :: reacc(:)
!hpf$ distribute reacc(block)
      integer,allocatable :: gcfmask1(:),gcfmask2(:),gcfmask3(:)
!hpf$ distribute gcfmask1(block)
!hpf$ distribute gcfmask2(block)
!hpf$ distribute gcfmask3(block)


      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)
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)

      allocate(polar(gc_failures))
      allocate(latmax(gc_failures))
      allocate(latmin(gc_failures))
      allocate(lonmax(gc_failures))
      allocate(lonmin(gc_failures))
      allocate(londif(gc_failures))
      allocate(gcfmask1(gc_failures))
      allocate(gcfmask2(gc_failures))
      allocate(gcfmask3(gc_failures))
      allocate(dateline_case(gc_failures))
      allocate(sumofweights(gc_failures))
      allocate(analyzed_value(gc_failures))
      allocate(value_sow(gc_failures))
      allocate(reacc(gc_failures))

      print *, 'dattypw'
      print *, dattypw

      gcflat = latitude(gcfpoint)
      gcflon = longitude(gcfpoint)
      gcftypw=dattypw(gcfpoint)
      gcfslp=datslp(gcfpoint)
      gcfikx=int(kx(gcfpoint))
      gcffcerr=datfcerr(gcfpoint)

      print *, gcfpoint
      print *, gcfikx

      polar = 0
      WHERE(gcflat .GT. 84.9)
          polar = 1
          latmax = 90.
          latmin = 70.
      ENDWHERE
      WHERE (gcflat .LT. -84.9)
          polar = 1
          latmax = -70.
          latmin = -90.
      ENDWHERE
      WHERE((.not.(gcflat.GT.84.9)).and.(.not.(gcflat.LT.-84.9)))
          latmax = gcflat + 10.
          latmin = gcflat - 10.
      ENDWHERE
      WHERE(polar.eq.0) 
          londif = atan2(.3139,cos(gcflat * drc)) * rdc * .5
          lonmax = gcflon + londif
          lonmin = gcflon - londif
      end where
C  Take care of points near the international date line.
      gcfmask1=0
      gcfmask2=0
      gcfmask3=0
      where(polar.eq.0)
           gcfmask1=1
      end where
      where(lonmax .GT. 180.)
        gcfmask2=1
      end where
      where(lonmin .LT. -180.)
        gcfmask3=1
      endwhere
      WHERE (gcfmask1.eq.1.and.gcfmask2.eq.1)
            lonmax = lonmax - 360.
            dateline_case = 1
      end where
      WHERE (gcfmask1.eq.1.and.gcfmask2.eq.0.and.gcfmask3.eq.1)
            lonmin = lonmin + 360.
            dateline_case = 1
      end where
      where(gcfmask1.eq.1.and.gcfmask2.eq.0.and.gcfmask3.eq.0)
            dateline_case = 0
      ENDWHERE
      print *, 'dateline_case'
      print *, dateline_case

      sumofweights = 0.
      value_sow = 0.
c      print *, sumofweights

      work1( 1:nobs ) = 0.
      work3( 1:nobs ) = 0.

      DO j = 1 , gc_failures
c        FORALL(i=1:nobs, latitude(i).le.latmax(j) .and. latitude(i).ge.
c     &                   latmin(j) .and. dattypw(i) .eq. 1. and.
c     &           (( dateline_case(j).eq.0).and. longitude(i).gt.lonmin(j)
c     &                                 .and. longitude(i).lt.lonmax(j))
c     &    .or.  (dateline_case(j).eq.1.and. longitude(i).gt.lonmin(j)
c     &                                 .or.  longitude(i).lt.lonmax(j))
c     &           .or.   polar(j).eq.1 )) 

         do i=1,nobs
c       print *, 'ready to test', i, nobs
c       print *, latitude(i),dattypw(i),longitude(i)
c       print *, latmax(j),latmin(j),dateline_case(j),lonmax(j),polar(j)

        if(latitude(i).le.latmax(j) .and. latitude(i).ge.
     &              latmin(j) .and. dattypw(i) .eq. 1. and.
     &      (( dateline_case(j).eq.0.and. longitude(i).gt.lonmin(j)
     &                            .and. longitude(i).lt.lonmax(j))
     &  .or.  (dateline_case(j).eq.1.and. longitude(i).gt.lonmin(j)
     &                               .or.  longitude(i).lt.lonmax(j))
     &        .or.   polar(j).eq.1 )) then
             print *,'gcf   ',i,j, latmin(j),lonmin(j),lonmax(j),
     &                            latmax(j)
             print *,'obs   ',i,j, latitude(i),longitude(i)

c            print *, 'ready to do clat1 ', i,latitude(i),drc
            clat1(i) = cos(latitude(i) * drc)
c            print *, 'finished clat1 ', i
c            print *, 'ready to do clat2 ', i,gcflat(j),drc
            clat2(i) = cos(gcflat(j) * drc)
c            print *, 'finished clat2 ', i
c            print *, 'ready to do clond ',i,longitude(i),gcflon(j),drc
            clond(i) = cos((longitude(i) - gcflon(j)) * drc)
c            print *, 'finished clond ', i
c            print *, 'ready to do cdif1 ',i,latitude(i),gcflat(j),drc
            cdif1(i) = cos((latitude(i) - gcflat(j)) * drc)
c            print *, 'finished cdif1 ', i
c            print *, 'ready to do temp ',i,cdif1(i),clat1(i),clat2(i),
c     &                clond(i)
            temp(i) = 1.- cdif1(i)+clat1(i)+clat2(i)*(1.-clond(i))
c            print *, 'finished temp ', i
c         print *,i,j
c         print *, clat1(i),clat2(i),clond(i),cdif1(i),temp(i)
c            print *, 'ready to do distance ',i,temp(i)
            distance(i) = sqrt(2. * temp(i)) * rad_e
c            print *, 'finished distance ', i
c            print *, 'ready to do weight',i,distance(i),alpha
            weight(i) = (alpha - distance(i)) / (alpha + distance(i))
c            print *, 'finished weight ', i
c            print *, 'ready to do work1',i,weight(i)
            work1(i) = weight(i)
c            print *, 'finished work1 ', i
c            print *, 'ready to do work3',i,weight(i),datslp(i)
            work3(i) = weight(i) * datslp(i)
c            print *, 'finished work3 ', i
c         print *, distance(i),weight(i),work1(i),work3(i)
c         END FORALL
         END if
        end do
          sumofweights(j)=sum(work1(1:nobs),dim=1)
          value_sow(j)=sum(work3(1:nobs),dim=1)
      ENDDO
      reacc=0
      print *, 'sumofweights'
      print *, sumofweights
      print *, 'obserr1'
      print *, obserr1(gcfikx)
      print *, 'gcfslp'
      print *, gcfslp
      WHERE (sumofweights .GT. 0.)
            analyzed_value = value_sow / sumofweights
      end where
      print *, 'analyzed_value'
      print *, analyzed_value
      print *, (gcfslp - analyzed_value) ** 2 .LT.
     &        (obserr1(gcfikx) ** 2 + gcffcerr**2)
     &           * slp_tolerance
      print *, 'gcftypw'
      print *, gcftypw
      WHERE (((gcfslp - analyzed_value) ** 2 .LT.
     &        (obserr1(gcfikx) ** 2 + gcffcerr**2)
     &           * slp_tolerance) 
     &          .and.
     &      (sumofweights .GT. 0.)
     &          .and.
     &      (gcftypw .EQ. 0.1))
                 gcftypw = 1.2
                 reacc =  1
      end where
      WHERE (((gcfslp - analyzed_value) ** 2 .LT.
     &        (obserr1(gcfikx) ** 2 + gcffcerr**2)
     &           * slp_tolerance)
     &          .and.
     &      (sumofweights .GT. 0.)
     &          .and.
     &      (gcftypw .EQ. 0.3))
                 gcftypw = 1.3
                 reacc = 1
      ENDWHERE
      reaccepts=sum(reacc)
C  Find high quality neighbouring observations.
      neighbours = 0
      WRITE( 6,'(//," QCSLP: ",i8," pressure reports reaccepted by",    
     &  " the buddy check.")' ) reaccepts
C  Report quality control results.
      data_mm(gcfpoint,ind_type) = gcftypw
c      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)
     
C  Perform a similar quality control for the winds.
C  Obtain first-guess sea-level wind values at the locations of the
C  observations.

      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
      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
      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.
      deallocate(gcflat)
      deallocate(gcflon)
      deallocate(gcftypw)
      deallocate(gcfslp)
      deallocate(gcfikx)
      deallocate(gcffcerr)
      deallocate(gcfpoint)

      deallocate(polar)
      deallocate(latmax)
      deallocate(latmin)
      deallocate(lonmax)
      deallocate(lonmin)
      deallocate(londif)
      deallocate(dateline_case)
      deallocate(sumofweights)
      deallocate(analyzed_value)
      deallocate(value_sow)
      deallocate(reacc)

      neighbours = 0

      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)

      allocate(polar(gc_failures))
      allocate(latmax(gc_failures))
      allocate(latmin(gc_failures))
      allocate(lonmax(gc_failures))
      allocate(lonmin(gc_failures))
      allocate(londif(gc_failures))
      allocate(dateline_case(gc_failures))
      allocate(sumofweights(gc_failures))
      allocate(delta_slwfail(gc_failures))
      allocate(analyzed_value(gc_failures))
      allocate(testval(gc_failures))
      allocate(value_sow(gc_failures))
      allocate(reacc(gc_failures))


      datslu=data_mm(1:nobs,ind_uwnd)
      datslv=data_mm(1:nobs,ind_vwnd)

      gcflat = latitude(gcfpoint)
      gcflon = longitude(gcfpoint)
      gcftypw=dattypw(gcfpoint)
      gcfslu=datslu(gcfpoint)
      gcfslv=datslv(gcfpoint)
      gcfikx=int(kx(gcfpoint))
      gcffcerr=datfcerr(gcfpoint)

      polar = 0
      latmax = gcflat + 10.
      latmin = gcflat - 10.
      WHERE(gcflat .GT. 84.9)
          polar = 1
          latmax = 90.
          latmin = 70.
      end where
      WHERE (gcflat .LT. -84.9)
          polar = 0
          latmax = -70.
          latmin = -90.
      end where
      WHERE(polar.eq.0) 
          londif = atan2(.3139,cos(gcflat * drc)) * rdc * .5
          lonmax = gcflon + londif
          lonmin = gcflon - londif
      end where
C  Take care of points near the international date line.
      WHERE ((lonmax .GT. 180.).and.(polar.eq.0))
            lonmax = lonmax - 360.
            dateline_case = 1
      end where
      WHERE ((lonmin .LT. -180.).and.(polar.eq.0))
            lonmin = lonmin + 360.
            dateline_case = 1
      end where
      wHERE((polar.eq.0).and.(lonmin .ge. -180.).and.(lonmax .le. 180.))
            dateline_case = 0
      ENDWHERE

      sumofweights = 0.
      value_sow = 0.
      work1 = 0.
      work3 = 0.

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

      DO j = 1 , gc_failures
          WHERE(latitude.LE.latmax(j).AND.latitude.GE.latmin(j).AND.
     &    search_typ.GE.2. .AND.(dateline_case(j).eq.1.AND.longitude
     &    .GT.lonmin(j).OR.longitude.LT.lonmax(j).OR.polar(j).eq.1))
            clat1 = cos(latitude * drc)
            clat2 = cos(gcflat(j) * drc)
            clond = cos((longitude - gcflon(j)) * drc)
            cdif1 = cos((latitude - gcflat(j)) * 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
          ENDWHERE
          sumofweights(j)=sum(work1,dim=1)
          value_sow(j)=sum(work3,dim=1)
      ENDDO

      reacc=0
      WHERE(sumofweights .GT. 0.) 
          analyzed_value = value_sow / sumofweights
          delta_slwfail = sqrt(gcfslu**2 + gcfslv**2)
      end where
      WHERE(sumofweights .GT. 0.and.
     &    ((delta_slwfail - analyzed_value) ** 2 .LT. 
     &    (obserr3(gcfikx) ** 2 + gcffcerr** 2)*slw_tolerance))
                testval = gcftypw
      end where
      WHERE(sumofweights .GT. 0.and.
     &     ((int(testval*100.)-int(testval*10.)*10.).EQ.3.))
              gcftypw = gcftypw + 2.
              reacc =  1
      ELSEWHERE
              gcftypw = gcftypw + 2.01
              reacc = 1
      ENDWHERE
      reaccepts=sum(reacc)
      PRINT *, reaccepts,' wind reports reaccepted by buddy check.'
c      CALL sqcrpt(data_mm,nobs,2,nymd,nhms)
      data_mm(gcfpoint,itypw) = gcftypw
      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
      forall(i=1:nobs) nze(i)=sum(nz(1:i),dim=1)
      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000) 
     &             ze(nze(i))=zze(i)
      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000) 
     &             zelat(nze(i))=zzelat(i)
      forall(i=1:nobs,warn(i).AND.nze(i).lt.20000) 
     &             zelon(nze(i))=zzelon(i)
      maxnze=maxval(nze)
      IF (maxnze .GE. 20000) THEN
          message = ' Excessive GEN 1000MB heights. First 20000 used.'
          CALL flash(3,message,1)
          exit = .TRUE.
      ENDIF

c  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
C  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)

      deallocate(polar)
      deallocate(latmax)
      deallocate(latmin)
      deallocate(lonmax)
      deallocate(lonmin)
      deallocate(londif)
      deallocate(dateline_case)
      deallocate(sumofweights)
      deallocate(analyzed_value)
      deallocate(value_sow)
      deallocate(reacc)

      RETURN
      END
