!--------------------------------------------------------------------
! Program to read DATEM format measured and calculated air 
! concentration files and perform some elementary statistical 
! analyses for model comparison. Procedures based upon ETEX analyses.
!--------------------------------------------------------------------
! Unit 10 - measurement data
! Unit 20 - model calculation results
! Unit 30 - statistical results summary
! Unit 40 - optional merged measured and calculated output file
! Unit 50 - optional station selection file for statistics
! Unit 60 - supplemental appended output file
! Unit 90 - diagnostic message file
!--------------------------------------------------------------------
! Last Revised: 19 Jul 2001 - original version
!               02 Apr 2002 - change in rank calculation
!               16 Dec 2003 - more flexible format
!               14 May 2004 - included optional station file
!               06 Jul 2004 - expanded zero testing
!               13 Sep 2006 - minor code change
!               12 Jun 2008 - simplified directories
!               09 Feb 2009 - included FAR, TS, POD statistics
!               25 Sep 2009 - simplified directories to command line
!               31 Mar 2010 - correction to model string length
!               22 Aug 2011 - spatial averaging station -1 now period 
!               27 Sep 2011 - skip identical concentrations for KSP test
!               04 Nov 2011 - supplemental appended output file
!               11 Dec 2011 - supplemental file name change
!               18 Oct 2012 - bootstrap resampling for correlation
!               31 Oct 2012 - index correction in bootstrap resampling
!--------------------------------------------------------------------

PROGRAM STATMAIN

  IMPLICIT NONE

  INTEGER, ALLOCATABLE :: SELECT(:)             ! station select list
  INTEGER, ALLOCATABLE :: SITE(:)               ! station label
  INTEGER, ALLOCATABLE :: IYR(:),IMO(:),IDA(:)  ! sample date
  INTEGER, ALLOCATABLE :: IHR(:),IDUR(:)        ! time & duration 
  REAL,    ALLOCATABLE :: SLAT(:),SLON(:)       ! location
  REAL,    ALLOCATABLE :: CALC(:), MEAS(:)      ! concentrations
  REAL,    ALLOCATABLE :: CALCW(:),MEASW(:)     ! working array 

  LOGICAL :: SIGR,ZEXCLD,ZMODEL,ZPLUME,BOOT
  REAL    :: RANK,FMS,FMT
  REAL    :: FRAC,BIAS,CINT,CFACT
  REAL    :: RSUM,RCNT,RDEL,RPCT
  REAL    :: MBAR,CBAR,CORL,TEQV,RATIO,NMSE
  REAL    :: FOEX,FA2,FA5,FAP,FZC,FZM,FAR,POD,TS,CLVL
  REAL    :: ZERO,KSP,MSORT(99),CSORT(99),MM,CC
  INTEGER :: K,N,NUMB,KAVG,KRET,KOUT,KOBS,NUMOBS,NUMCAL,NUMAVG,NUMSTN

  INTEGER       :: DATE_TIME(8)
  CHARACTER(12) :: REAL_CLOCK(3)

  REAL          :: OLAT,OLON
  LOGICAL       :: FTEST
  CHARACTER(20) :: MFILE
  COMMON /INFO/    FTEST,OLAT,OLON,MFILE

!--------------------------------------------------------------------
! subroutine interface definitions
  INTERFACE            
    SUBROUTINE DLOAD (KUNIT,SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON, &
                      MEAS,CALC,NUMSTN,SELECT,CFACT)
    IMPLICIT NONE
    INTEGER, INTENT(IN)    :: KUNIT
    INTEGER, INTENT(INOUT) :: SITE(:)             
    INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:)
    INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       
    REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)     
    REAL,    INTENT(INOUT) :: MEAS(:),CALC(:) 
    INTEGER, INTENT(IN)    :: NUMSTN
    INTEGER, INTENT(IN)    :: SELECT(:)
    REAL,    INTENT(IN)    :: CFACT    
    END SUBROUTINE DLOAD

    SUBROUTINE CUMDIST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,CSORT,MSORT)
    IMPLICIT NONE
    REAL,      INTENT(IN)  :: MEAS(:),CALC(:)
    REAL,      INTENT(IN)  :: ZERO
    LOGICAL,   INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME
    REAL,      INTENT(OUT) :: CSORT(99),MSORT(99) 
    END SUBROUTINE CUMDIST

    SUBROUTINE AVRGS (SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NEWOBS)
    IMPLICIT NONE
    INTEGER, INTENT(INOUT) :: SITE(:)              ! station label
    INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:) ! sample date
    INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       ! time duration 
    REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)      ! location
    REAL,    INTENT(INOUT) :: MEAS(:),CALC(:)      ! concentration
    INTEGER, INTENT(OUT)   :: NEWOBS               ! revised count
    END SUBROUTINE AVRGS

    SUBROUTINE AVRGT (SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NEWOBS)
    IMPLICIT NONE
    INTEGER, INTENT(INOUT) :: SITE(:)              ! station label
    INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:) ! sample date
    INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       ! time duration 
    REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)      ! location
    REAL,    INTENT(INOUT) :: MEAS(:),CALC(:)      ! concentration
    INTEGER, INTENT(OUT)   :: NEWOBS               ! revised count
    END SUBROUTINE AVRGT

    SUBROUTINE KSPARM (CSORT,MSORT,KSP)
    IMPLICIT NONE
    REAL,     INTENT(IN)  :: CSORT(99),MSORT(99)
    REAL,     INTENT(OUT) :: KSP
    END SUBROUTINE KSPARM

    SUBROUTINE FIGMST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,FMS)
    IMPLICIT NONE
    REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
    REAL,    INTENT(IN)  :: ZERO 
    LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME  
    REAL,    INTENT(OUT) :: FMS
    END SUBROUTINE FIGMST

    SUBROUTINE SPATIAL (MEAS,CALC,CLVL,FAR,POD,TS)
    IMPLICIT NONE
    REAL,    INTENT(IN)  :: MEAS(:),CALC(:),CLVL
    REAL,    INTENT(OUT) :: FAR,POD,TS
    END SUBROUTINE SPATIAL


    SUBROUTINE BIASES (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,BIAS,CINT)
    IMPLICIT NONE
    REAL, INTENT(IN)    :: MEAS(:),CALC(:)
    REAL, INTENT(IN)    :: ZERO
    LOGICAL, INTENT(IN) :: ZEXCLD,ZMODEL,ZPLUME
    REAL, INTENT(OUT)   :: BIAS,CINT
    END SUBROUTINE BIASES

    SUBROUTINE SCATER (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,  &
                       FOEX,FA2,FA5,FAP,FZC,FZM)
    IMPLICIT NONE
    REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
    REAL,    INTENT(IN)  :: ZERO 
    LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME  
    REAL,    INTENT(OUT) :: FOEX,FA2,FA5,FAP,FZC,FZM
    END SUBROUTINE SCATER

    SUBROUTINE CZERO(VALUE,ZERO)
    IMPLICIT NONE
    REAL,   INTENT(IN)    :: VALUE(:)
    REAL,   INTENT(INOUT) :: ZERO    
    END SUBROUTINE CZERO

    SUBROUTINE LINEAR (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,MBAR,CBAR,  &
                       NUMB,CORL,NMSE,SIGR,TEQV)
    IMPLICIT NONE
    REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
    REAL,    INTENT(IN)  :: ZERO
    LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME
    INTEGER, INTENT(OUT) :: NUMB   
    REAL,    INTENT(OUT) :: MBAR,CBAR,CORL,NMSE,TEQV
    LOGICAL, INTENT(OUT) :: SIGR
    END SUBROUTINE LINEAR

    SUBROUTINE SETOBS (kunit,numobs,numstn,select)
    IMPLICIT NONE
    INTEGER, INTENT(IN)    :: KUNIT
    INTEGER, INTENT(OUT)   :: NUMOBS
    INTEGER, INTENT(IN)    :: NUMSTN
    INTEGER, INTENT(INOUT) :: SELECT(:) 
    END SUBROUTINE SETOBS

    SUBROUTINE SETARG (ZERO,KAVG,KOUT,ZEXCLD,ZMODEL,ZPLUME,BOOT,NUMSTN,CFACT,CLVL)
    IMPLICIT NONE
    REAL,    INTENT(OUT) :: ZERO
    INTEGER, INTENT(OUT) :: KAVG
    INTEGER, INTENT(OUT) :: KOUT
    LOGICAL, INTENT(OUT) :: ZEXCLD,ZMODEL,ZPLUME
    LOGICAL, INTENT(OUT) :: BOOT
    INTEGER, INTENT(OUT) :: NUMSTN
    REAL,    INTENT(OUT) :: CFACT
    REAL,    INTENT(OUT) :: CLVL
    END SUBROUTINE SETARG
  END INTERFACE 
