C-----------------------------------------------------------------------
C  MAIN PROGRAM SURUFIT
!   Author: Suranjana Saha
C-----------------------------------------------------------------------
c   make sure when you change levels, you check pmandt and pmandb
C-----------------------------------------------------------------------
      PROGRAM SURUFIT

      PARAMETER (IDBUG=0,IPR=1)
      PARAMETER (NSTC=9)
      PARAMETER (NPLV=3)
      PARAMETER (NVAR=4)
      PARAMETER (NREG=7)
      PARAMETER (NSUB=1)
      PARAMETER (NBAK=2)
 
      CHARACTER*80 HDSTR,OBSTR,FCSTR,ANSTR,QMSTR,PSTR
      CHARACTER*8  SUBSET

      real(8)    HDR(14)
      real(8)    PSOB(4),PSPR(4)
      real(8)    BAK(10,255,NBAK)
      real(8)    OBS(10,255),QMS(10,255)

      real(8)    SPRS(NSTC,NPLV,NVAR,NREG,NSUB,NBAK)
      real(8)    CNTO,CNTN,RAT1,RAT2,WT1,WT2
      real(8)    PMANDB(NPLV),PMANDT(NPLV)
      real(8)    STC(NSTC,5,NBAK)

      real(4)    GDATA(NREG,NSUB)

      LOGICAL    MANDONLY,REGION
      INTEGER    INDEXV(NVAR)


      DATA HDSTR
     ./'SID XOB YOB DHR ELV TYP T29 ITP SQN RQM DUP PRG SRC RUD'/
      DATA PSTR /'POB PAN PFC PQM CAT=0'/
      DATA OBSTR/'POB QOB TOB ZOB UOB VOB'/
      DATA FCSTR/'PFC QFC TFC ZFC UFC VFC'/
      DATA ANSTR/'PAN QAN TAN ZAN UAN VAN'/
      DATA QMSTR/'PQM QQM TQM ZQM WQM CAT'/
 
      real(8) BMISS /10E10/

      DATA RMISS / -9.99E+33 /
      DATA LUBFR/11/
      data indexv/2,3,4,1/
c...  t,z,w,q
c
      DATA PMANDB / 1000,700,300/
      DATA PMANDT /  700,300,150/
c
      levt1=pmandt(1)
      levb1=pmandb(1)
      levt2=pmandt(2)
      levb2=pmandb(2)
      levt3=pmandt(3)
      levb3=pmandb(3)
c
       CALL OPENBF(LUBFR,'IN ',LUBFR)
c
C-----------------------------------------------------------------------
 
C  ZERO THE FIT ARRAYS
C  -------------------
 
      SPRS = 0.
      
      bmiss=10e10; call setbmiss(bmiss) ! this sets bufrlib missing value to 10e10

C  READ AND "SURU-FIT" THE PREPDA/BUFR RECORDS
C  -------------------------------------------
 
10    DO WHILE(IREADMG(LUBFR,SUBSET,IDATE).EQ.0)
c... check for subset...
      IF(ITYP(SUBSET).EQ.0) GOTO 10
11    DO WHILE(IREADSB(LUBFR).EQ.0)

c... check for aircraft only...
      CALL UFBINT(LUBFR,HDR,14,  1,NLEV,HDSTR)
c
c... check for region...
      XOB=hdr(2)
      YOB=hdr(3)
c
      IF(.NOT.REGION(XOB,YOB,0)) GOTO 11
 
C  READ THE DATA
C  -------------
C  GENERATE A PRESSURE LEVEL LOOKUP TABLE
      CALL UFBINT(LUBFR,OBS,10,255,NLEV,OBSTR)
      CALL UFBINT(LUBFR,BAK(1,1,1),10,255,NLFC,FCSTR)
      CALL UFBINT(LUBFR,BAK(1,1,2),10,255,NLAN,ANSTR)
      CALL UFBINT(LUBFR,QMS,10,255,NLQM,QMSTR)
c

C  CREATE AND ACCUMULATE THE STATISTICS ARRAY FOR EACH REALIZATION
C  ---------------------------------------------------------------
 
c... j  is level (21)
c... k  is variable (5)
c... ll is region (7)
c... m  is subset (1)
c... n  is background (2)
c
      POB = OBS(1,1)
      PQM = QMS(1,1)
      CAT = QMS(6,1)
C
      if(idbug.eq.1) print *,' pob ',pob,' pqm ',pqm,' cat ',cat
c
      if(pqm.le.3) then
      j=0
      if((pob.le.levb1).and.(pob.gt.levt1)) j=1
      if((pob.le.levb2).and.(pob.gt.levt2)) j=2
      if((pob.le.levb3).and.(pob.gt.levt3)) j=3

      if(j.gt.0) then
c
      STC = 0.
c... start forecast-loop
         DO IB=1,2
C
c... start variable-loop
         DO IQ=1,4
c
         IQQ = IQ+1
         IR = IQQ+1
C
         IF(QMS(IQQ,1).LE.3 .AND. IQ.LE.3 .AND. CAT.NE.4) THEN
