***********************************************************************
* PROGRAM: PRODUCE_RPSS
* LANGUAGE: FORTRAN 77
*---------------------------------------------------------------------
* PURPOSE: Calculates Ranked Probability Skill Scores
*---------------------------------------------------------------------
* USAGE: ./rpss.Linux.exe ${1} ${2} ${3} ${4} ${5}
*     where ${1} - record number corresponding to current valid period
*                  being scored
*           ${2} - Filename of Binary file containing gridded 
*                  probabilistic below normal probability values
*           ${3} - Filename of Binary file containing gridded 
*                  probabilistic above normal probability values
*           ${4} - Filename of Binary file containing gridded 
*                  observation categories
*           ${5} - Flag indicating observation source [H=Wei Shi Temp,
*                  Pingping Precip; U=Dave Unger Temperature]
*----------------------------------------------------------------*
* INPUT FILES:
*    mask_2x2.data - Land/Ocean mask for 2x2 degree CONUS grid
*    AFILE - Binary file containing gridded probabilistic above 
*            normal probability values
*    BFILE - Binary file containing gridded probabilistic below 
*            normal probability values
*    OFILE - Binary file containing gridded observation categories
* OUTPUT FILES: NONE (output written to screen/standard output)
*----------------------------------------------------------------*
* SUBROUTINES USED:
*    NONE
* FUNCTIONS USED:
*    NONE
*----------------------------------------------------------------*
* INPUT VARIABLES:
*    AFILE - Filename of Binary file containing gridded
*            probabilistic above normal probability values
*    BFILE - Filename of Binary file containing gridded
*            probabilistic below normal probability values
*    MREC - record number corresponding to current valid period
*           being scored
*    OBSTYPE - Flag indicating observation source [H=Wei Shi Temp,
*              Pingping Precip; U=Dave Unger Temperature]
*    OFILE - Filename of Binary file containing gridded
*            observation categories
* LOCAL VARIABLES:
*    C - Matrix containing Cumulative forecast probability by category
*    CL - Counter for class (i.e. below, near normal, above normal)
*    CLASS - Number of possible forecast/observation categories
*    CO - Matrix containing Cumulative observation probability by category
*    COT - Matrix containing Cumulative observation probability
*    CPER - Matrix containing Cumulative climatological probability by category
*    CPERT - Matrix containing Cumulative climatological probability
*    CT - Matrix containing Cumulative forecast probability 
*    DOMAIN - Flag indicating whether gird point falls within forecast
*             domain
*    FABV - gridded probabilistic above normal probability values
*    FBLW - gridded probabilistic below normal probability values
*    FORECASTS - Flag indicating whether a valid forecast exists for
*                a particular grid point
*    GRID_POINTS - Number of valid grid points
*    I - Counter for grid box location (X direction)
*    IREC - record number corresponding to current valid period
*           being scored
*    J - Counter for grid box location (Y direction)
*    OBS - gridded observation category values
*    OBSERVATIONS - Flag indicating whether valid observations exist for
*                   a particular grid point
*    OCL - Counter for observed category
*    P - Matrix containing forecast probability values for all 
*        three categories
*    PO - Matrix containing observation probability values for all
*        three categories
*    RPSC - Matrix containing climatological Ranked Probability Score
*           for each category and grid box
*    RPSF - Matrix containing forecast Ranked Probability Score
*           for each category and grid box
*    RPSCSUM - Array containing climatological Ranked Probability Score 
*              for all categories and for each category individually 
*    RPSFSUM - Array containing forecast Ranked Probability Score
*              for all categories and for each category individually
*    RPSS - Array containing Ranked Probability Skill Scores for
*           all categories as well as each category individually
*    XMASK - Matrix containing land/ocean mask data
*    XX - Number of grid boxes in X direction
*    YY - Number of grid boxes in Y direction
*----------------------------------------------------------------*
* AUTHOR(S): SCOTT HANDEL
*----------------------------------------------------------------*
* DATE: NOVEMBER 09, 2012
*----------------------------------------------------------------*
* MODIFICATIONS:
******************************************************************
      PROGRAM PRODUCE_RPSS
      IMPLICIT NONE
      CHARACTER*1 OBSTYPE
      CHARACTER*4 MREC
      CHARACTER*300 OFILE,BFILE,AFILE
      INTEGER I,J,IREC,XX,YY,CLASS
      INTEGER CL,OCL,GRID_POINTS
      PARAMETER(XX=36)
      PARAMETER(YY=19)
      PARAMETER(CLASS=3)
      REAL XMASK(XX,YY)
      REAL FBLW(XX,YY),FABV(XX,YY),OBS(XX,YY)
      REAL CPER(XX,YY,CLASS),CPERT(XX,YY),C(XX,YY,CLASS),CO(XX,YY,CLASS)
      REAL CT(XX,YY),COT(XX,YY),RPSS(0:CLASS)
      REAL RPSFSUM(0:CLASS),RPSCSUM(0:CLASS)
      REAL RPSC(XX,YY,CLASS),RPSF(XX,YY,CLASS)
      REAL P(XX,YY,CLASS),PO(XX,YY,CLASS)
      LOGICAL FORECASTS,OBSERVATIONS,DOMAIN