!--------------------------------------------------------------------

! diagnostic messages
  OPEN(90,FILE='message.txt')

! read any command line arguments and set default values
  CALL SETARG (ZERO,KAVG,KOUT,ZEXCLD,ZMODEL,ZPLUME,BOOT,NUMSTN,CFACT,CLVL)

! load station list if selection file available
  IF(NUMSTN.GT.0)THEN
     ALLOCATE (select(numstn),STAT=kret)
     CALL SETOBS(50,numobs,numstn,select)
     CLOSE(50)
     write(90,*)'(main) Allocated station select array: ',numstn
  END IF

! open measurement data file to determine array allocations
  IF(KOUT.NE.2)THEN
     CALL SETOBS(10,numobs,numstn,select)
     CALL SETOBS(20,numcal,numstn,select)
  ELSE
     CALL SETOBS(40,numobs,numstn,select)
     NUMCAL=NUMOBS
  END IF
  write(90,*)'(main) Number of measured and calculated: ',numobs,numcal

  ALLOCATE (site (numobs), STAT=kret)
  ALLOCATE (iyr  (numobs), STAT=kret)
  ALLOCATE (imo  (numobs), STAT=kret)
  ALLOCATE (ida  (numobs), STAT=kret)
  ALLOCATE (ihr  (numobs), STAT=kret)
  ALLOCATE (idur (numobs), STAT=kret)
  ALLOCATE (slat (numobs), STAT=kret)
  ALLOCATE (slon (numobs), STAT=kret)
  ALLOCATE (calc (numobs), STAT=kret)
  ALLOCATE (meas (numobs), STAT=kret)
  ALLOCATE (calcw(numobs), STAT=kret)
  ALLOCATE (measw(numobs), STAT=kret)

! set all concentration data to missing
  calc=-1.0
  meas=-1.0

!--------------------------------------------------------------------
! load measurement and model data

  IF(KOUT.NE.2)THEN
     CALL DLOAD (10,site,iyr,imo,ida,ihr,idur,slat,slon,  &
                 meas,calc,numstn,select,cfact)
     write(90,*)'(main) Finished loading measured data'
     CALL DLOAD (20,site,iyr,imo,ida,ihr,idur,slat,slon,  &
                 meas,calc,numstn,select,cfact)
     write(90,*)'(main) Finished loading calculated data'
  ELSE
     CALL DLOAD (40,site,iyr,imo,ida,ihr,idur,slat,slon,  &
                 meas,calc,numstn,select,cfact)
     write(90,*)'(main) Finished loading measured-calculated data'
  END IF

!--------------------------------------------------------------------
! optional temporal or spatial averaging

  IF(KAVG.EQ.0)THEN
!    no averaging
     NUMAVG=NUMOBS
     WRITE(30,'(I10,A)') NUMAVG,' Unaveraged data points for processing'

  ELSEIF(KAVG.EQ.1) THEN
!    temporal averaging
     CALL AVRGT  &
          (SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NUMAVG)
     WRITE(30,'(I10,A)') NUMAVG,' Data points after temporal averaging'

  ELSEIF(KAVG.EQ.2) THEN
!    spatial averaging
     CALL AVRGS  &
          (SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NUMAVG)
     WRITE(30,'(I10,A)') NUMAVG,' Data points after spatial averaging'

  ELSEIF(KAVG.LT.0) THEN
!    temporal low-pass filter
     WRITE(90,*) 'Low-pass filter not yet available'
  END IF

!--------------------------------------------------------------------
! set the zero concentration level

  WRITE(30,'(F10.2,A)') ZERO,' Percentile input for zero measured'
  CALL CZERO (MEAS,ZERO)
  WRITE(30,'(F10.2,A)') ZERO,' Zero measured concentration value'

!--------------------------------------------------------------------
! save measured calculated merged file for other processing programs

  IF(KOUT.EQ.1)THEN
     WRITE(40,'(A5,A8,A7,A8,4A5,2A10)')'Num','  Site  ',  &
          'Lat','Lon','Yr','Mo','Da','Hr','Meas','Calc'

     K=0
     DO N=1,NUMOBS

!       output non-missing pairs (especially applies to averaged data)
        IF(MEAS(N).GE.0.0.AND.CALC(N).GE.0.0)THEN

           MM=MEAS(N)
           CC=CALC(N)
           IF(MM.LE.ZERO)MM=0.0
           IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

!          skip output when both zero
           IF(MM.GT.0.0.OR.CC.GT.0.0)THEN
              K=K+1  
              WRITE(40,'(I5,I8,F7.2,F8.2,4I5,2F10.2)')        &
                    K,SITE(N),SLAT(N),SLON(N),                &             
                    IYR(N),IMO(N),IDA(N),IHR(N),MM,CC
           END IF
        END IF
     END DO
     CLOSE(40)
  END IF

!--------------------------------------------------------------------
! linear averages and correlation

  IF(BOOT)THEN
     RDEL=0.20           ! fraction of data set deleted each pass
     KOBS=SIZE(MEAS,1)   ! number of iterations equals data points
     RSUM=0.0            ! variance summation
     RCNT=0.0            ! summation counter

     CALL DATE_AND_TIME (REAL_CLOCK(1),REAL_CLOCK(2),REAL_CLOCK(3),DATE_TIME)
     RPCT=DATE_TIME(8)/999.0 ! initialize random number with milliseconds

     DO K=1,KOBS
!       reload working data array before each pass through the data
        MEASW=MEAS
        CALCW=CALC
 
        DO N=1,KOBS
!          set RDEL fraction of the observations to missing
           CALL RANDOM_NUMBER(RPCT) 
           IF(RPCT.LT.RDEL)THEN
              MEASW(N)=-1.0
              CALCW(N)=-1.0
           END IF
        END DO
        CALL LINEAR (MEASW,CALCW,ZERO,ZEXCLD,ZMODEL,ZPLUME,MBAR,CBAR,NUMB,CORL,NMSE,SIGR,TEQV)
        RSUM=RSUM+CORL*CORL
        RCNT=RCNT+1.0
     END DO
!    compute linear regression statistics based upon the complete data set
     CALL LINEAR (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,MBAR,CBAR,NUMB,CORL,NMSE,SIGR,TEQV)
!    replace the correlation with the average value from the resampling
     CORL=SQRT(RSUM/RCNT)

  ELSE
     CALL LINEAR (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,MBAR,CBAR,NUMB,CORL,NMSE,SIGR,TEQV)
  END IF

  WRITE(30,'(A)')      ' '
  IF(SIGR)THEN
!    regression significant at the 99% level
     WRITE(30,'(F10.2,A)') CORL,     ' Correlation coefficient (P=99%)'
  ELSE
     WRITE(30,'(F10.2,A)') CORL,     ' Correlation coefficient'
  END IF
  WRITE(30,'(F10.2,A)') TEQV,     ' T-value (|Slope|/Standard Error)'
  WRITE(30,'(F10.2,A)') MBAR,     ' Average measured concentration'
  WRITE(30,'(F10.2,A)') CBAR,     ' Average calculated concentration'
  WRITE(30,'(F10.2,A)') CBAR/MBAR,' Ratio of calculated/measured'
  WRITE(30,'(F10.2,A)') NMSE,     ' Normalized mean square error'
  WRITE(30,'(I10,A)')   NUMB,     ' Number of pairs analyzed'

!--------------------------------------------------------------------
! simple bias and confidence interval

  CALL BIASES (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,BIAS,CINT)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)')  BIAS,       ' Average bias [(C-M)/N]'
  WRITE(30,'(F10.2,A)') (BIAS-CINT), ' Lo 99 % confidence interval'
  WRITE(30,'(F10.2,A)') (BIAS+CINT), ' Hi 99 % confidence interval'
  FRAC=2.0*BIAS/(CBAR+MBAR)
  WRITE(30,'(F10.2,A)') FRAC, ' Fractional bias [2B/(C+M)]'

!--------------------------------------------------------------------
! spatial statistics
  
  CALL SPATIAL (MEAS,CALC,CLVL,FAR,POD,TS)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') FAR, ' False Alarm Rate [fa/(fa+hit)] '
  WRITE(30,'(F10.2,A)') POD, ' Probability of Detection [hit/(hit+miss)] '
  WRITE(30,'(F10.2,A)') TS,  ' Threat Score [hit/(fa+hit+miss)] '
   

!--------------------------------------------------------------------
! figures of merit in space

  CALL FIGMST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,FMS)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') FMS, ' Fig of merit in space (%)'

