c---------------------------------------------------------------------c
c         N A S A   /  G M S B   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    Program Name: SQCRPT ( Surface Quality Control Report )          c
c    Programmer  : Michael Seablom, NASA/GSFC Code 911                c
c                  Greenbelt, MD 20771  (301)286-7360 (FTS)888-7360   c
c                                                                     c
c    Date        : August 22, 1991                                    c
c    Computer    : Cray YMP 464                                       c
c    Language    : Cray Fortran (CFT77 Version 4.0)                   c
c                                                                     c
c  Program SQCRPT prints out surface quality control statistics.      c
c  This subroutine is patterned after QCREPORT which performs similar c
c  duties in the upper air analysis.                                  c
c                                                                     c
c---------------------------------------------------------------------c
C**********************************************************************
C  Data parallel implementation (CM Fortran) : M. Makivic (NPAC)      *
C  Date:       4/94                                                   *
C**********************************************************************


      SUBROUTINE sqcrpt(data_mm,nobs,type,nymd,nhms)
      include 'params.h'
      REAL data_mm(nobs,NRECORDS)
      INTEGER nobs, statistics (mxqobs, 10), type
      CHARACTER * 6 analysis_time, date_temp, time_temp
      CHARACTER * 8 dataname, obsname, analysis_date
      LOGICAL qobs, offlev, onestage
      COMMON / scan / maxobs, rscan, onestage, MaxData
      COMMON / dataptrs / ikx, ilat, ilon, iprs, ityp, idel, ifcer
      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)

!hpf$ distribute data_mm(block,block)
!hpf$ distribute hh(block,block,block)
!hpf$ distribute obserr(block,block,block)
!hpf$ distribute offlev(block)
!hpf$ distribute qobs(block,block)
!hpf$ distribute statistics(block,block)
!hpf$ distribute vcoroe(block,block,block)
!hpf$ distribute ww(block,block,block)

      integer,dimension(nobs):: stat
!hpf$ distribute stat(block)
      integer,dimension(nobs):: itype
!hpf$ distribute itype(block)
      integer,dimension(nobs):: testval
!hpf$ distribute TESTVAL(block)
      real,dimension(nobs):: fe_acc
!hpf$ distribute fe_acc(block)
      real,dimension(nobs):: fe_f
!hpf$ distribute fe_f(block)
      real,dimension(nobs) :: delta_acc
!hpf$ distribute delta_acc(block)
      real,dimension(nobs) :: delta_f
!hpf$ distribute delta_f(block)
      integer,dimension(nobs):: number_acc
!hpf$ distribute number_acc(block)
      integer,dimension(nobs):: number_f
!hpf$ distribute number_f(block)
      integer, parameter:: ind_del=ind_slp

      WRITE( 6,1001)
1001  FORMAT('--------------------------------------------------------')
1000  FORMAT (//,5x,"Data",8x,"Total",6x,"Number",4x," Suspect Observat"
     &,"ions",2x,"Gross Check  Reaccepted ",4x,"Mean Fcst Error",10x,"RM
     &S - Delta",/,5x,"Type",7x,"Reports",4x,"Accepted",5x,"Total",6x,"A
     &ccepted",4x,"Failures",2x,"Buddy Check",5x,"Pass",8x,"Fail",8x,"Pa
     &ss",8x,"Fail",/)
      WRITE( 6,1001)
      WRITE( 6,1002)
1002  FORMAT(14x,"Quality Control, Version 2.1, Autumn, 1991")
      WRITE( 6,1001)
      WRITE( 6,1001)
      WRITE( date_temp,'(i6)' ) nymd
      analysis_date = date_temp(3:4) // '/' // date_temp(5:6) // '/' 
     &    // date_temp(1:2)
      WRITE( time_temp,'(bz,i6)' ) nhms
      analysis_time = time_temp(1:2) // ':' // time_temp(3:4) // 'Z'
      IF (nhms .EQ. 0) THEN
        analysis_time = ' 0:00Z'
      ENDIF
      WRITE( 6,'(///)' ) 
      IF (type .EQ. 1) THEN
        WRITE( 6,1003)
1003  FORMAT(" SEA LEVEL PRESSURE STATISTICS VARIABLE TOLERANCES")
        WRITE( 6,1004) analysis_date,analysis_time
1004  FORMAT(" Analysis date: ",a8,3x,"Analysis time: ",bz,a6)
      ENDIF
      IF (type .NE. 1) THEN
        WRITE( 6,1005)
1005  FORMAT(" OCEAN SURFACE WIND STATISTICS VARIABLE TOLERANCES")
        WRITE( 6,1004) analysis_date,analysis_time
      ENDIF
      WRITE( 6,1000 ) 
      statistics = 0

      DO j = 1 , 200
          stat=0
          statistics(j,1)=0
c Tally the total number of observations for each data type.
          WHERE (data_mm(:,ind_kx ) .EQ. float(j)) 
            stat = 1
          ELSEWHERE
            stat = 0
          ENDWHERE
          statistics(j,1)=sum(stat)
c Total the number of accepted reports.
          WHERE ((int(data_mm(:,ind_type)) .EQ. type .OR. 
     &          int(data_mm(:,ind_type)) .EQ. 3) .AND. data_mm(:,ind_kx)
     &          .EQ. float(j))
            stat=1
          ELSEWHERE
            stat = 0
          ENDWHERE
          statistics(j,2)=0
          statistics(j,2) = sum(stat)
c          print *, 'passed 1,2'
c Number of suspect observations. If type = 1, pressure observations
c are being tallied, otherwise tally wind observations.
c Suspicious pressure observations.
          if ( type .eq. 1 ) then
c Suspicious pressure observations.
            WHERE ( data_mm(:,ind_kx) .eq. float(j) )
              testval = int( data_mm(:, ind_type) * 10. )
            elsewhere
              testval=-100000
            ENDWHERE

            WHERE ( testval .eq. 13 .or. testval .eq. 3 )
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,3)=0
            statistics(j,3)=sum(stat)

            WHERE ( testval .eq. 13 )
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,4)=0
            statistics(j,4) = sum(stat)