***********************************************************************
* GET COMMAND LINE ARGUMENTS
***********************************************************************
      CALL GETARG (1,MREC)
      CALL GETARG (2,BFILE)
      CALL GETARG (3,AFILE)
      CALL GETARG (4,OFILE)
      CALL GETARG (5,OBSTYPE)
***********************************************************************
* OPEN INPUT FILES
***********************************************************************
      OPEN (40,FILE='../library/mask_2x2.data',ACCESS='DIRECT',
     X  FORM='UNFORMATTED',RECL=XX*YY*4,STATUS='OLD')
      OPEN (51,FILE=BFILE,ACCESS='DIRECT',
     X  FORM='UNFORMATTED',RECL=XX*YY*4,STATUS='OLD')
      OPEN (52,FILE=AFILE,ACCESS='DIRECT',
     X  FORM='UNFORMATTED',RECL=XX*YY*4,STATUS='OLD')
      OPEN (60,FILE=OFILE,
     X       ACCESS='DIRECT',FORM='UNFORMATTED',RECL=XX*YY*4,
     X       STATUS='UNKNOWN')
***********************************************************************
* READ 2X2 MASK DATA
***********************************************************************
      READ(40,REC=1)((XMASK(I,J),I=1,XX),J=1,YY)
***********************************************************************
* CONVERT RECORD NUMBER TO INTEGER
***********************************************************************
      READ(MREC,'(I4)') IREC
***********************************************************************
* READ FORECAST DATA
***********************************************************************
      READ(51,REC=IREC)FBLW
      READ(52,REC=IREC)FABV
      IF (OBSTYPE .EQ. 'H') THEN
         READ(60,REC=IREC)OBS
      ELSEIF (OBSTYPE .EQ. 'U') THEN
         READ(60,REC=1)OBS
      ENDIF
***********************************************************************
* TRANSFORM FORECAST PROBABILITY MATRIX
***********************************************************************
      DO I=1,XX
         DO J=1,YY
            P(I,J,1)=FBLW(I,J)
            P(I,J,3)=FABV(I,J)
            P(I,J,2)=100.-(FBLW(I,J)+FABV(I,J))
         ENDDO
      ENDDO
*****************************************************************
*     ESTABLISH OBS PROBABILITIES                               *
*****************************************************************
      DO CL=1,CLASS
         IF (OBSTYPE .EQ. 'H') THEN
            OCL=CL-2
         ELSEIF (OBSTYPE .EQ. 'U') THEN
            OCL=CL
         ENDIF
         DO I=1,XX
            DO J=1,YY
               IF (OBS(I,J) .EQ. OCL) THEN
                  PO(I,J,CL)=1.0
               ELSE
                  PO(I,J,CL)=0.0
               ENDIF
            ENDDO
         ENDDO
      ENDDO
***********************************************************************
* INITIALIZE RPSS PARAMETERS
***********************************************************************
      DO CL=1,CLASS
         DO I=1,XX
            DO J=1,YY
               RPSF(I,J,CL)=0.0
               RPSC(I,J,CL)=0.0
            ENDDO
         ENDDO
         RPSFSUM(CL)=0.0
         RPSCSUM(CL)=0.0
         RPSS(CL)=-99.00
      ENDDO
***********************************************************************
* INITIALIZE NUMBER OF VALID GRID POINTS
***********************************************************************
      GRID_POINTS=0
***********************************************************************
* LOOP THROUGH GRID POINTS
***********************************************************************
      DO I=1,XX
         DO J=1,YY