!--------------------------------------------------------------------
! scatter diagram statistics FOEX, FA5, FA2

  CALL SCATER (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,   &
               FOEX,FA2,FA5,FAP,FZC,FZM)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') FOEX,' Factor exceeding [N(C>M)/N-0.5]'
  WRITE(30,'(F10.2,A)') FA2, ' Percent C/M '//CHAR(177)//' 2'
  WRITE(30,'(F10.2,A)') FA5, ' Percent C/M '//CHAR(177)//' 5'
  WRITE(30,'(F10.2,A)') FAP, ' Percent M>0 and C>0'
  WRITE(30,'(F10.2,A)') FZC, ' Percent M>0 and C=0'
  WRITE(30,'(F10.2,A)') FZM, ' Percent M=0 and C>0'
 
!--------------------------------------------------------------------
! cumulative frequency distribution

  CALL CUMDIST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,CSORT,MSORT)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') MSORT(95),' Measured   95-th percentile'
  WRITE(30,'(F10.2,A)') MSORT(90),' Measured   90-th percentile'
  WRITE(30,'(F10.2,A)') MSORT(75),' Measured   75-th percentile'
  WRITE(30,'(F10.2,A)') MSORT(50),' Measured   50-th percentile'

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') CSORT(95),' Calculated 95-th percentile'
  WRITE(30,'(F10.2,A)') CSORT(90),' Calculated 90-th percentile'
  WRITE(30,'(F10.2,A)') CSORT(75),' Calculated 75-th percentile'
  WRITE(30,'(F10.2,A)') CSORT(50),' Calculated 50-th percentile'

  CALL KSPARM (CSORT,MSORT,KSP)

  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') KSP,' Kolmogorov-Smirnov Parameter'

!--------------------------------------------------------------------
! final rank

  RANK=CORL*CORL+(1.0-ABS(FRAC/2.0))+(FMS/100.0)+(1.0-KSP/100.0) 
  WRITE(30,'(A)')      ' '
  WRITE(30,'(F10.2,A)') RANK,' Final rank (C,FB,FMS,KSP)'

!-------------------------------------------------------------------
! supplemental appended output file

  IF(FTEST) WRITE(60,'(A,1X,8F10.2)') MFILE,OLAT,OLON,CORL,NMSE,FRAC,FMS,KSP,RANK

END PROGRAM STATMAIN


!--------------------------------------------------------------------
! Data filter is used to apply spatial average to the data before 
! statistical analysis. Only temporal variability remains.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               16 Dec 2003
!               22 Aug 2011
!--------------------------------------------------------------------

SUBROUTINE AVRGS(SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NEWOBS)

  IMPLICIT NONE

  INTEGER, INTENT(INOUT) :: SITE(:)              ! station label
  INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:) ! sample date
  INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       ! time duration 
  REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)      ! location
  REAL,    INTENT(INOUT) :: MEAS(:),CALC(:)      ! concentration
  INTEGER, INTENT(OUT)   :: NEWOBS               ! revised count

  INTEGER, ALLOCATABLE   :: KTIME(:)
  INTEGER, ALLOCATABLE   :: KYR(:),KMO(:),KDA(:) 
  INTEGER, ALLOCATABLE   :: KHR(:),KDUR(:)        

  REAL,    ALLOCATABLE   :: MM(:),CC(:),NUM(:)

  INTEGER  :: N,K,KRET,NSITE,NUMOBS,OBTIME 

! number of observation records
  NUMOBS=SIZE(MEAS,1)

! determine number of unique time periods and durations
  ALLOCATE (KTIME(NUMOBS),KDUR(NUMOBS), STAT=KRET)
  ALLOCATE (KYR(NUMOBS),KMO(NUMOBS),KDA(NUMOBS),KHR(NUMOBS),STAT=KRET)

! initialize station counter and clock    
  NSITE=0
  OBTIME=(IYR(1)*365+IMO(1)*31+IDA(1))*24+IHR(1)

! go through each observation
  obs : DO N=1,NUMOBS

!    check list for match with each unique time period & duration
     DO K=1,NSITE
!       simple formula to compute serial time
        OBTIME=(IYR(N)*365+IMO(N)*31+IDA(N))*24+IHR(N)
        IF(OBTIME.EQ.KTIME(K).AND.IDUR(N).EQ.KDUR(K))CYCLE obs
     END DO

!    no match add to list
     NSITE=NSITE+1
     KTIME(NSITE)=OBTIME
     KDUR(NSITE) =IDUR(N)
     KYR(NSITE)  =IYR(N)
     KMO(NSITE)  =IMO(N)
     KDA(NSITE)  =IDA(N)
     KHR(NSITE)  =IHR(N)
  END DO obs

! averaging variables
  ALLOCATE (MM(NSITE),CC(NSITE),NUM(NSITE), STAT=KRET)
  MM = 0.0
  CC = 0.0
  NUM= 0.0

! go through each observation
  sum : DO N=1,NUMOBS

!    skip missing data
     IF(MEAS(N).LT.0.0.OR.CALC(N).LT.0.0)CYCLE sum

!    check list for match
     stn : DO K=1,NSITE
        OBTIME=(IYR(N)*365+IMO(N)*31+IDA(N))*24+IHR(N)
        IF(OBTIME.EQ.KTIME(K).AND.IDUR(N).EQ.KDUR(K))EXIT stn
     END DO stn

!    accumulate for that time
     MM(K) =MM(K)+MEAS(N)
     CC(K) =CC(K)+CALC(N)
     NUM(K)=NUM(K)+1.0
  END DO sum

! set all data to missing 
  SITE = -1
  MEAS = -1.0
  CALC = -1.0
  SLAT = -1.0
  SLON = -1.0

! and then fill with average
  N=0
  avg : DO K=1,NSITE
     IF(NUM(K).LE.0.0)CYCLE avg
     N=N+1
     SITE(N)=N
     MEAS(N)=MM(K)/NUM(K)
     CALC(N)=CC(K)/NUM(K)
     IDUR(N)=KDUR(K)
     IYR(N) =KYR(K)
     IMO(N) =KMO(K)
     IDA(N) =KDA(K)
     IHR(N) =KHR(K)
  END DO avg
  NEWOBS=N

  DEALLOCATE (KTIME,KDUR,KYR,KMO,KDA,KHR,MM,CC,NUM)

END SUBROUTINE AVRGS


!--------------------------------------------------------------------
! Data filter is used to apply temporal average to the data before 
! statistical analysis. Only spatial variability remains.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               16 Dec 2003
!--------------------------------------------------------------------

SUBROUTINE AVRGT(SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,MEAS,CALC,NEWOBS)

  IMPLICIT NONE

  INTEGER, INTENT(INOUT) :: SITE(:)              ! station label
  INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:) ! sample date
  INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       ! time duration 
  REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)      ! location
  REAL,    INTENT(INOUT) :: MEAS(:),CALC(:)      ! concentration
  INTEGER, INTENT(OUT)   :: NEWOBS               ! revised count

  INTEGER, ALLOCATABLE   :: KSITE(:)
  REAL,    ALLOCATABLE   :: MM(:),CC(:),NUM(:)
  REAL,    ALLOCATABLE   :: TLAT(:),TLON(:)

  INTEGER  :: N,K,KRET,NSITE,NUMOBS 

! number of observation records
  NUMOBS=SIZE(MEAS,1)

! determine number of unique stations
  ALLOCATE (KSITE(NUMOBS), STAT=KRET)
  ALLOCATE (TLAT(NUMOBS),TLON(NUMOBS), STAT=KRET)
     
  NSITE=0
! go through each observation
  obs : DO N=1,NUMOBS

!    check list for match, building up the number of entries
!    with each new data observation
     DO K=1,NSITE
        IF(KSITE(K).EQ.SITE(N))CYCLE obs
     END DO

!    no match add to list
     NSITE=NSITE+1
     KSITE(NSITE)=SITE(N)
     TLAT (NSITE)=SLAT(N)
     TLON (NSITE)=SLON(N)

  END DO obs

! averaging variables
  ALLOCATE (MM(NSITE),CC(NSITE),NUM(NSITE), STAT=KRET)
  MM = 0.0
  CC = 0.0
  NUM= 0.0

! go through each observation
  sum : DO N=1,NUMOBS

!    skip missing data
     IF(MEAS(N).LT.0.0.OR.CALC(N).LT.0.0)CYCLE sum

!    check list and exit loop with match
     stn : DO K=1,NSITE
        IF(KSITE(K).EQ.SITE(N)) EXIT stn
     END DO stn

!    accumulate for that station
     MM(K) =MM(K)+MEAS(N)
     CC(K) =CC(K)+CALC(N)
     NUM(K)=NUM(K)+1.0

  END DO sum

! set all time data to missing 
  MEAS = -1.0
  CALC = -1.0
  IDUR = -1
  IYR  = -1
  IMO  = -1
  IDA  = -1
  IHR  = -1

