                                                                        
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
!                                                                     .
!  Program:    QCSLP                                                  .
!  Language:   Fortran 77 (with Fortran 90 enhancements)              .
!  Computer:   Cray YMP 864                                           .
!                                                                     .
!  Programmer: M. Seablom                                             .
!  Date:       December 16, 1991
!                                                                     .
!  Program QCSLP  computes sea-level pressure and wind deltas, and    .
!  performs version 2.1 quality control.                              .
!                                                                     .
!                                                                     .
!......................................................................
 
 
      subroutine qcslp  ( phis, slp, slu, slv, data, nobs, t, ze, zelat,
     &                    zelon, nze, fgslp, suspect, fcerr )
 
 
      real           phis(im,jnp), slp(im,jnp), slu(im,jnp),
     &               slv( im,jnp), data(1500000      ), zelat(20000),
     &               zelon(20000), t(im,jnp,mlev), fgslp( im,jnp ),
     &               pm( nobs ), work1( nobs ), work2( nobs ), 
     &               observation_error, fcerr( im,jnp ), work3( nobs ),
     &               latmax, latmin, lonmax, lonmin, distance, weight,
     &               alpha, ze(20000), phis1(im,jnp)
 
      integer        gc_failures, gcfpoint( 30000 ), reaccepts
 
      common /sfcdrag/  drag(13104), iwaterflag(13104)
      common /dataptrs/ ikx,ilat,ilon,iprs,ityp,islp,ifcerr
      common /cortabs/  hcorfe ( 1000, 3, 20, 2 )
      common /scan/     maxobs, rscan, onestage, max_data
      common /gridxx/   im, jnp, mlev
      common / obs /    maxdex, nfgge, lun, dataname(5), hh(13,13,20),
     &                  vcoroe(18,18,4), obserr(200,4,19), ww(13,13,20),
     &                  obsname(200), qobs(200,5), offlev(200)
 
 
      real           indx,kx,lat,lon,xdat1,xdat2,diff,type
 
      logical        exit, onestage, suspect(150000,2),    polar, obs,
     &               offlev, dateline_case
      character *60  message
      character * 8  dataname, obsname
 
      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
 
 
!  Define additional pointers to the DATA array.
 
 
      islu  = islp + nobs
      islv  = islu + nobs
      ifcerr= islv + nobs
      itypw = ityp
 
      call getconi ( 'NYMD    ',nymd )
      call getconi ( 'NHMS    ',nhms )
      call getconh ( 'RHOBAR  ',rhobar )
      call getconh ( 'RADEI   ',radei )
      call getconh ( 'OMEGA2  ',omega2 )
 
      do i = 1,nobs
        data ( itypw+i-1 ) = 0.
      end do
 
 
!  Obtain first-guess pressure values at the locations of the
!  observations.
 
      call terpv ( data(ilat), data(ilon), work1, nobs, slp(1,1) )
 
!  Interpolate forecast errors to observation locations.
 
      call terpv ( data(ilat), data(ilon), data(ifcerr), nobs, fcerr )
 
!  Establish loop to compute deltas and perform the gross check of
!  the sea-level pressure data.
 
      gc_failures = 0
 
      do i = 1,nobs    !  *** Begin SLP gross check loop *******
 
        if ( data( islp+i-1 ) .lt. fill ) then
            data( itypw+i-1 ) = 1.
        else
            data( itypw+i-1 ) = 0.4
        end if
 
        latitude = data( ilat+i-1 )
        longitude= data( ilon+i-1 )
        kx       = data( ikx +i-1 )
 
        if ( latitude .gt. 90.  .or.  latitude .lt. -90. .or.
     &       longitude .gt. 180. .or. longitude .lt. -180. ) then
 
           message =  'QCSLP:  Impossible data location; rejected'
     &               // 'point.'
           data( itypw+i-1 ) = 0.4
           call flash ( 3,message,1 )
 
        end if
 
        if ( data( itypw+i-1 ) .eq. 1. ) data( islp+i-1 ) =
     &         data( islp+i-1 ) - work1( i )
 
!  Correction for erroneous ship and buoy observations.
 
c       if ( data( itypw+i-1 ) .eq. 1. .and. data( ikx+i-1 ) .ge. 3.
c    &     .and. data( ikx+i-1 ) .le. 6. .and. data( islp+i-1 ) .ge.
c    &     50. ) data( islp+i-1 ) = data( islp+i-1 ) - 100.
 