c Gross check failures.

            WHERE ( testval .eq. 12 .or. testval .eq. 1 )
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,5)=0
            statistics(j,5) = sum(stat)

c Number reaccepted by buddy check.
            WHERE ( testval .eq. 12)
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,6)=0
            statistics(j,6) = sum(stat)

          ELSE
c Suspicious wind observations.


c            print *, 'inside type 2'
            WHERE ( data_mm(:,ind_kx) .eq. float(j) )
             testval = int( data_mm(:, ind_type) * 100. ) -
     &                int( data_mm( :,ind_type ) * 10.)* 10
            else where
             testval=-10000
            ENDWHERE
c            print *, 'computed testval'

            WHERE ( testval .eq. 3 ) 
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE

            statistics(j,3)=0
            statistics(j,3) = sum(stat)
c            print *, 'computed 3'

            WHERE ( (testval .eq. 3) .and.
     &           (int(data_mm(:,ind_type)) .eq. 2 .or.
     &           int(data_mm(:,ind_type)) .eq. 3) )
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,4)=0
            statistics(j,4) = sum(stat)
c            print *, 'computed 4'

c Gross check failures.

            WHERE ( testval .eq. 2 .or. testval .eq. 1 )
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,5)=0
            statistics(j,5) = sum(stat)
c            print *, 'computed 5'

c Number reaccepted by buddy check.
            WHERE ( testval .eq. 2)
              stat=1
            ELSEWHERE
              stat=0
            ENDWHERE
            statistics(j,6)=0
            statistics(j,6) = sum(stat)
c            print *, 'computed 6'

          ENDIF

c          print *, j, statistics(j,:)
c Average forecast errors and deltas for observations which passed and
c failed the quality control.

       if ( type .eq. 1 ) then

          number_acc =  0
          fe_acc = 0.
          delta_acc = 0.
          number_f =  0
          fe_f = 0.
          delta_f = 0.


          WHERE ( data_mm(:,ind_kx) .eq. float(j) )
             itype = int( data_mm( :,ind_type ) )
          else where
             itype=-100000
          ENDWHERE
 
          WHERE ( itype .eq. 1 .or. itype .eq. 3 )
 
              number_acc =  1
              fe_acc = data_mm(:, ind_fcer)
              delta_acc = (data_mm( :,ind_del))**2

          end where
          WHERE
     &       ((.not.(itype.eq.1.or.itype.eq.3.or.itype.eq.-100000))
     &                .and. 
     &       ( (int( data_mm(:,ind_type)*10.)).ne.4))

              number_f =  1
              fe_f = data_mm(:, ind_fcer )
              delta_f = data_mm(:, ind_del)**2

          ENDWHERE

          fe_accept = sum(fe_acc)
          fe_fail   = sum(fe_f)
          delta_accept = sum(delta_acc)
          delta_fail   = sum(delta_f)
          number_accept = sum(number_acc)
          number_fail   =  sum(number_f)

 
       else
          number_acc =  0
          fe_acc = 0.
          delta_acc = 0.
          number_f =  0
          fe_f = 0.
          delta_f = 0.

 
          WHERE ( data_mm(:,ind_kx) .eq. float(j) ) 
             itype = int( data_mm(:,ind_type) )
          else where
             itype = -10000
          ENDWHERE
c          print *,'computed itype'
 
          WHERE ( itype .eq. 2 .or. itype .eq. 3 )
              number_acc = 1
              fe_acc = data_mm(:, ind_fcer )
              delta_acc = 
     &         (sqrt(data_mm(:,ind_uwnd)**2  +
     &               data_mm(:,ind_vwnd)**2))**2
          end where
c          print *,'computed acc'

          WHERE 
     &          ((.not.( itype .eq. 2 .or. itype .eq. 3 ))
     &          .and.
     &          ((int( data_mm(:,ind_type) * 100. ) -
     &           int( data_mm(:, ind_type) * 10.)* 10.).ne.4))
 
              number_f =  1
              fe_f = data_mm(:, ind_fcer)
              delta_f = 
     &                 (sqrt( data_mm(:,ind_uwnd) **2 +
     &                 data_mm(:,ind_vwnd) ** 2 ))**2
 
          ENDWHERE
c          print *,'computed f'
          fe_accept = sum(fe_acc)
          fe_fail = sum(fe_f)
          delta_accept=sum(delta_acc)
          delta_fail=sum(delta_f)
          number_accept=sum(number_acc)
          number_fail=sum(number_f)
c          print *,'computed accept, fail'
       end if
 
         if ( number_fail .ne. 0 ) then
           fe_fail = fe_fail / float( number_fail )
           delta_fail = sqrt(delta_fail / float( number_fail ))
         end if
         if ( number_accept .ne. 0 ) then
           fe_accept = fe_accept / float( number_accept )
           delta_accept = sqrt(delta_accept / float( number_accept ))
         end if
         if ( statistics(j,1) .ne. 0 ) write(6,'(1x,a8,6i12,4e12.4)')
     &      obsname(j), ( statistics(j,ii),ii=1,6), fe_accept, fe_fail,
     &      delta_accept, delta_fail
 
      ENDDO

      RETURN
      END