! fill the concentration data arrays with the average
  N=0
  avg : DO K=1,NSITE
     IF(NUM(K).LE.0.0) CYCLE avg
     N=N+1
     MEAS(N)=MM(K)/NUM(K)
     CALC(N)=CC(K)/NUM(K)
     SITE(N)=KSITE(K)
     SLAT(N)=TLAT(K)
     SLON(N)=TLON(K)
  END DO avg
  NEWOBS=N

  DEALLOCATE (KSITE,TLAT,TLON,MM,CC,NUM)

END SUBROUTINE AVRGT


!--------------------------------------------------------------------
! Computes simple bias and confidence interval at the 5% level using
! tabulated student-t test values. 
!--------------------------------------------------------------------
! Last Revised: 02 Aug 2001
!               06 Jul 2004 - expanded zero testing
!--------------------------------------------------------------------

SUBROUTINE BIASES (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,BIAS,CINT)

  IMPLICIT NONE

  REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
  REAL,    INTENT(IN)  :: ZERO  
  LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME
  REAL,    INTENT(OUT) :: BIAS,CINT  ! returns bias and 
                                     ! confidence interval

  LOGICAL :: pick
  INTEGER :: k,numobs
  REAL    :: mm,cc,summ,sumk

! 99% confidence interval using single sided student-t test (P=0.005)
  REAL    :: tval(8)

! for degrees:   2     4     8    16    32    64   128   inf
  DATA TVAL / 9.92, 4.60, 3.36, 2.92, 2.75, 2.66, 2.62, 2.58/

! number of observation records
  NUMOBS=SIZE(MEAS,1)

!-----------------------------------
! simple bias 

  SUMM=0.0
  SUMK=0.0

  DO K=1,NUMOBS

     CC=CALC(K)
     MM=MEAS(K)

!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN

!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.PICK).AND.ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))PICK=.TRUE.

        IF(PICK)THEN
!          summation of differences
           SUMM=SUMM+(CC-MM)
           SUMK=SUMK+1.0
        END IF

     END IF
  END DO

  BIAS=0.0
  IF(SUMK.GT.0.0) BIAS=SUMM/SUMK

!--------------------------------------
! confidence interval 

  SUMM=0.0
  SUMK=0.0

  DO K=1,NUMOBS

     CC=CALC(K)
     MM=MEAS(K)

!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN

!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF(ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0)) PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))     PICK=.TRUE.

        IF(PICK)THEN
!          summation of square of difference from mean
           SUMM=SUMM+(CC-MM-BIAS)**2
           SUMK=SUMK+1.0
        END IF

     END IF

  END DO

! find nearest t-value from table to compute confidence interval
  IF(SUMK.GT.0.0)THEN
    K=MAX(1,MIN(8,NINT(ALOG(SUMK)/ALOG(2.0))))
    CINT=TVAL(K)*SQRT(SUMM/(SUMK*(SUMK-1.0)))
  ELSE
    CINT=0.0
  END IF

END SUBROUTINE BIASES


!---------------------------------------------------------------------
! Compute cumulative concentration distribution independently for
! measured and calculated for those pairs when both are non-zero.
!---------------------------------------------------------------------
! last revised: 02 Aug 2001
!               06 Jul 2004 - expanded zero testing
!---------------------------------------------------------------------

SUBROUTINE CUMDIST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,CSORT,MSORT)

  IMPLICIT NONE

  REAL,      INTENT(IN)  :: MEAS(:),CALC(:)
  REAL,      INTENT(IN)  :: ZERO
  LOGICAL,   INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME
  REAL,      INTENT(OUT) :: CSORT(99),MSORT(99) 

  REAL,      ALLOCATABLE :: TEMP(:)  ! temp working array
  INTEGER :: NUMVAL,NUMOBS,KRET,K
  LOGICAL :: PICK
  REAL    :: CC,MM

! subroutine interface definitions
  INTERFACE
    SUBROUTINE SORT (NUMB,VALUE,DSORT)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: NUMB
    REAL,    INTENT(IN)  :: VALUE(:)
    REAL,    INTENT(OUT) :: DSORT(99)
    END SUBROUTINE SORT
  END INTERFACE

! number of observation records
  NUMOBS=SIZE(MEAS,1)

! number of non-zero pairs
  ALLOCATE (TEMP(NUMOBS), STAT=KRET)

! measured distribution

  NUMVAL=0
  DO K=1,NUMOBS
     CC=CALC(K)
     MM=MEAS(K)
!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN
!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.PICK).AND.ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))PICK=.TRUE.

        IF(PICK)THEN
           NUMVAL=NUMVAL+1            
           TEMP(NUMVAL)=MM
        END IF
     END IF
  END DO
  CALL SORT(NUMVAL,TEMP,MSORT)
  
! calculated distribution

  NUMVAL=0
  DO K=1,NUMOBS
     CC=CALC(K)
     MM=MEAS(K)
!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN
!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF(ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0)) PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))     PICK=.TRUE.

        IF(PICK)THEN
           NUMVAL=NUMVAL+1            
           TEMP(NUMVAL)=CC
        END IF
     END IF
  END DO
  CALL SORT(NUMVAL,TEMP,CSORT)
  
  DEALLOCATE (TEMP)

END SUBROUTINE CUMDIST


!---------------------------------------------------------------------
! Convert percentile zero concentration level specified on the command
! line to actual value. Zero is input as cumulative percentage level 
! of all non-zero concentrations. The percentage is converted to air 
! concentration units at that percentile after ascending order data 
! sort. 
!---------------------------------------------------------------------
! last revised: 18 Jul 2001
!---------------------------------------------------------------------

SUBROUTINE CZERO(VALUE,ZERO)

  IMPLICIT NONE

  REAL,   INTENT(IN)    :: VALUE(:)    ! input data array
  REAL,   INTENT(INOUT) :: ZERO        ! in percent; out concentration

  REAL,     ALLOCATABLE :: TEMP(:)     ! temp working array

  INTEGER :: N,KRET,KMIN,NUMALL,NUMVAL ! counters
  REAL    :: DSORT(99)                 ! sortted output percentiles

! subroutine interface definitions
  INTERFACE
    SUBROUTINE SORT (NUMB,VALUE,DSORT)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: NUMB
    REAL,    INTENT(IN)  :: VALUE(:)
    REAL,    INTENT(OUT) :: DSORT(99)
    END SUBROUTINE SORT
  END INTERFACE

! test for true zero
  IF(ZERO.EQ.0.0)RETURN

! input array size
  NUMALL=SIZE(VALUE,1)

! determine number of non-zero concentrations
  NUMVAL=0
  DO N=1,NUMALL
     IF(VALUE(N).GT.0.0)NUMVAL=NUMVAL+1
  END DO
  ALLOCATE (TEMP(NUMVAL), STAT=KRET)

! copy non-zero values into dummy array 
  NUMVAL=0
  DO N=1,NUMALL
     IF(VALUE(N).GT.0.0)THEN
        NUMVAL=NUMVAL+1
        TEMP(NUMVAL)=VALUE(N)
     END IF
  END DO
  
! ascending order sort
  CALL SORT(NUMVAL,TEMP,DSORT)
  write(90,*)'(czero) 95th percentile concentration: ',dsort(95)
  write(90,*)'(czero) 90th percentile concentration: ',dsort(90)
  write(90,*)'(czero) 75th percentile concentration: ',dsort(75)
  write(90,*)'(czero) 50th percentile concentration: ',dsort(50)

! find concentration of selected percentile
  ZERO=DSORT(INT(ZERO))

  DEALLOCATE (TEMP)

END SUBROUTINE CZERO


!--------------------------------------------------------------------
! Load measurement followed by model data according to unit number. 
! On the first pass the measurement data are loaded and on the 
! second pass the model data are loaded to match each measurement 
! data record. All data loaded are saved in memory.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               16 Dec 2003
!               17 May 2004 - station selection option
!--------------------------------------------------------------------

SUBROUTINE DLOAD (KUNIT,SITE,IYR,IMO,IDA,IHR,IDUR,SLAT,SLON,  &
                  MEAS,CALC,NUMSTN,SELECT,CFACT)

  IMPLICIT NONE

  INTEGER, INTENT(IN)    :: KUNIT
  INTEGER, INTENT(INOUT) :: SITE(:)              ! station label
  INTEGER, INTENT(INOUT) :: IYR(:),IMO(:),IDA(:) ! sample date
  INTEGER, INTENT(INOUT) :: IHR(:),IDUR(:)       ! time duration 
  REAL,    INTENT(INOUT) :: SLAT(:),SLON(:)      ! location
  REAL,    INTENT(INOUT) :: MEAS(:),CALC(:)      ! concentration
  INTEGER, INTENT(IN)    :: NUMSTN               ! number of stations
  INTEGER, INTENT(IN)    :: SELECT(:)            ! station listing
  REAL,    INTENT(IN)    :: CFACT                ! normalization factor

  LOGICAL  :: PICK
  INTEGER  :: KSITE,KRET
  INTEGER  :: N,K,KYR,KMO,KDA,KHR,KDUR,NUMOBS 
  REAL     :: TLAT,TLON,TCON,XCON

  CHARACTER(80) :: LABEL