c       if ( data( itypw+i-1 ) .eq. 1. .and. data( ikx+i-1 ) .ge. 3.
c    &     .and. data( ikx+i-1 ) .le. 6. .and. data( islp+i-1 ) .le.
c    &    -50. ) data( islp+i-1 ) = data( islp+i-1 ) + 100.
 
 
        delta = data( islp+i-1 )
        observation_error = obserr( int(kx), 1, 19 )
        forecast_error = data( ifcerr+i-1 )
 
        if ( delta**2 .gt. ( observation_error**2 + forecast_error**2 )
     &      * slp_tolerance .and. data(itypw+i-1) .ne. .4) then
 
            data( itypw + i - 1 ) = 0.1
            work3( i ) = 1.
 
        end if
 
!  Suspect observations automatically fail the gross check.
 
        if ( suspect(i,1) .and. data( islp+i-1 ) .lt. fill ) then
            data( itypw + i - 1 ) = 0.3
            work3( i ) = 1.
        end if
 
      end do           !  ***  End SLP gross check loop ********
 
      call wheneq ( nobs, work3, 1, 1., gcfpoint, gc_failures )
 
      if ( gc_failures .gt. 30000 ) then
        message = 'QCSLP:  Excessive gross check failures.'
        call flash( 5,message,1 )
      end if
 
      write( 6,* ) ' QCSLP:  ', gc_failures, ' sea-level pressure',
     &             ' observations have failed the gross check.'
 
      reaccepts = 0
 
!  Analyze data surrounding an observation which failed the gross
!  check. Search for observations within 1000km of the suspect point.
 
      do 2 i = 1,gc_failures   !  ***  Begin SLP buddy check loop ****
 
        polar = .false.
        ifail = gcfpoint( i )
        gcflat = data( ifail+ilat-1 )
        gcflon = data( ifail+ilon-1 )
 
 
        if ( gcflat .gt. 84.9 ) then
            polar = .true.
            latmax = 90.
            latmin = 70.
        else if ( gcflat .lt. -84.9 ) then
            polar = .true.
            latmax = -70.
            latmin = -90.
        else
            latmax = gcflat  + 10.
            latmin = gcflat  - 10.
        end if
 
        if ( .not. polar ) then
 
            londif = atan2( .3139, cos( gcflat*drc ) ) * rdc * .5
      print *, i,polar,gcflat,gcflon,drc,rdc,cos( gcflat*drc),
     >   atan2( .3139, cos( gcflat*drc ) ),londif
            lonmax = gcflon + londif
            lonmin = gcflon - 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
 
          search_lat = data( ilat+j-1 )
          search_lon = data( ilon+j-1 )
          search_slp = data( islp+j-1 )
          search_typ = data( itypw+j-1 )
 
          if ( search_lat .le. latmax .and. search_lat .ge.
     &         latmin .and. search_typ .eq. 1. ) then
 
             if (( .not.(dateline_case).and. search_lon.gt.lonmin
     &                                 .and. search_lon.lt.lonmax)
     &           .or.  (dateline_case  .and. search_lon.gt.lonmin
     &                                 .or.  search_lon.lt.lonmax)
     &           .or.   polar ) then

 
                  clat1 = cos ( search_lat * drc )
                  clat2 = cos ( data(ifail+ilat-1)*drc )
                  clond = cos ((search_lon - data(ifail+ilon-1))*drc)
                  cdif1 = cos ((search_lat - data(ifail+ilat-1))*drc)
                  temp  = 1.-cdif1+ clat1+clat2*(1.-clond)
                  distance =  sqrt ( 2.*temp ) * rad_e
 
                  weight = ( alpha - distance ) / ( alpha + distance )
                  work1(j) = weight
                  work3(j) = weight * search_slp
 
             end if
 
          end if
 
        end do
 
        do isum = 1,nobs
 
          sumofweights = sumofweights + work1( isum )
          value_sow    = value_sow    + work3( isum )
 
        end do
 
        if ( sumofweights .gt. 0. ) then
 
          analyzed_value = value_sow / sumofweights
          print *, i,data( ifail+islp-1 ),analyzed_value
          print *, i,obserr( int(data( ifail+ikx-1 )),1,19 )
          print *, i,data(ifcerr+ifail-1),slp_tolerance
 
          if (  (data( ifail+islp-1 ) - analyzed_value)**2 .lt.
     &          (obserr( int(data( ifail+ikx-1 )),1,19 )**2 +
     &           data(ifcerr+ifail-1)**2) * slp_tolerance ) then
 
 
               if ( data( ifail+itypw-1 ) .eq. 0.1 ) then
                 data( ifail+itypw-1 ) = 1.2
                 reaccepts = reaccepts + 1
               else if ( data( ifail+itypw-1) .eq. 0.3 ) then
                 data( ifail+itypw-1 ) = 1.3
                 reaccepts = reaccepts + 1
               end if
 
          end if
          print *, i, data( ifail+itypw-1), reaccepts
 
        end if
 