c
            IF(IQ.EQ.1) THEN
            IF(OBS(IQQ,1).EQ.BMISS) GO TO 2234
            IF(IB.EQ.1) OBS(IQQ,1)=OBS(IQQ,1)*1.E-3
            IF(BAK(IQQ,1,IB).EQ.BMISS) GO TO 2234
            BAK(IQQ,1,IB)=BAK(IQQ,1,IB)*1.E-3
            ENDIF
c  count
            STC(1,IQ,IB) = 1.
c  f
            STC(2,IQ,IB) = BAK(IQQ,1,IB)
c  o
            STC(3,IQ,IB) = OBS(IQQ,1)
c  f * o
            STC(4,IQ,IB) = BAK(IQQ,1,IB)*OBS(IQQ,1)
c  f**2
            STC(5,IQ,IB) = BAK(IQQ,1,IB)**2
c  o**2
            STC(6,IQ,IB) = OBS(IQQ,1)**2
C
         ELSEIF(QMS(IQQ,1).LE.3 .AND. IQ.EQ.4) THEN
c  count
            STC(1,IQ,IB) = 1.
c  uf
            STC(2,IQ,IB) = BAK(IQQ,1,IB)
c  vf
            STC(3,IQ,IB) = BAK(IR,1,IB)
c  uo
            STC(4,IQ,IB) = OBS(IQQ,1)
c  vo
            STC(5,IQ,IB) = OBS(IR,1)
c  uf*uo + vf*vo
            STC(6,IQ,IB)=BAK(IQQ,1,IB)*OBS(IQQ,1)+BAK(IR,1,IB)*OBS(IR,1)
c  uf**2 + vf**2
            STC(7,IQ,IB) = BAK(IQQ,1,IB)**2 + BAK(IR,1,IB)**2
c  uo**2 + vo**2
            STC(8,IQ,IB) = OBS(IQQ,1)**2 + OBS(IR,1)**2
c  sqrt (uf**2 + vf**2) - sqrt (uo**2 + vo**2)
            STC(9,IQ,IB) = SQRT(STC(7,IQ,IB)) - SQRT(STC(8,IQ,IB))
C
         ENDIF
C
         if(idbug.eq.1) then
         if(j.eq.3) print *,' ib ',ib,' iq ',iq,' stc ',
     *              (stc(k,iq,ib),k=1,3)
         endif
c
 2234    continue
c
c... end variable-loop
         ENDDO
c... end forecast-loop
         ENDDO
 
         M = ITYP(SUBSET)
 
         if (m > 0) then
         DO N=1   ,NBAK
         DO LL=1  ,NREG
C
         IF(REGION(XOB,YOB,LL)) THEN
C
            DO K=1,NVAR
            nstat=6
            if(k.eq.nvar) nstat=9

c        print *,'j is ',j,' k is ',k,' ll is ',ll,' m is ',m,' n is ',n
            cnto = SPRS(1,J,K,LL,M,N) 
            SPRS(1,J,K,LL,M,N) = SPRS(1,J,K,LL,M,N) + STC(1,K,N)
            cntn = SPRS(1,J,K,LL,M,N) 

            if(cntn.gt.cnto) then
            wt1 = cnto/cntn
            wt2 = 1.-wt1
c
            DO I=2,nstat
c
            sprso = SPRS(I,J,K,LL,M,N)
            rat1 = wt1*SPRS(I,J,K,LL,M,N)
            rat2 = wt2*STC(I,K,N)
c
            SPRS(I,J,K,LL,M,N) = rat1 + rat2
            sprsn = SPRS(I,J,K,LL,M,N)
c
            if((idbug.eq.1).and.(i.eq.2)) then
            if((j.eq.4).and.(k.eq.2)) then
            if((m.eq.1).and.(n.eq.1)) then
            write(6,3000) ll,cnto,cntn,wt1,wt2,sprso,STC(I,K,N),
     *      rat1,rat2,sprsn
            endif
            endif
            endif
c
            ENDDO          ! ... end stat-loop
c
            endif

c... end variable-loop
            ENDDO
C
         ENDIF
c... end region-loop
         ENDDO
C
c... end forecast-loop
         ENDDO
         endif     ! if itype > 0
C
c...  only if correct qc flag
      ENDIF
c...  only if correct level
      ENDIF
 
      ENDDO
c... end subset-loop
      ENDDO
 
       CALL CLOSBF(LUBFR)

C  FINISH UP
C  ---------
C   write out grads data file...
 
      iw=50
      do ibak=1,nbak
      iw=iw+1
c
      do ivarx=1,nvar
      ivar=indexv(ivarx)
      nstat=6
      if(ivar.eq.nvar) nstat=9
      if(idbug.eq.1) print *,ibak,' ivar ',ivar,' nstat ',nstat
c
      do nst=1,nstat
c
      do iplv=1,nplv
c
      do isub=1,nsub
      do ireg=1,nreg
c
      gdata(ireg,isub)=sprs(nst,iplv,ivar,ireg,isub,ibak)
c
c... end region-loop
      enddo
c... end subset-loop
      enddo