! number of observation records
  NUMOBS=SIZE(SITE,1)

! skip first two identification records
  READ(KUNIT,'(A80)')LABEL
  READ(KUNIT,'(A80)')
  write(90,*)'(dload) Processing file: ',label(:40)
  write(90,*)'(dload) On unit number: ',kunit

  n=0
  kret=0

  DO WHILE (kret.EQ.0)

     IF(KUNIT.NE.40)THEN
        READ(KUNIT,*,IOSTAT=kret) KYR,KMO,KDA,KHR,KDUR,TLAT,TLON,TCON,KSITE

     ELSE
!       reads file of measured-calculated pairs created by this program
        READ(KUNIT,'(I5,I8,F7.2,F8.2,4I5,2F10.2)',IOSTAT=kret)   &
             k,ksite,tlat,tlon,kyr,kmo,kda,khr,tcon,xcon
        IF(N.LT.NUMOBS.AND.KRET.EQ.0)THEN
           n=n+1
           iyr(n)=kyr
           imo(n)=kmo
           ida(n)=kda
           ihr(n)=khr
           idur(n)=9999
           slat(n)=tlat
           slon(n)=tlon
           meas(n)=tcon*cfact
           calc(n)=xcon*cfact
           site(n)=ksite
        ELSEIF(KRET.EQ.0)THEN
           write(90,*)'(dload) Record count exceeding limit: ',n,numobs
        END IF
     END IF

!    check for station match 
     IF(numstn.GT.0.AND.kret.EQ.0)THEN
        PICK=.FALSE.
        stnloop : DO k=1,numstn
           IF(select(k).EQ.ksite)THEN
              PICK=.TRUE.
              EXIT stnloop
           END IF
        END DO stnloop

     ELSEIF(KRET.NE.0)THEN
        PICK=.FALSE.
     ELSE
        PICK=.TRUE.
     END IF

     IF(PICK.AND.KUNIT.EQ.10)THEN
        N=N+1
        IF(N.GT.NUMOBS)THEN
           WRITE(*,*)'Measurement data count exceeding array: ',numobs
           STOP
        END IF

!       measurement data always loaded first
        iyr(n)=kyr
        imo(n)=kmo
        ida(n)=kda
        ihr(n)=khr
        idur(n)=kdur
        slat(n)=tlat
        slon(n)=tlon
        meas(n)=tcon*cfact
        site(n)=ksite

     ELSEIF(PICK.AND.KUNIT.EQ.20)THEN
        DO K=1,NUMOBS
!          match model results to observation record, both start
!          time and sample duration must match
           IF(KYR .EQ.IYR(K) .AND.KMO  .EQ.IMO(K).AND.  &
              KDA .EQ.IDA(K) .AND.KHR  .EQ.IHR(K).AND.  &
              KDUR.EQ.IDUR(K).AND.KSITE.EQ.SITE(K))THEN
              n=n+1
              calc(k)=tcon*cfact
           END IF
        END DO
     END IF

  END DO
  write(90,*)'(dload) Loaded # data points: ',n

  CLOSE(KUNIT)

END SUBROUTINE DLOAD


!--------------------------------------------------------------------
! Figures of Merit in Space (FMS) is the ratio of area covered by 
! both calculated and measured patterns divided by the total 
! calculated and measured area (Intersection over Union). Values are
! not contoured to determine area, but as a simplification each 
! station is assumed to cover the same are as all other stations and
! hence station counts are used for the calculation. 
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               06 Jul 2004 - additional zero tests
!--------------------------------------------------------------------

SUBROUTINE FIGMST (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,FMS)

  IMPLICIT NONE

  REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
  REAL,    INTENT(IN)  :: ZERO
  LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME   
  REAL,    INTENT(OUT) :: FMS  

  LOGICAL :: pick
  INTEGER :: k,numobs
  REAL    :: mm,cc,sumi,sumu,summax,summin

! number of observation records
  NUMOBS=SIZE(MEAS,1)

!---------------------------------
! figure of merit in space

  SUMU=0.0        ! union count
  SUMI=0.0        ! intersection count

  DO K=1,NUMOBS

     CC=CALC(K)
     MM=MEAS(K)

!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN

!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0
 
        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.PICK).AND.ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))PICK=.TRUE.

        IF(PICK)THEN
!          number in union
           IF(MM.GT.0.0.OR. CC.GT.0.0) SUMU=SUMU+1.0
!          number in intersection
           IF(MM.GT.0.0.AND.CC.GT.0.0) SUMI=SUMI+1.0
        END IF

     END IF
  END DO

  FMS=0.0
  IF(SUMU.GT.0.0) FMS=100.0*SUMI/SUMU
 
END SUBROUTINE FIGMST


!--------------------------------------------------------------------
! False Alarm Rate (FAR), Probability of Detection (POD), and 
! Threat Score (TS) Calculation
!--------------------------------------------------------------------
! Last Revised: 09 Feb 2009
!--------------------------------------------------------------------

SUBROUTINE SPATIAL (MEAS,CALC,CLVL,FAR,POD,TS)

  IMPLICIT NONE

  REAL,    INTENT(IN)  :: MEAS(:),CALC(:),CLVL
  REAL,    INTENT(OUT) :: FAR,POD,TS  

  INTEGER :: K,NUMOBS,HIT,MISS,FALSE
  REAL    :: MM,CC

! number of observation records
  NUMOBS=SIZE(MEAS,1)

! initialize values
  HIT = 0
  MISS = 0
  FALSE = 0

  DO K=1,NUMOBS

     CC=CALC(K)
     MM=MEAS(K)

     IF(MM .GE. CLVL .AND. CC .GE. CLVL)HIT = HIT + 1
     IF(MM .GE. CLVL .AND. CC .LT. CLVL)MISS = MISS + 1
     IF(MM .LT. CLVL .AND. CC .GE. CLVL)FALSE = FALSE + 1

  END DO
  
  IF(FLOAT(FALSE)+FLOAT(HIT).EQ.0.0)THEN
     FAR=0.0
  ELSE
     FAR = 100.*(FLOAT(FALSE) / (FLOAT(FALSE) + FLOAT(HIT)))  
  END IF
  TS  = 100.*(FLOAT(HIT) / (FLOAT(FALSE) + FLOAT(HIT) + FLOAT(MISS)))
  POD = 100.*(FLOAT(HIT) / (FLOAT(HIT) + FLOAT(MISS)))

END SUBROUTINE SPATIAL


!--------------------------------------------------------------------
! Kolmogorov-Smirnov parameter is defined as the absolute value of
! the maximum difference between two cumulative distributions. Values
! are returned as the percent difference. 
!-------------------------------------------------------------------- 
! Last Revised: 11 Jul 2001
!               11 May 2011 - changed range from 1-99% to 5-95%
!               27 Sep 2011 - skip identical concentrations
!--------------------------------------------------------------------

SUBROUTINE KSPARM (CSORT,MSORT,KSP)

  IMPLICIT NONE

  REAL,     INTENT(IN)  :: CSORT(99),MSORT(99)
  REAL,     INTENT(OUT) :: KSP

  INTEGER   :: KM,KC

  KSP=0.0

! find probabilities for each concentration pair by going through
! each distribution, previously sorted into ascending order, until
! the concentrations match.

  DO KM=5,95

     IF(MSORT(KM).NE.MSORT(KM-1))THEN
     cloop : DO KC=5,95
        IF(CSORT(KC).EQ.CSORT(KC-1)) CYCLE cloop
        IF(CSORT(KC).GE.MSORT(KM))THEN
           KSP=MAX(KSP,ABS(FLOAT(KC-KM)))
           EXIT cloop
        END IF
     END DO cloop
     END IF

  END DO

END SUBROUTINE KSPARM


!--------------------------------------------------------------------
! Linear statistics for mean values, correlation, and normalized mean
! square error.
!--------------------------------------------------------------------
! Last Revised: 02 Aug 2001
!               06 Jul 2004 - expanded zero testing
!               28 Sep 2011 - test for SER=0 when computing T-value
!--------------------------------------------------------------------

SUBROUTINE LINEAR (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,MBAR,CBAR,   &
                   NUMB,CORL,NMSE,SIGR,TEQV)

  IMPLICIT NONE

  REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
  REAL,    INTENT(IN)  :: ZERO
  LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME
  INTEGER, INTENT(OUT) :: NUMB   
  REAL,    INTENT(OUT) :: MBAR,CBAR,CORL,NMSE,TEQV
  LOGICAL, INTENT(OUT) :: SIGR

  LOGICAL :: pick
  INTEGER :: k,numobs
  REAL    :: sxy,sxx,syy,var,ser,slope
  REAL    :: mm,cc,summ,sumc,sumk,sqrm,sqrc,sqrb,sqrd