2     end do                   ! *** End SLP buddy check loop ********
 
      write(6,'(//," QCSLP: ",i8," pressure reports reaccepted by",
     &      " the buddy check.")' ) reaccepts
 
!  Report quality control results.
      open(unit=33,file='data1.out',status='unknown')
      write(33,*) data(1 + 4 * nobs:5*nobs)

 
      call sqcrpt ( data, nobs, 1, nymd, nhms )
 
      do i = 1,nobs
        data( itypw+i-1 ) = int( data( itypw+i-1 ) )
      end do
 
      call rptlly ( data(ikx), data(itypw), nobs, 1 )
 
!  Perform a similar quality control for the winds.
 
!  Obtain first-guess sea-level wind values at the locations of the
!  observations.
 
      call terpv ( data(ilat), data(ilon), work1, nobs, slu(1,1) )
      call terpv ( data(ilat), data(ilon), work2, nobs, slv(1,1) )

!  Reset ocean points in the PHIS1 array to zero.

c     do i = 1,im*jnp
c         phis1(i) = phis(i) * float( iwaterflag(i) )
c     end do
 
!  Interpolate the PHIS1 (terrain) array to observation locations.
 
c     call terpv ( data(ilat), data(ilon), work3, nobs, phis1(1,1) )
 
!  Filter out fill values and wind observations over land.
 
      gc_failures = 0
 
      do i = 1,nobs 
 
        if ( data( islu+i-1 ) .lt. fill .and. data( islv+i-1 )
     &       .lt. fill ) then
 
            data( itypw+i-1 ) = data( itypw+i-1 ) + 2.
        else
            data( itypw+i-1 ) = data( itypw+i-1 ) + 0.04
 
        end if
 
      end do

!  Establish loop to compute deltas and perform the gross check.

      work3(1:nobs) = 0.

      do i = 1,nobs

        latitude = data( ilat+i-1 )
        longitude= data( ilon+i-1 )
        kx       = data( ikx +i-1 )
 
        if ( latitude .gt. 90.  .or.  latitude .lt. -90. .or.
     &       longitude .gt. 180. .or. longitude .lt. -180. ) then
 
           message =  'QCSLP : Impossible data location; rejected'
     &               // ' point.'
           data( itypw+i-1 ) = 0.4
           call flash ( 3,message,1 )
 
        end if
 
        if ( data( itypw+i-1 ) .ge. 2. ) then
            data( islu+i-1 ) = data( islu+i-1 ) - work1( i )
            data( islv+i-1 ) = data( islv+i-1 ) - work2( i )
        delta = sqrt( data( islu+i-1 )**2 + data( islv+i-1 )**2 )
        end if
 
 
        observation_error = obserr( int(kx), 3, 19 )
 
!  Calculate the wind forecast error.
 
 
        if ( abs( latitude ) .lt. 1.e-5 ) latitude = 1.e-5
        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( ifcerr+i-1 ) * sqrt(-hcorfe( 1,3,1,1 ))
 
        if ( delta**2 .gt.(( observation_error**2 + forecast_error**2 )
     &      * slw_tolerance) .and. data(itypw+i-1) .ge. 2. .and.
     &      (.not.suspect(i,2)) ) then
 
            data( itypw + i - 1 ) = data( itypw+i-1 ) - 1.99
            work3( i ) = 1.
        else
            work3( i ) = 0.
 
        end if
 
!  Suspect observations automatically fail the gross check.
 
        if ( suspect(i,2) .and. data( islu+i-1 ) .lt. fill ) then
            data( itypw + i - 1 ) = data( itypw+i-1 ) - 1.97
            work3( i ) = 1.
        end if
 
      end do           !  ***  End SLP gross check loop ********
 
      call wheneq ( nobs, work3, 1, 1., gcfpoint, gc_failures )
 
      if ( gc_failures .gt. 30000 ) then
        message = 'QCSLP:  Excessive gross check failures.'
        call flash( 5,message,1 )
      end if
 
      write( 6,* ) ' QCSLP:  ', gc_failures, ' sea-level wind',
     &             ' observations have failed the gross check.'
 
      reaccepts = 0
 
