************************************
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
      real,dimension(1:nobs):: clat1,clat2,clond,cdif1,distance,weight
!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)
      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)
      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?     ***********
      integer :: londif
      real :: sumofweights
      real :: value_sow
      real :: analyzed_value
      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(reacc(gc_failures))

      gcflat = latitude(gcfpoint)
      print *, 'gcflat'
      print *, gcflat
      gcflon = longitude(gcfpoint)
      print *, 'gcflon'
      print *, gcflon
      gcftypw=dattypw(gcfpoint)
      print *, 'gcftypw'
      print *, gcftypw
      gcfslp=datslp(gcfpoint)
      print *, 'gcfslp'
      print *, gcfslp
      gcfikx=int(kx(gcfpoint))
      print *, 'gcfikx'
      print *, gcfikx
      gcffcerr=datfcerr(gcfpoint)
      print *, 'gcffcerr'
      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
        print *, i,polar,gcflat(i),gcflon(i),drc,rdc,
     c    cos( gcflat(i)*drc ),atan2( .3139, cos( gcflat(i)*drc ) ),
     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.
 
      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(j)
         end where
 
         sumofweights = sum(work1)
         value_sow    = sum(work3)
 
        if ( sumofweights .gt. 0. ) then

          analyzed_value = value_sow / sumofweights
          print *, i,gcfslp(i),analyzed_value
          print *, i,obserr1( gcfikx(i))
          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
          print *, i, gcftypw(i), reaccepts


 
        end if
 
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
      open(unit=33,file='data2.out',status='unknown')
      write(33,*) data_mm(1:nobs,ind_type)
 
!  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.

      CALL terpv(latitude,longitude,work1,nobs,slu)
      CALL terpv(latitude,longitude,work2,nobs,slv)

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


      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.


      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(reacc(gc_failures))

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

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

      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.
 
        do j = 1,nobs
 
 
          if ( latitude(j) .le. latmax .and. latitude(j) .ge.
     &         latmin .and. search_typ(j) .ge. 2. ) then
 
 
             if (( .not.(dateline_case).and. longitude(j).gt.lonmin
     &                                 .and. longitude(j).lt.lonmax)
     &           .or.  (dateline_case  .and. longitude(j).gt.lonmin
     &                                 .or.  longitude(j).lt.lonmax)
     &           .or.   polar ) then

                  clat1(j) = cos ( search_lat * drc )
                  clat2(j) = cos ( data(ifail+ilat-1)*drc )
                  clond(j) = cos ((search_lon - data(ifail+ilon-1))*drc)
                  cdif1(j) = cos ((search_lat - data(ifail+ilat-1))*drc)
                  temp(j)  = 1.-cdif1+ clat1+clat2*(1.-clond)
                  distance(j) =  sqrt ( 2.*temp ) * rad_e 
 
                  weight(j)=(alpha-distance(j))/(alpha+distance(j))
                  work1(j) = weight(j)
                  work3(j) = weight(j) * search_slw(j)
 
             end if
 
          end if
 
        end do
 
        sumofweights = sum(work1)
        value_sow    = sum(work3)
 
        if ( sumofweights .gt. 0. ) then
 
          analyzed_value = value_sow / sumofweights
          delta_slwfail = sqrt( gcfslu(i)**2 + gcfslv(i)**2)
 
          if (  (delta_slwfail - analyzed_value)**2 .lt.
     &          (obserr3(gcfikx(i) )**2 +
     &           gcffcerr(i)**2) * slw_tolerance ) then
 
               testval = gcftypw(i )
               if ((int(testval*100.)-int(testval*10.)*10.).eq.3.)then
                 gcftypw(i ) = gcftypw(i) + 2.
                 reaccepts = reaccepts + 1
               else
                 gcftypw(i ) = gcftypw(i)  + 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.'
      data_mm(gcfpoint,itypw) = gcftypw

      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
      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
  Write out to DEL dataset.
      WRITE( 91 ) nobs,nymd,nhms,'SLP'
      DO i = 1 , nobs
        WRITE( 91 ) data_mm(i,ind_lat),data_mm(i,ind_lon),
     &  data_mm(i,ind_type),data_mm(i,ind_slp),data_mm(i,ind_kx),
     &  data_mm(i,islu),data_mm(i,islv),data_mm(i,ind_slp)
      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)

      RETURN
      END