! 99% confidence interval using single sided student-t test (P=0.01)
  REAL    :: tval(8)

! for degrees:   2     4     8    16    32    64   128   inf
  DATA TVAL / 6.96, 3.75, 2.90, 2.58, 2.46, 2.39, 2.36, 2.33/

! number of observation records
  NUMOBS=SIZE(MEAS,1)

! standard sums
  SUMM=0.0
  SUMC=0.0
  SUMK=0.0

! squared sums
  SQRM=0.0
  SQRC=0.0
  SQRB=0.0
  SQRD=0.0

  DO K=1,NUMOBS

     CC=CALC(K)
     MM=MEAS(K)

!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN
             
!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.PICK).AND.ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))PICK=.TRUE.

        IF(PICK)THEN
           SUMM=SUMM+MM
           SUMC=SUMC+CC
           SUMK=SUMK+1.0
           SQRM=SQRM+MM*MM
           SQRC=SQRC+CC*CC
           SQRB=SQRB+MM*CC
           SQRD=SQRD+(CC-MM)**2
         END IF

     END IF

  END DO

  SIGR=.FALSE.
  IF(SUMK.GT.2.0)THEN
!    mean values
     MBAR=SUMM/SUMK
     CBAR=SUMC/SUMK
     NUMB=INT(SUMK)

!    normalized mean square error
     IF(MBAR.EQ.0.0.OR.CBAR.EQ.0.0)THEN
        NMSE=999.99
     ELSE
        NMSE=SQRD/MBAR/CBAR/SUMK
     END IF

!    variance about the regression
     SYY=SQRC-SUMC*SUMC/SUMK
     SXX=SQRM-SUMM*SUMM/SUMK
     SXY=SQRB-SUMC*SUMM/SUMK
     VAR=SYY-SXY*SXY/SXX

!    slope and correlation coefficient
     SLOPE=SXY/SXX
     IF(SXX.EQ.0.0.OR.SYY.EQ.0.0)THEN
        CORL=0.0
     ELSE
        CORL=SXY/SQRT(SXX*SYY)  
     END IF

!    standard error of the regression coefficient
     SER=SQRT(VAR/SXX/(SUMK-2.0))

!    equivalent t value
     TEQV=999.99
     IF(SER.NE.0.0)TEQV=ABS(SLOPE)/SER

!    significant if TEQV greater than t-value from table
     K=MAX(1,MIN(8,NINT(ALOG(SUMK-2.0)/ALOG(2.0))))
     IF(TEQV.GT.TVAL(K))SIGR=.TRUE.
 
  ELSE
     CORL=0.0
     NMSE=0.0
     MBAR=0.0
     CBAR=0.0
     NUMB=0
  END IF
  
END SUBROUTINE LINEAR


!--------------------------------------------------------------------
! Scatter diagram statistics to compute factor exceeding, within a
! factor of two and five, and number of zero calculated and measured.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               06 Jul 2004 - additional zero options
!--------------------------------------------------------------------

SUBROUTINE SCATER (MEAS,CALC,ZERO,ZEXCLD,ZMODEL,ZPLUME,  &
                   FOEX,FA2,FA5,FAP,FZC,FZM)

  IMPLICIT NONE

  REAL,    INTENT(IN)  :: MEAS(:),CALC(:)
  REAL,    INTENT(IN)  :: ZERO 
  LOGICAL, INTENT(IN)  :: ZEXCLD,ZMODEL,ZPLUME  
  REAL,    INTENT(OUT) :: FOEX,FA2,FA5,FAP,FZC,FZM

  LOGICAL :: pick
  INTEGER :: k,numb,nump,numx,num2,num5,nzc,nzm,numobs
  REAL    :: mm,cc,ratio

  NUMOBS=SIZE(MEAS,1)   ! number of observation records
  NUMX=0                ! number of calculated exceeding measured
  NUM2=0                ! number within a factor of 2
  NUM5=0                ! number within a factor of 5
  NUMP=0                ! number of plume cases
  NUMB=0                ! number of non-zero pairs
  NZC=0                 ! number of zero calculated
  NZM=0                 ! number of zero measured

  DO K=1,NUMOBS
     CC=CALC(K)
     MM=MEAS(K)

!    both values zero or positive (not missing)        
     IF(MM.GE.0.0.AND.CC.GE.0.0)THEN

!       adjust for zero
        IF(MM.LE.ZERO)MM=0.0
        IF(ZMODEL.AND.CC.LE.ZERO)CC=0.0

        PICK=.FALSE.
        IF(ZPLUME.AND.(MM.GT.0.0.AND.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.PICK).AND.ZEXCLD.AND.(MM.GT.0.0.OR.CC.GT.0.0))PICK=.TRUE.
        IF((.NOT.ZEXCLD).AND.(.NOT.ZPLUME))PICK=.TRUE.

!       skip statistics if both values = 0.0
        IF(PICK)THEN

!          total number over/under prediction
           IF(CC.GT.MM)NUMX=NUMX+1

!          number of non 0,0 pairs
           NUMB=NUMB+1

           IF(MM.EQ.0.0)THEN
!             measured zero
              NZM=NZM+1

           ELSEIF(CC.EQ.0.0)THEN
!             calculated zero
              NZC=NZC+1

           ELSE
!             ratios within a factor of 2 or factor of 5
              RATIO=CC/MM
              IF(RATIO.GE.0.5.AND.RATIO.LE.2.0)THEN
                 NUM2=NUM2+1
                 NUM5=NUM5+1
                 NUMP=NUMP+1
              ELSEIF(RATIO.GE.0.2.AND.RATIO.LE.5.0)THEN
                 NUM5=NUM5+1
                 NUMP=NUMP+1
              ELSE
                 NUMP=NUMP+1
              END IF
           END IF

        END IF

     END IF
  END DO

! scatter results as percentage over/under, factor-2, factor-5
  FOEX= 100.0*NUMX/NUMB-50.0
  FA2 = 100.0*NUM2/NUMB
  FA5 = 100.0*NUM5/NUMB
  FAP = 100.0*NUMP/NUMB
  FZC = 100.0*NZC /NUMB
  FZM = 100.0*NZM /NUMB

END SUBROUTINE SCATER


!--------------------------------------------------------------------
! Process command line arguments and open all files for either input
! or output.
!--------------------------------------------------------------------
! Last Revised: 03 Aug 2001 - original version
!               28 May 2002 - change in default model variation
!               18 May 2004 - optional station selection file
!               06 Jul 2004 - more zero options
!               05 Apr 2005 - different output file using merged data 
!               09 Feb 2009 - optional contingency level for spatial
!                             statistics
!               03 Nov 2011 - supplemental appended output file
!               17 Oct 2012 - added bootstrap resampling option
!--------------------------------------------------------------------

SUBROUTINE SETARG (ZERO,KAVG,KOUT,ZEXCLD,ZMODEL,ZPLUME,BOOT,NUMSTN,CFACT,CLVL)

  IMPLICIT NONE

! percentile level in the measured data that is used to define 
! the threshold air concentration for statistics
  REAL,    INTENT(OUT) :: ZERO  
! indicator for temporal (=1) or spatial (=2) averaging                                
  INTEGER, INTENT(OUT) :: KAVG
! indicator to create merged output file                                
  INTEGER, INTENT(OUT) :: KOUT
! flag to exclude zero-zero pairs from statistics
  LOGICAL, INTENT(OUT) :: ZEXCLD, ZMODEL, ZPLUME
! flag for bootstrap resampling
  LOGICAL, INTENT(OUT) :: BOOT
! number of stations to analyze
  INTEGER, INTENT(OUT) :: NUMSTN
! concentration normalization factor
  REAL,    INTENT(OUT) :: CFACT
