***********************************************************************
* PROGRAM: PRODUCE_HEIDKE
* LANGUAGE: FORTRAN 77
*---------------------------------------------------------------------
* PURPOSE: Calculates Heidke Skill Scores
*---------------------------------------------------------------------
* USAGE: ./heidke.Linux.exe ${1} ${2} ${3} ${4}
*     where ${1} - record number corresponding to current valid period
*                  being scored
*           ${2} - Filename of Binary file containing gridded 
*                  categorical forecast values (full path)
*           ${3} - Filename of Binary file containing gridded 
*                  observation categories
*           ${4} - 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
*    FFILE - Binary file containing gridded categorical forecast 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:
*    FFILE - Filename of Binary file containing gridded categorical
*            forecast values (full path)
*    MREC - record number corresponding to current valid period
*           being scored
*    OFILE - Filename of Binary file containing gridded observation
*            categories 
*    OBSTYPE - Flag indicating observation source [H=Wei Shi Temp,
*              Pingping Precip; U=Dave Unger Temperature]
* LOCAL VARIABLES:
*    ALLCOR - Number of forecasts correct (including EC) 
*    ALLEXP - Expected number of forecasts correct (including EC)
*    ALLHEI - Heidke skill score of all forecasts (including EC)
*    CVG - Percent coverage
*    EXPECT - Expected number of non-EC forecasts correct
*    EXPNOF - Expected number of EC forecasts correct
*    HEIDKE - Heidke skill score of non-EC forecasts
*    I - Counter for grid X location
*    ICOUNT - Total number of forecasts
*    IREC - record number corresponding to current valid period
*           being scored
*    IRIGHT - Total number of correct forecasts
*    ITOT - Total number of non-EC forecasts
*    J - Counter for grid Y location
*    KREC - record number corresponding to current valid period
*           being scored
*    NOFCST - Number of EC forecasts
*    T - Matrix containing categorical observation values
*    XMASK - Matrix of gridded mask data
*    Y - Matrix contaning categorical forecast values
*----------------------------------------------------------------*
* AUTHOR(S): SCOTT HANDEL
*----------------------------------------------------------------*
* DATE: NOVEMBER 09, 2012
*----------------------------------------------------------------*
* MODIFICATIONS:
******************************************************************
       PROGRAM PRODUCE_HEIDKE
       IMPLICIT NONE
       CHARACTER*1 OBSTYPE
       CHARACTER*4 MREC
       CHARACTER*300 FFILE,OFILE
       INTEGER NOFCST,ICOUNT,IRIGHT,ITOT,IREC,KREC,I,J
       REAL ALLHEI,ALLEXP,ALLCOR,EXPNOF,HEIDKE,CVG,EXPECT
       REAL T(36,19)
       REAL Y(36,19),XMASK(36,19)
***********************************************************************
* INITIALIZE DATA MATRICIES
***********************************************************************
       DATA Y/684*-99./,T/684*-99./
***********************************************************************
*      GET COMMAND LINE ARGUMENTS
***********************************************************************
       CALL GETARG (1,MREC)
       CALL GETARG (2,FFILE)
       CALL GETARG (3,OFILE)
       CALL GETARG (4,OBSTYPE)
***********************************************************************
* OPEN INPUT FILES
***********************************************************************
       OPEN (40,FILE='../library/mask_2x2.data',ACCESS='DIRECT',
     X   FORM='UNFORMATTED',RECL=36*19*4,STATUS='OLD')
       OPEN (50,FILE=FFILE,ACCESS='DIRECT',
     X   FORM='UNFORMATTED',RECL=36*19*4,STATUS='OLD')
       OPEN (60,FILE=OFILE,
     X       ACCESS='DIRECT',FORM='UNFORMATTED',RECL=36*19*4,
     X       STATUS='UNKNOWN')
***********************************************************************
* READ FORECAST MASK
***********************************************************************
       READ(40,REC=1)((XMASK(I,J),I=1,36),J=1,19)
***********************************************************************
* READ RECORD NUMBER
***********************************************************************
       READ(MREC,'(I4)') IREC
***********************************************************************
* DETERMINE CORRECT OBSERVATION NUMBER BASED ON WHICH OBS YOU ARE
* READING IN (RECORD NUMBER IS 1 FOR UNGER'S OBS AND IREC FOR HALPERTS
* OBS
***********************************************************************
       IF (OBSTYPE .EQ. 'H' ) THEN
          KREC=IREC
       ELSE
          KREC=1
       ENDIF 
***********************************************************************
* READ FORECASTS AND OBS
***********************************************************************
       READ(50,REC=IREC)Y
       READ(60,REC=KREC)T
***********************************************************************
*    COMPARE FORECAST TO VERIFICATION FILE
***********************************************************************
        ITOT=0
        IRIGHT=0
        ICOUNT=0
        DO 100 I=1,36
           DO 100 J=1,19
              IF(XMASK(I,J).LT.0.0) GO TO 100
              IF(T(I,J).LT.-50.) GO TO 100
              ICOUNT=ICOUNT+1
              IF(Y(I,J).EQ.0.0) GO TO 100
              IF ( OBSTYPE .EQ. 'H' ) THEN
                 IF((Y(I,J)-2).EQ.T(I,J)) IRIGHT=IRIGHT+1
              ELSE
                 IF(Y(I,J).EQ.T(I,J)) IRIGHT=IRIGHT+1
              ENDIF
              ITOT=ITOT+1 
100     CONTINUE
***********************************************************************
* CALCULATE NON-EC HEIDKE
***********************************************************************
        EXPECT=FLOAT(ITOT)/3.0
        CVG=(FLOAT(ITOT)/FLOAT(ICOUNT))*100.0
        HEIDKE = ((FLOAT(IRIGHT)-EXPECT)/(FLOAT(ITOT)-EXPECT))*100.0
***********************************************************************
* CALCULATE ALL DOMAIN HEIDKE
***********************************************************************
        NOFCST=ICOUNT-ITOT
        EXPNOF=FLOAT(NOFCST)/3.0
        ALLCOR=FLOAT(IRIGHT)+EXPNOF
        ALLEXP=FLOAT(ICOUNT)/3.0
        ALLHEI=((ALLCOR-ALLEXP)/(FLOAT(ICOUNT)-ALLEXP))*100.
***********************************************************************
* WRITE OUT HEIDKE SCORES
***********************************************************************
        WRITE(*,2001)HEIDKE,ALLHEI,CVG
2001    FORMAT(F7.2,F7.2,F7.2)
        RETURN
        END