***********************************************************************
* CHECK TO SEE IF GRID POINT FALLS WITHIN LAND/OCEAN MASK
***********************************************************************
            IF (XMASK(I,J) .LT. 0.0) THEN
               DOMAIN=.FALSE.
            ELSE
               DOMAIN=.TRUE.
            ENDIF
***********************************************************************
* CHECK TO SEE IF OBSERVATION DATA EXISTS FOR A PARTICULAR GRID POINT
***********************************************************************
            IF (OBSTYPE .EQ. 'H') THEN
               IF (OBS(I,J) .LT. -1.0 .OR. OBS(I,J) .GT. 1.0) THEN
                  OBSERVATIONS=.FALSE.
               ELSE
                  OBSERVATIONS=.TRUE.
               ENDIF
            ELSEIF (OBSTYPE .EQ. 'U') THEN
               IF (OBS(I,J) .LT. 1.0 .OR. OBS(I,J) .GT. 3.0) THEN
                  OBSERVATIONS=.FALSE.
               ELSE
                  OBSERVATIONS=.TRUE.
               ENDIF
            ENDIF
***********************************************************************
* CHECK TO SEE IF FORECAST DATA EXISTS FOR A PARTICULAR GRID POINT
***********************************************************************
            FORECASTS=.TRUE.
            DO CL=1,CLASS
               IF (P(I,J,CL) .LT. 0.0 .OR. P(I,J,CL) .GT. 100.0)
     $         FORECASTS=.FALSE.
            ENDDO
***********************************************************************
* BEGIN CALCULATIONS FOR VALID GRID POINTS
***********************************************************************
            IF (DOMAIN .AND. OBSERVATIONS .AND. FORECASTS) THEN
***********************************************************************
* COUNT NUMBER OF VALID GRID POINTS
***********************************************************************
               GRID_POINTS=GRID_POINTS+1
*****************************************************************
* CONVERT PDFS TO CDFS                                          *
*****************************************************************
               CT(I,J)=0
               COT(I,J)=0
               CPERT(I,J)=0
               DO CL=1,CLASS
                  CT(I,J)=CT(I,J)+P(I,J,CL)
                  COT(I,J)=COT(I,J)+PO(I,J,CL)
                  CPERT(I,J)=CPERT(I,J)+0.3333
                  C(I,J,CL)=CT(I,J)
                  CO(I,J,CL)=COT(I,J)
                  CPER(I,J,CL)=CPERT(I,J)
               ENDDO
*****************************************************************
* DO FCST AND CLIMO RPS FOR RPSS SCORES                         *
*****************************************************************
               DO CL=1,CLASS
                  RPSF(I,J,CL)=((C(I,J,CL)/100.0)-CO(I,J,CL))**2
                  RPSC(I,J,CL)=(CPER(I,J,CL)-CO(I,J,CL))**2
                  RPSFSUM(0)=RPSFSUM(0)+RPSF(I,J,CL)
                  RPSFSUM(CL)=RPSFSUM(CL)+RPSF(I,J,CL)
                  RPSCSUM(0)=RPSCSUM(0)+RPSC(I,J,CL)
                  RPSCSUM(CL)=RPSCSUM(CL)+RPSC(I,J,CL)
               ENDDO
            ENDIF
         ENDDO
      ENDDO
*****************************************************************
* CALCULATE RPSS                                                *
*****************************************************************
      IF ( GRID_POINTS .GT. 0 ) THEN
         DO CL=0,CLASS
            IF ( RPSFSUM(CL) .LE. 0.0 .OR. RPSCSUM(CL) .LE. 0.0) THEN
               RPSS(CL)=-99.00
            ELSE
               RPSS(CL)=1-(RPSFSUM(CL)/RPSCSUM(CL))
            ENDIF
         ENDDO
***********************************************************************
* WRITE RPSS
***********************************************************************
         WRITE(*,35) RPSS(0),GRID_POINTS
35       FORMAT(F8.2,10x,I5)
***********************************************************************
* WRITE ERROR MESSAGE IF NO VALID GRID POINTS FOUND
***********************************************************************
      ELSE
         WRITE(*,*) 'No valid grid points found'
      ENDIF
************************************************************************
* CLOSE FILES
************************************************************************
      CLOSE(40)
      CLOSE(50)
      CLOSE(51)
      CLOSE(52)
      CLOSE(60)
      END