! contingency level for spatial statistics
  REAL,    INTENT(OUT) :: CLVL

  CHARACTER(80) :: DNAME,RNAME,SNAME,LABEL,MODEL,SITES,SUFFIX,TITLE
  CHARACTER(1)  :: TRACER
  LOGICAL       :: FTEST
  INTEGER       :: KRET,KLEN,KLEN0,KLEN1,KLEN2,KLEN3,KLEN4,KLEN5,KARG,IARGC

  KARG=IARGC()
  IF(KARG.LT.1)THEN
     WRITE(*,*)'Program to read DATEM format measured and calculated air' 
     WRITE(*,*)'concentration files and perform some elementary statistical' 
     WRITE(*,*)'analyses for model comparison. Procedures based upon ETEX.'
     WRITE(*,*)'Input and output file names are generated automatically when'
     WRITE(*,*)'the tracer character ID is set to a non-zero value.'
     WRITE(*,*)
     WRITE(*,*)'USAGE: statmain [arguments]'
     WRITE(*,*)'-a[averaging: space (s), time (t), low-pass filter (#)]'
     WRITE(*,*)'-b[bootstrap resampling to compute correlation]'
     WRITE(*,*)'-c[concentration normalization factor]'
     WRITE(*,*)'-d[measurement directory or file name (when -t0)]'
     WRITE(*,*)'-f[station selection file suffix]'
     WRITE(*,*)'-l[contingency level for spatial statistics (1)]'
     WRITE(*,*)'-m[model variation string (when -t<>0]'
     WRITE(*,*)'-o[write (1) merged data or read (2) merged file]'
     WRITE(*,*)'-p[plume statistics (both meas & calc gt zero)]'
     WRITE(*,*)'-r[calculation results directory or file name (when -t0)]'
     WRITE(*,*)'-s[supplemental appended output file name]'
     WRITE(*,*)'-t[tracer character number: (0),1,2,3, etc]'
     WRITE(*,*)'-x[exclude 0,0 pairs]'
     WRITE(*,*)'-y[set calculated to zero when below measured zero]'
     WRITE(*,*)'-z[percentile level of non-zero for zero measured]'
     STOP
  END IF

! concentration normalization factor
  CFACT=1.0

! contingency level for spatial statistics
  CLVL = 1.

! default values if not set on command line
  ZERO=0.0

! tracer identification
  TRACER='0'

! experimental measurement data
  KLEN0=0
  DNAME=''

! results directory
  KLEN1=0
  RNAME=''

! optional model variations string
  KLEN2=0
  MODEL=''

! station selection file
  KLEN3=0
  SITES=''

! output file label suffix
  KLEN4=0
  SUFFIX=''

! supplemental appended output file
  KLEN5=0
  SNAME=''

! bootstrap resampling
  BOOT=.FALSE.

  KAVG=0
  KOUT=0

  ZEXCLD=.FALSE.
  ZMODEL=.FALSE.
  ZPLUME=.FALSE.
 
! go through each argument
  DO WHILE (KARG.GT.0)

     CALL GETARG(KARG,LABEL)

     SELECT CASE (LABEL(1:2))

!    experiment directory or measurement data file
     CASE ('-d','-D')
        KLEN0=INDEX(LABEL,' ')-1
        DNAME=LABEL(3:KLEN0)
        KLEN0=INDEX(DNAME,' ')-1
        write(90,*)'(setarg) Set experiment directory: ',dname(:klen0)

!    results directory or results file
     CASE ('-r','-R')
        KLEN1=INDEX(LABEL,' ')-1
        RNAME=LABEL(3:KLEN1)
        KLEN1=INDEX(RNAME,' ')-1
        write(90,*)'(setarg) Set results directory: ',rname(:klen1)

!    model variation or input file
     CASE ('-m','-M')
        KLEN2=INDEX(LABEL,' ')-1
        MODEL=LABEL(3:KLEN2)
        KLEN2=INDEX(MODEL,' ')-1
        write(90,*)'(setarg) Set model variation: ',model(:klen2)

!    station selection file
     CASE ('-f','-F')
        KLEN3=INDEX(LABEL,' ')-1
        SITES=LABEL(3:KLEN3)
        KLEN3=INDEX(SITES,' ')-1
        write(90,*)'(setarg) Set station selection: ',sites(:klen3)

!    supplemental appended output file
     CASE ('-s','-S')
        KLEN5=INDEX(LABEL,' ')-1
        SNAME=LABEL(3:KLEN5)
        KLEN5=INDEX(SNAME,' ')-1
        write(90,*)'(setarg) Set supplemental appended file: ',sname(:klen5)

!    tracer number
     CASE ('-t','-T')
        TRACER=LABEL(3:3)
        write(90,*)'(setarg) Set tracer number: ',tracer

!    percentile level to define zero measured concentrations
     CASE ('-z','-Z')
        READ(LABEL(3:),'(F10.0)')ZERO
        ZERO=MAX(0.0,MIN(99.0,ZERO))
        write(90,*)'(setarg) Set tracer zero percentile: ',zero

!    temporal or spatial averaging
     CASE ('-a','-A')
        SELECT CASE (LABEL(3:3))
        CASE ('t','T')
           KAVG=1
           write(90,*)'(setarg) Set temporal average'
        CASE ('s','S')
           KAVG=2
           write(90,*)'(setarg) Set temporal average'
        CASE DEFAULT
           READ(LABEL(3:),'(I5)')KAVG
           write(90,*)'(setarg) Set temporal filter',KAVG
           KAVG=-KAVG
        END SELECT

!    merged data output file
     CASE ('-o','-O')
        SELECT CASE (LABEL(3:3))
        CASE ('1')
           KOUT=1
        CASE ('2')
           KOUT=2
        CASE DEFAULT
           KOUT=1
        END SELECT

        IF(KOUT.EQ.1)THEN
           write(90,*)'(setarg) Set file for output of measured-calculated pairs'
        ELSE
           write(90,*)'(setarg) Set file for input of measured-calculated pairs'
        END IF

!    exclude 0,0 pairs
     CASE ('-x','-X')
        ZEXCLD=.TRUE.
        write(90,*)'(setarg) Exclude 0,0 pairs from statistics'

!    set measured zero to value other than zero
     CASE ('-y','-Y')
        ZMODEL=.TRUE.
        write(90,*)'(setarg) Set calculated zero same as measured zero'

!    include only plumes (meas and calc gt zero) in statistics
     CASE ('-p','-P')
        ZPLUME=.TRUE.
        write(90,*)'(setarg) Compute only plume statistics'

!    concentration normalization factor
     CASE ('-c','-C')
        READ(LABEL(3:),'(F10.0)')CFACT
        write(90,*)'(setarg) Set concentration factor: ',cfact

!    contingency level is used as a threshold for detection statistics  
     CASE ('-l','-L')
        READ(LABEL(3:),'(F10.3)')CLVL
        write(90,*)'(setarg) Set contingency level for spatial: ',CLVL

!    turn on bootstrap resampling for correlation
     CASE ('-b','-B')
        BOOT=.TRUE.
        write(90,*)'(setarg) Turn on bootstrap resampling for correlation'

     END SELECT
     KARG=KARG-1
  END DO

! consistency checks
  IF(ZEXCLD.AND.ZPLUME)ZEXCLD=.FALSE.

  IF(KOUT.EQ.0.OR.KOUT.EQ.1)THEN
!    open the measurement data file
     IF(TRACER.EQ.'0')THEN
        LABEL=DNAME(1:KLEN0)
     ELSE
        LABEL=DNAME(1:KLEN0)//'/meas-t'//TRACER//'.txt'
     END IF
     INQUIRE (FILE=LABEL,EXIST=FTEST)
     IF(FTEST)THEN
        OPEN(10,FILE=LABEL,ACTION='READ')
        READ(10,'(A)')TITLE
        BACKSPACE(10)
        write(90,*)'(setarg) Opened measurement file: ',label(:40)
     ELSE
        WRITE(*,*)'File not found:',LABEL
        READ(*,*)
        STOP
     END IF

!    open the model calculations file
     IF(TRACER.EQ.'0')THEN
        LABEL=RNAME(1:KLEN1)
     ELSE
        LABEL=RNAME(1:KLEN1)//'/ct'//TRACER//'_'//MODEL(1:KLEN2)//'.txt'
     END IF
     INQUIRE (FILE=LABEL,EXIST=FTEST)
     IF(FTEST)THEN
        OPEN(20,FILE=LABEL,ACTION='READ')
        write(90,*)'(setarg) Opened calculated file: ',label(:40)
     ELSE
        WRITE(*,*)'File not found:',LABEL
        READ(*,*)
        STOP
     END IF

!    open the supplemental appended output file
     IF(KLEN5.GT.0)THEN
        LABEL=SNAME(1:KLEN5)
        INQUIRE (FILE=LABEL,EXIST=FTEST)
        IF(FTEST)THEN
           OPEN(60,FILE=LABEL,ACTION='READWRITE')
           write(90,*)'(setarg) Opened existing supplemental appended file: ',label(:40)
           KRET=0
           DO WHILE (KRET.EQ.0)
              READ(60,'(A)',IOSTAT=KRET) 
           END DO
        ELSE
           OPEN(60,FILE=LABEL,ACTION='WRITE')
           WRITE(60,'(A20,1X,8A10)') 'File                ','Lat','Lon','Corr',   &
                                     'NMSE','FB','FMS','KSP','Rank'
           write(90,*)'(setarg) Opened new supplemental appended file: ',label(:40)
        END IF
     END IF

  END IF