!  Analyze data surrounding an observation which failed the gross
!  check. Search for observations within 1000km of the suspect point.
 
      do 3 i = 1,gc_failures   !  ***  Begin SLW buddy check loop ****
 
        polar = .false.
        ifail = gcfpoint( i )
        gcflat = data( ifail+ilat-1 )
        gcflon = data( ifail+ilon-1 )
 
 
        if ( gcflat .gt. 84.9 ) then
            polar = .true.
            latmax = 90.
            latmin = 70.
        else if ( gcflat .lt. -84.9 ) then
            polar = .true.
            latmax = -70.
            latmin = -90.
        else
            latmax = gcflat  + 10.
            latmin = gcflat  - 10.
        end if
 
        if ( .not. polar ) then
 
            londif = atan2( .3139, cos( gcflat*drc ) ) * rdc * .5
            lonmax = gcflon + londif
            lonmin = gcflon - 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
 
          search_lat = data( ilat+j-1 )
          search_lon = data( ilon+j-1 )
          search_slw = sqrt(data( islu+j-1 )**2 + data(islv+j-1)**2 )
          search_typ = data( itypw+j-1 )
 
          if ( search_lat .le. latmax .and. search_lat .ge.
     &         latmin .and. search_typ .ge. 2. ) then
 
 
             if (( .not.(dateline_case).and. search_lon.gt.lonmin
     &                                 .and. search_lon.lt.lonmax)
     &           .or.  (dateline_case  .and. search_lon.gt.lonmin
     &                                 .or.  search_lon.lt.lonmax)
     &           .or.   polar ) then

                  clat1 = cos ( search_lat * drc )
                  clat2 = cos ( data(ifail+ilat-1)*drc )
                  clond = cos ((search_lon - data(ifail+ilon-1))*drc)
                  cdif1 = cos ((search_lat - data(ifail+ilat-1))*drc)
                  temp  = 1.-cdif1+ clat1+clat2*(1.-clond)
                  distance =  sqrt ( 2.*temp ) * rad_e 
 
                  weight = ( alpha - distance ) / ( alpha + distance )
                  work1(j) = weight
                  work3(j) = weight * search_slw
 
             end if
 
          end if
 
        end do
 
        do isum = 1,nobs
 
          sumofweights = sumofweights + work1( isum )
          value_sow    = value_sow    + work3( isum )
 
        end do
 
        if ( sumofweights .gt. 0. ) then
 
          analyzed_value = value_sow / sumofweights
          delta_slwfail = sqrt( data(ifail+islu-1)**2 +
     &                       data(ifail+islv-1)**2 )
 
          if (  (delta_slwfail - analyzed_value)**2 .lt.
     &          (obserr( int(data( ifail+ikx-1 )),3,19 )**2 +
     &           data(ifail+ifcerr-1)**2) * slw_tolerance ) then
 
               testval = data( ifail+itypw-1 )
               if ((int(testval*100.)-int(testval*10.)*10.).eq.3.)then
                 data( ifail+itypw-1 ) = data( ifail+itypw-1 ) + 2.
                 reaccepts = reaccepts + 1
               else
                 data( ifail+itypw-1 ) = data(ifail+itypw-1) + 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, nobs, 2, nymd, nhms )
 
      do i = 1,nobs
        data( itypw+i-1 ) = int( data( itypw+i-1 ) )
      end do
 
      work1 ( 1:nobs ) = 0.
      work2 ( 1:nobs ) = 0.
 
!  Create 1000mb heights based upon lowest level first guess temperature
!  and surface pressure.
 
      R_d = 287.
      g = 9.81
 
!  Interpolate temperatures to locations of observations.
 
 
      call terpv ( data(ilat), data(ilon), work1, nobs, t(1,1,mlev) )
 
      call terpv ( data(ilat), data(ilon), work2, nobs, fgslp)
 
      nze = 0
      exit = .false.
      i = 1
 
      do while ( .not. exit )
 
        if ( data(ikx+i-1) .ge. 3.  .and.  data(ikx+i-1) .le. 6.
     &     .and. data(islp+i-1) .lt. 1.e10 ) then
 
           nze = nze + 1
 
           ze( nze ) = ( R_d * work1(i) / g ) * ((1000./work2(i)) *
     &                  ( data(islp+i-1) / work2(i) ) )
 
           zelat( nze ) = data( ilat+i-1 )
           zelon( nze ) = data( ilon+i-1 )
 
 
        end if
 
        if ( nze .ge. 20000 ) then
         message = ' Excessive GEN 1000MB heights. First 20000 used.'
         call flash ( 3,message,1 )
         exit = .true.
        end if
 
        if ( i .ge. nobs ) exit = .true.
        i = i + 1
 
      end do
 
c  Write out to DEL dataset.
 
      write(91) nobs, nymd, nhms, 'SLP'
      do i = 1,nobs
 
        write(91) data(ilat+i-1), data(ilon+i-1), data(itypw+i-1),
     &           data(islp+i-1), data(ikx+i-1),  data(islu+i-1),
     &           data(islv+i-1), data(islp+i-1)
 
      end do
 
 
!  Tally all observations.
 
      call rptlly ( data(ikx), data(itypw), nobs, 1 )
      call rptlly ( data(ikx), data(itypw), nobs, 2 )
 
      return
 
      end