c
      write(iw) gdata
c
      if(ipr.eq.1) then
      if(ibak.eq.1) then
      if(iplv.eq.1) then
      ilevt=pmandt(iplv)
      ilevb=pmandb(iplv)
      if(ivar.eq.1)
     *write(6,1231) ibak,nst,ilevt,ilevb,(gdata(ireg,1),ireg=1,nreg)
      if(ivar.eq.2)
     *write(6,1232) ibak,nst,ilevt,ilevb,(gdata(ireg,1),ireg=1,nreg)
      if(ivar.eq.3)
     *write(6,1233) ibak,nst,ilevt,ilevb,(gdata(ireg,1),ireg=1,nreg)
      if(ivar.eq.4)
     *write(6,1234) ibak,nst,ilevt,ilevb,(gdata(ireg,1),ireg=1,nreg)
      endif
      endif
      endif
c
c... end level-loop
      enddo
c... end stat-loop
      enddo
c
c... end variable-loop
      enddo
c
      close(iw)
c... end forecast-loop
      enddo
c
 1231  format('q fcs=',i2,2x,'stat = ',i2,2x,'lev = ',2i6,2x,7f12.2)
 1232  format('t fcs=',i2,2x,'stat = ',i2,2x,'lev = ',2i6,2x,7f12.2)
 1233  format('z fcs=',i2,2x,'stat = ',i2,2x,'lev = ',2i6,2x,7f12.2)
 1234  format('w fcs=',i2,2x,'stat = ',i2,2x,'lev = ',2i6,2x,7f12.2)
c
 3000  format('reg ',i2,' cnto ',f4.0,' cntn ',f4.0,
     *        ' wt1 ',f5.2,' wt2 ',f5.2,2x,5f12.2)
c
      PRINT'("SURUFIT PROCESSING COMPLETED")'
      STOP
      END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      LOGICAL FUNCTION REGION(X,Y,I)
 
      PARAMETER(NREG=7)
 
      DIMENSION X1(NREG),X2(NREG),Y1(NREG),Y2(NREG)
 
C.... X goes from east to west
C.... Y goes from south to north
      DATA X1(1),X2(1),Y1(1),Y2(1) /0.,360.,-90.,90./
      DATA X1(2),X2(2),Y1(2),Y2(2) /0.,360.,20.,80./
      DATA X1(3),X2(3),Y1(3),Y2(3) /0.,360.,-80.,-20./
      DATA X1(4),X2(4),Y1(4),Y2(4) /0.,360.,-20.,20./
      DATA X1(5),X2(5),Y1(5),Y2(5) /235.,295.,25.,55./
      DATA X1(6),X2(6),Y1(6),Y2(6) /350.,25.,35.,70./
      DATA X1(7),X2(7),Y1(7),Y2(7) /65.,145.,5.,45./

C-----------------------------------------------------------------------
 
C  CHECK FOR A VALID REGION INDEX AND BE OPTIMISTIC
C  ------------------------------------------------
 
      IF(I.LT.0.OR.I.GT.NREG) THEN
         REGION = .FALSE.
         RETURN
      ELSE
         REGION = .TRUE.
      ENDIF
 
C  SETUP THE SEARCH PARAMETERS
C  ---------------------------
 
      IF(I.EQ.0) THEN
         I1 = 1
         I2 = NREG
      ELSE
         I1 = I
         I2 = I
      ENDIF
 
C  LOOK FOR A REGION MATCH
C  -----------------------
 
      DO I0=I1,I2
      IF(Y.GE.Y1(I0) .AND. Y.LE.Y2(I0)) THEN
         IF(X1(I0).LE.X2(I0)) THEN
            IF(X.GE.X1(I0) .AND. X.LE.X2(I0)) RETURN
         ELSEIF(X1(I0).GT.X2(I0)) THEN
            IF(X.GE.X1(I0) .OR.  X.LE.X2(I0)) RETURN
         ENDIF
      ENDIF
      ENDDO
 
C  IF NO MATCH, RETURN FALSE
C  -------------------------
 
      REGION = .FALSE.
      RETURN
      END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      INTEGER FUNCTION ITYP(SUBSET)
 
      PARAMETER(NSUB=1)
 
      CHARACTER*8 SUBSET,SUBTYP(1)
      DATA SUBTYP /'AIRCFT'/
C
C     CHARACTER*8 SUBSET,SUBTYP(15)
C     DATA SUBTYP /'ADPUPA','AIRCAR','AIRCFT','SATWND','PROFLR',
C    .             'VADWND','SATBOG','SATEMP','ADPSFC','SFCSHP',
C    .             'SFCBOG','SPSSMI','SYNDAT','ERS1DA','GOESND'/
 
C  LOOK FOR A MATCH TO RETURN NON-ZERO
C  -----------------------------------
 
      DO I=1,NSUB
      ITYP = I
      IF(SUBSET.EQ.SUBTYP(I)) RETURN
      ENDDO
 
C  IF NO MATCH, RETURN ZERO
C  ------------------------
 
      ITYP = 0
      RETURN
      END