! open the optional station selection file (except when reading
! the merged data file (kout=2)
  IF(KLEN3.NE.0.AND.KOUT.NE.2)THEN
     IF(TRACER.EQ.'0')THEN
        LABEL=SITES(1:KLEN3)
     ELSE
        LABEL='stn'//TRACER//'_'//SITES(1:KLEN3)//'.txt'
        SUFFIX=SITES
        KLEN4=KLEN3
     END IF

     INQUIRE(FILE=LABEL,EXIST=FTEST)
     IF(FTEST)THEN
        OPEN(50,FILE=LABEL)
        write(90,*)'(setarg) Opened station select file: ',label(:40) 

        READ(50,'(a)')
        READ(50,'(a)')
        numstn=0
        DO WHILE (kret.EQ.0)
           READ(50,'(a)',IOSTAT=kret) 
           numstn=numstn+1
        END DO
        numstn=numstn-1
        REWIND(50)
        write(90,*)'(setarg) Selected # stations: ',numstn
     ELSE
        numstn=0
     END IF

  ELSE
!    no station select file then process data from all sites
     NUMSTN=0
     SITES='All'
     KLEN3=3
  END IF

! define the output file suffix
  IF(TRACER.EQ.'0')THEN
     SUFFIX='.txt' 
     KLEN4=4
  ELSE
     SUFFIX=TRACER//'_'//MODEL(1:KLEN2)//'.txt'
     KLEN4=KLEN2+6
  END IF
 
! open the file for the merged data
  IF(KOUT.EQ.1.OR.KOUT.EQ.2)THEN
     IF(KAVG.EQ.0)THEN
        LABEL='dataA'//SUFFIX(1:KLEN4)
     ELSEIF(KAVG.EQ.1)THEN
        LABEL='dataT'//SUFFIX(1:KLEN4)
     ELSEIF(KAVG.EQ.2)THEN
        LABEL='dataS'//SUFFIX(1:KLEN4)
     END IF
     KLEN=INDEX(LABEL,' ')
     OPEN(40,FILE=LABEL(1:KLEN))
     write(90,*)'(setarg) Opened file for merged data: ',label(1:klen)

     IF(KOUT.EQ.1)WRITE(40,'(A)')TITLE           
  END IF 

! open the statistical results file
  IF(KAVG.EQ.0)THEN
     LABEL='statA'//SUFFIX(1:KLEN4)
  ELSEIF(KAVG.EQ.1)THEN
     LABEL='statT'//SUFFIX(1:KLEN4)
  ELSEIF(KAVG.EQ.2)THEN
     LABEL='statS'//SUFFIX(1:KLEN4)
  END IF

! reading merged data file use special output file name
  IF(KOUT.EQ.2)THEN
     LABEL='statF'//SUFFIX(1:KLEN4)
  END IF
  write(90,*)'(setarg) Opened statistical results file: ',label(1:20)

  KLEN=INDEX(LABEL,' ')
  OPEN(30,FILE=LABEL(1:KLEN))
  WRITE(30,'(A)')' ' 
  WRITE(30,'(11X,A)')DNAME(1:KLEN1)//' results file: '//LABEL(1:KLEN)
  WRITE(30,'(3A)')' Model variation: '//MODEL(1:KLEN2),   &
                '   Tracer number: '//TRACER,             &
                '   Station select: '//SITES(1:KLEN3)
  WRITE(30,'(11X,A)')'----------------------------------- '
  WRITE(30,'(A)')' ' 

END SUBROUTINE setarg


!--------------------------------------------------------------------
! Determine the number of data records in standard format DATEM file
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001 - original version
!               17 May 2004 - added optional station file
!               03 Nov 2011 - check for supplemental information
!--------------------------------------------------------------------

SUBROUTINE SETOBS (kunit,numobs,numstn,select)

  IMPLICIT NONE

  INTEGER, INTENT(IN)    :: KUNIT
  INTEGER, INTENT(OUT)   :: NUMOBS
  INTEGER, INTENT(IN)    :: NUMSTN
  INTEGER, INTENT(INOUT) :: SELECT(:) ! station select list

  INTEGER  :: K,KRET,KSITE
  INTEGER  :: IDUM(5)
  REAL     :: RDUM(3)

  REAL          :: OLAT,OLON
  LOGICAL       :: FTEST
  CHARACTER(20) :: MFILE
  COMMON /INFO/    FTEST,OLAT,OLON,MFILE

  numobs=0
  kret=0

! check for supplemental information in model prediction file
  IF(KUNIT.EQ.20)THEN
     FTEST=.FALSE.
     INQUIRE(UNIT=60,OPENED=FTEST)
     IF(FTEST)THEN
        READ(kunit,'(2F10.0,1X,A)') OLAT,OLON,MFILE
     ELSE
        READ(kunit,'(a)')
     END IF
  ELSE
     READ(kunit,'(a)')
  END IF
  READ(kunit,'(a)')

  DO WHILE (kret.EQ.0)
     IF(KUNIT.EQ.10.OR.KUNIT.EQ.20)THEN
!       load sampling or measurement data
        READ(kunit,*,IOSTAT=kret) IDUM,RDUM,KSITE

     ELSEIF(KUNIT.EQ.50)THEN
!       load station selection data 
        READ(kunit,*,IOSTAT=kret) IDUM(1),KSITE
        IF(KRET.EQ.0) SELECT(numobs+1)=KSITE

     ELSEIF(KUNIT.EQ.40)THEN
!       load data from previously merged (meas & calc) file 
        READ(kunit,*,IOSTAT=kret)
     END IF

!    check for station match
!    numsta>0 only if unit 50 previously read for station listing
     IF(numstn.GT.0.AND.kret.EQ.0)THEN
        stnloop : DO k=1,numstn
           IF(select(k).EQ.ksite)THEN
              numobs=numobs+1
              EXIT stnloop
           END IF
        END DO stnloop

     ELSEIF(kret.EQ.0)THEN
        numobs=numobs+1

     ELSE
        REWIND(kunit)
     END IF

  END DO

END SUBROUTINE setobs


!-------------------------------------------------------------------- 
! Ascending order binary sort which returns the sorted values in
! array DSORT by cumulative percent from 1 to 99.  For instance, 
! DSORT(50) returns the median value, while DSORT(25 & 75) return the
! quartile values. Routine requres 2**N passes through the data, 
! where N equals the number of data points.
!--------------------------------------------------------------------
! Last Revised: 18 Jul 2001
!               14 Aug 2002
!--------------------------------------------------------------------

SUBROUTINE SORT(NUMPTS,VALUE,DSORT)

  IMPLICIT NONE

  INTEGER,  INTENT(IN)  :: NUMPTS
  REAL,     INTENT(IN)  :: VALUE(:)
  REAL,     INTENT(OUT) :: DSORT(99)

  INTEGER, ALLOCATABLE :: TEMP (:)
  INTEGER, ALLOCATABLE :: INDX (:)

  INTEGER :: I,I1,I2,I3,L1,L2,KRET,NLEN,LEN1,LEN2,INDX1,INDX2

  IF(NUMPTS.LE.1)THEN
     WRITE(*,*)'Too few points to sort: ',NUMPTS
     RETURN
  END IF

  ALLOCATE (TEMP(NUMPTS), STAT=KRET)
  ALLOCATE (INDX(NUMPTS), STAT=KRET)

! initialize index array
  DO I=1,NUMPTS
     INDX(I)=I
  END DO

    NLEN=1
100 IF(NLEN.GE.NUMPTS)GO TO 800

    I1=1
    I2=I1+NLEN
    I3=1

200 L1=0
    L2=0
    LEN1=NUMPTS-I1+1
    IF(LEN1.GT.NLEN)LEN1=NLEN
    LEN2=NUMPTS-I2+1
    IF(LEN2.GT.NLEN)LEN2=NLEN

250 INDX1=INDX(I1+L1)
    INDX2=INDX(I2+L2)
    IF(VALUE(INDX1).GT.VALUE(INDX2))GO TO 400

300 TEMP(I3)=INDX(I1+L1)
    L1=L1+1
    I3=I3+1
    GO TO 500

400 CONTINUE
    TEMP(I3)=INDX(I2+L2)
    L2=L2+1
    I3=I3+1

500 IF(I3.GT.NUMPTS)GO TO 550
    IF(L1.LT.LEN1.AND.L2.LT.LEN2)GO TO 250
    IF(L1.EQ.LEN1.AND.L2.LT.LEN2)GO TO 400
    IF(L1.LT.LEN1.AND.L2.EQ.LEN2)GO TO 300
    I1=I3
    I2=I1+NLEN
    IF(I2.LE.NUMPTS)GO TO 200
    L1=0
    GO TO 300

550 DO I=1,NUMPTS
       INDX(I)=TEMP(I)
    END DO
    NLEN=NLEN*2
    GO TO 100

800 CONTINUE

! return percentile values at 1% increments
  DO I=1,99
     DSORT(I)=VALUE(INDX(MAX(1,INT(I*NUMPTS/100.0))))
  END DO

DEALLOCATE (temp,indx)

END SUBROUTINE SORT
