!########################################################################
!########################################################################
!###### ######
!###### PROGRAM QPFSTATS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!########################################################################
!########################################################################
PROGRAM QPFSTATS,19
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads in two HDF files created by INTQPF (a forecast and a
! verification) and a binary bit map created by QPFMASK. The precip.
! fields in the HDF files extracted, bias and equitable threat score
! are calculated, and the statistics are output.
!
! AUTHOR: Eric Kemp, March 2000
!
! MODIFICATION HISTORY:
! Eric Kemp, 31 March 2000.
! Added lat/lon coordinates of four corners for NCL plotting.
!
!-----------------------------------------------------------------------
!
! Variable declarations
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
! Forecast variables
!
!-----------------------------------------------------------------------
INTEGER :: fdate(6)
INTEGER,ALLOCATABLE :: ftimesec(:)
INTEGER :: fmapproj
REAL :: fscale
REAL :: ftrulat1, ftrulat2, ftrulon
REAL :: fdx,fdy
CHARACTER*8 :: fmodel
CHARACTER*4 :: fgrid
CHARACTER*4,ALLOCATABLE,DIMENSION(:) :: fid_p,fid_sfc
CHARACTER*32,ALLOCATABLE,DIMENSION(:) :: fname_p,funit_p,fname_sfc, &
funit_sfc
REAL,ALLOCATABLE, DIMENSION(:,:,:,:,:) :: fvar_p
REAL,ALLOCATABLE, DIMENSION(:,:,:,:) :: fvar_sfc
REAL,ALLOCATABLE, DIMENSION(:) :: fpressure
REAL,ALLOCATABLE, DIMENSION(:,:) :: flat2d,flon2d
REAL :: fscswlat,fscswlon
REAL :: fmapswlat,fmapswlon,fmapnwlat,fmapnwlon, &
fmapselat,fmapselon,fmapnelat,fmapnelon
INTEGER :: fflag_p,fflag_sfc
REAL,ALLOCATABLE, DIMENSION(:,:,:) :: fprecip
INTEGER :: fnx,fny,fntime
INTEGER :: fnpreslevel,fnvar_p,fnvar_sfc
!-----------------------------------------------------------------------
!
! Verification variables
!
!-----------------------------------------------------------------------
INTEGER :: vdate(6)
INTEGER,ALLOCATABLE :: vtimesec(:)
INTEGER :: vmapproj
REAL :: vscale
REAL :: vtrulat1, vtrulat2, vtrulon
REAL :: vdx,vdy
CHARACTER*8 :: vmodel
CHARACTER*4 :: vgrid
CHARACTER*4,ALLOCATABLE,DIMENSION(:) :: vid_p,vid_sfc
CHARACTER*32,ALLOCATABLE,DIMENSION(:) :: vname_p,vunit_p,vname_sfc, &
vunit_sfc
REAL,ALLOCATABLE, DIMENSION(:,:,:,:,:) :: vvar_p
REAL,ALLOCATABLE, DIMENSION(:,:,:,:) :: vvar_sfc
REAL,ALLOCATABLE, DIMENSION(:) :: vpressure
REAL,ALLOCATABLE, DIMENSION(:,:) :: vlat2d,vlon2d
REAL :: vscswlat,vscswlon
REAL :: vmapswlat,vmapswlon,vmapnwlat,vmapnwlon, &
vmapselat,vmapselon,vmapnelat,vmapnelon
INTEGER :: vflag_p,vflag_sfc
REAL,ALLOCATABLE, DIMENSION(:,:,:) :: vprecip
INTEGER :: vnx,vny,vntime
INTEGER :: vnpreslevel,vnvar_p,vnvar_sfc
!-----------------------------------------------------------------------
!
! Bitmap variables
!
!-----------------------------------------------------------------------
INTEGER :: mdate(6)
INTEGER,ALLOCATABLE :: mtimesec(:)
INTEGER :: mmapproj
REAL :: mscale
REAL :: mtrulat1, mtrulat2, mtrulon
REAL :: mdx,mdy
CHARACTER*8 :: mmodel
CHARACTER*4 :: mgrid
CHARACTER*4,ALLOCATABLE,DIMENSION(:) :: mid_p,mid_sfc
CHARACTER*32,ALLOCATABLE,DIMENSION(:) :: mname_p,munit_p,mname_sfc, &
munit_sfc
REAL,ALLOCATABLE, DIMENSION(:,:,:,:,:) :: mvar_p
REAL,ALLOCATABLE, DIMENSION(:,:,:,:) :: mvar_sfc
REAL,ALLOCATABLE, DIMENSION(:) :: mpressure
REAL,ALLOCATABLE, DIMENSION(:,:) :: mlat2d,mlon2d
REAL :: mscswlat,mscswlon
REAL :: mmapswlat,mmapswlon,mmapnwlat,mmapnwlon, &
mmapselat,mmapselon,mmapnelat,mmapnelon
INTEGER :: mflag_p,mflag_sfc
INTEGER,ALLOCATABLE, DIMENSION(:,:) :: mask
INTEGER :: mnx,mny,mntime
INTEGER :: mnpreslevel,mnvar_p,mnvar_sfc
INTEGER :: found_mask
!-----------------------------------------------------------------------
!
! Other variables
!
! NOTE: CTA = array of "yes" forecast, "yes" observation grid points
! CTB = array of "yes" forecast, "no" observation grid points
! CTC = array of "no" forecast, "yes" observation grid points
! CTD = array of "no" forecast, "no" observation grid points
!
! Reference: Wilks, D. S., 1995: Statistical Methods in the
! Atmospheric Sciences. Academic Press, 476pp.
!
!-----------------------------------------------------------------------
REAL,ALLOCATABLE :: ets(:,:),bias(:,:)
INTEGER,ALLOCATABLE :: CTA(:,:),CTB(:,:),CTC(:,:),CTD(:,:)
INTEGER :: status,i,j,k,l
REAL :: factor
REAL,PARAMETER :: missing = -9999.
INTEGER,PARAMETER :: nstat = 2
CHARACTER*4,ALLOCATABLE,DIMENSION(:) :: ctstatsid
CHARACTER*32,ALLOCATABLE,DIMENSION(:) :: ctstatsname
REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ctstats
CHARACTER*32 :: threshunit
CHARACTER*72,PARAMETER :: line72 = &
'------------------------------------------------------------------------'
REAL :: vtrulat(2)
REAL,ALLOCATABLE :: vxsc(:),vysc(:)
REAL :: vx0,vy0
!-----------------------------------------------------------------------
!
! Namelists
!
!-----------------------------------------------------------------------
CHARACTER*256 :: finfilename
NAMELIST /hdf_finput/ finfilename
CHARACTER*256 :: vinfilename
NAMELIST /hdf_vinput/ vinfilename
CHARACTER*256 :: minfilename
NAMELIST /bitmap/ minfilename
INTEGER :: precipunit
INTEGER :: numthresh
INTEGER,PARAMETER :: maxnumthresh = 10
REAL :: thresh(maxnumthresh)
NAMELIST /threshold/ precipunit,numthresh,thresh
CHARACTER*256 :: hdfoutfilename,ascoutfilename
NAMELIST /output/ hdfoutfilename,ascoutfilename
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
fmodel=''
vmodel=''
WRITE(6,'(//5x,a)') &
&'###################################################################'
WRITE(6,'(5x,a,/5x,a)') &
&'# #',&
&'# Welcome to QPFSTATS, a program that reads in two precip. HDF #'
WRITE(6,'(5x,a)') &
&'# files from INTQPF and a binary bitmap file from QPFMASK, and #'
WRITE(6,'(5x,a)') &
&'# calculates equitable threat score and bias statistics. #', &
&'# #',&
&'###################################################################'
WRITE(6,*)
!-----------------------------------------------------------------------
!
! Read in namelist
!
!-----------------------------------------------------------------------
PRINT *, 'Reading NAMELIST hdf_finput'
READ(5,hdf_finput)
PRINT *, 'Reading NAMELIST hdf_vinput'
READ(5,hdf_vinput)
PRINT *, 'Reading NAMELIST bitmap'
READ(5,bitmap)
PRINT *, 'Reading NAMELIST threshold'
READ(5,threshold)
PRINT *, 'Reading NAMELIST output'
READ(5,output)
WRITE(6,*)'EMK: thresh = ',thresh
IF (numthresh > maxnumthresh) THEN
WRITE(6,*)'ERROR: numthresh > maxnumthresh'
WRITE(6,*)'Change maxnumthresh and recompile.'
WRITE(6,*)'numthresh = ',numthresh,' maxnumthresh = ', &
maxnumthresh
WRITE(6,*)'Aborting...'
STOP
END IF
IF (precipunit /= 1 .AND. precipunit /= 2) THEN
WRITE(6,*)'ERROR: Invalid precipunit value'
WRITE(6,*)'precipunit = ',precipunit
WRITE(6,*)'Aborting...'
STOP
ENDIF
!-----------------------------------------------------------------------
!
! Read forecast precipitation data.
!
!-----------------------------------------------------------------------
CALL rd_verif_dims(fnx,fny,fnpreslevel,fntime,fnvar_p,fnvar_sfc, &
finfilename)
ALLOCATE (ftimesec(fntime), fpressure(fnpreslevel), &
fvar_p(fnx,fny,fnpreslevel,fntime,fnvar_p), &
fvar_sfc(fnx,fny,fntime,fnvar_sfc), &
flat2d(fnx,fny),flon2d(fnx,fny), &
fprecip(fnx,fny,fntime), &
fid_p(fnvar_p),fname_p(fnvar_p),funit_p(fnvar_p), &
fid_sfc(fnvar_sfc),fname_sfc(fnvar_sfc), &
funit_sfc(fnvar_sfc), &
STAT=status)
IF (status /= 0) CALL alloc_fail
(status, 'f1')
CALL rdverif ( fnx,fny,fnpreslevel,fntime,fnvar_p,fnvar_sfc, &
finfilename,fmodel,fgrid,fdate,ftimesec,fpressure, &
fmapproj,fscale,ftrulat1,ftrulat2,ftrulon,fdx,fdy, &
fscswlat,fscswlon, &
fmapswlat,fmapswlon,fmapnwlat,fmapnwlon, &
fmapselat,fmapselon,fmapnelat,fmapnelon, &
fflag_p,fvar_p, fid_p, fname_p, funit_p, &
fflag_sfc,fvar_sfc, fid_sfc, fname_sfc, funit_sfc )
IF (fflag_sfc == 0) THEN
WRITE(6,*)'ERROR: Could not find surface data in HDF file.'
WRITE(6,*)'Aborting...'
STOP
ENDIF
!-----------------------------------------------------------------------
!
! Read verification precipitation data.
!
!-----------------------------------------------------------------------
CALL rd_verif_dims(vnx,vny,vnpreslevel,vntime,vnvar_p,vnvar_sfc, &
vinfilename)
ALLOCATE (vtimesec(vntime), vpressure(vnpreslevel), &
vvar_p(vnx,vny,vnpreslevel,vntime,vnvar_p), &
vvar_sfc(vnx,vny,vntime,vnvar_sfc), &
vlat2d(vnx,vny),vlon2d(vnx,vny), &
vprecip(vnx,vny,vntime), &
vid_p(vnvar_p),vname_p(vnvar_p),vunit_p(vnvar_p), &
vid_sfc(vnvar_sfc),vname_sfc(vnvar_sfc), &
vunit_sfc(vnvar_sfc), &
STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
CALL rdverif ( vnx,vny,vnpreslevel,vntime,vnvar_p,vnvar_sfc, &
vinfilename,vmodel,vgrid,vdate,vtimesec,vpressure, &
vmapproj,vscale,vtrulat1,vtrulat2,vtrulon,vdx,vdy, &
vscswlat,vscswlon, &
vmapswlat,vmapswlon,vmapnwlat,vmapnwlon, &
vmapselat,vmapselon,vmapnelat,vmapnelon, &
vflag_p,vvar_p, vid_p, vname_p, vunit_p, &
vflag_sfc,vvar_sfc, vid_sfc, vname_sfc, vunit_sfc )
IF (vflag_sfc == 0) THEN
WRITE(6,*)'ERROR: Could not find surface data in HDF file.'
WRITE(6,*)'Aborting...'
STOP
ENDIF
!-----------------------------------------------------------------------
!
! Check dimensions, map projection, etc. of forecast and verification.
!
!-----------------------------------------------------------------------
IF (fnx /= vnx .OR. fny /= vny .OR. fntime /= vntime .OR. &
fscale /= vscale .OR. &
ftrulat1 /= vtrulat1 .OR. ftrulat2 /= vtrulat2 .OR. &
ftrulon /= vtrulon .OR. fdx /= vdx .OR. fdy /= vdy .OR. &
fscswlat /= vscswlat .OR. fscswlon /= vscswlon) THEN
WRITE(6,*)'ERROR: Grid information for forecast and verification'
WRITE(6,*)'data do not match!'
WRITE(6,*)'fnx = ',fnx,' vnx = ',vnx
WRITE(6,*)'fny = ',fny,' vny = ',vny
WRITE(6,*)'fscale = ',fscale,' vscale = ',vscale
WRITE(6,*)'ftrulat1 = ',ftrulat1,' vtrulat1 = ',vtrulat1
WRITE(6,*)'ftrulat2 = ',ftrulat2,' vtrulat2 = ',vtrulat2
WRITE(6,*)'ftrulon = ',ftrulon,' vtrulon = ',vtrulon
WRITE(6,*)'fdx = ',fdx,' vdx = ',vdx
WRITE(6,*)'fdy = ',fdy,' vdy = ',vdy
WRITE(6,*)'fscswlat = ',fscswlat,' vscswlat = ',vscswlat
WRITE(6,*)'fscswlon = ',fscswlon,' vscswlon = ',vscswlon
WRITE(6,*)
WRITE(6,*)'Aborting...'
STOP
END IF
DO i=1,6
IF (fdate(i) /= vdate(i)) THEN
PRINT *, 'FATAL: Forecast and analysis start dates do not match'
PRINT *, 'fdate= ', fdate
PRINT *, 'vdate= ', vdate
STOP
END IF
END DO
!-----------------------------------------------------------------------
!
! Read bitmap file.
!
!-----------------------------------------------------------------------
CALL rd_verif_dims(mnx,mny,mnpreslevel,mntime,mnvar_p,mnvar_sfc, &
minfilename)
ALLOCATE (mtimesec(mntime), mpressure(mnpreslevel), &
mvar_p(mnx,mny,mnpreslevel,mntime,mnvar_p), &
mvar_sfc(mnx,mny,mntime,mnvar_sfc), &
mlat2d(mnx,mny),mlon2d(mnx,mny), &
mid_p(mnvar_p),mname_p(mnvar_p),munit_p(mnvar_p), &
mid_sfc(mnvar_sfc),mname_sfc(mnvar_sfc), &
munit_sfc(mnvar_sfc), &
STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
CALL rdverif ( mnx,mny,mnpreslevel,mntime,mnvar_p,mnvar_sfc, &
minfilename,mmodel,mgrid,mdate,mtimesec,mpressure, &
mmapproj,mscale,mtrulat1,mtrulat2,mtrulon,mdx,mdy, &
mscswlat,mscswlon, &
mmapswlat,mmapswlon,mmapnwlat,mmapnwlon, &
mmapselat,mmapselon,mmapnelat,mmapnelon, &
mflag_p,mvar_p, mid_p, mname_p, munit_p, &
mflag_sfc,mvar_sfc, mid_sfc, mname_sfc, munit_sfc )
IF (mflag_sfc == 0) THEN
WRITE(6,*)'ERROR: Could not find bitmap data in HDF file.'
WRITE(6,*)'Aborting...'
STOP
ENDIF
!-----------------------------------------------------------------------
!
! Check dimensions, map projection, etc. of bitmap and verification
!
!-----------------------------------------------------------------------
IF (mnx /= vnx .OR. mny /= vny .OR. mntime /= vntime .OR. &
mscale /= vscale .OR. &
mtrulat1 /= vtrulat1 .OR. mtrulat2 /= vtrulat2 .OR. &
mtrulon /= vtrulon .OR. mdx /= vdx .OR. mdy /= vdy .OR. &
mscswlat /= vscswlat .OR. mscswlon /= vscswlon) THEN
WRITE(6,*)'ERROR: Grid information for bitmap and verification'
WRITE(6,*)'data do not match!'
WRITE(6,*)'mnx = ',mnx,' vnx = ',vnx
WRITE(6,*)'mny = ',mny,' vny = ',vny
WRITE(6,*)'mscale = ',mscale,' vscale = ',vscale
WRITE(6,*)'mtrulat1 = ',mtrulat1,' vtrulat1 = ',vtrulat1
WRITE(6,*)'mtrulat2 = ',mtrulat2,' vtrulat2 = ',vtrulat2
WRITE(6,*)'mtrulon = ',mtrulon,' vtrulon = ',vtrulon
WRITE(6,*)'mdx = ',mdx,' vdx = ',vdx
WRITE(6,*)'mdy = ',mdy,' vdy = ',vdy
WRITE(6,*)'mscswlat = ',mscswlat,' vscswlat = ',vscswlat
WRITE(6,*)'mscswlon = ',mscswlon,' vscswlon = ',vscswlon
WRITE(6,*)
WRITE(6,*)'Aborting...'
STOP
END IF
!-----------------------------------------------------------------------
!
! Dump HDF data into precip arrays, changing units as necessary.
!
!-----------------------------------------------------------------------
DO i = 1,fnvar_sfc
IF (fid_sfc(i) == 'APCP') THEN
IF (funit_sfc(i) == 'mm' .AND. precipunit == 1) THEN
fprecip = fvar_sfc(:,:,:,i)
ELSE IF (funit_sfc(i) == 'mm' .AND. precipunit == 2) THEN
fprecip = fvar_sfc(:,:,:,i)
factor = 39.37/1000.
CALL cvtpcpunit(fnx,fny,fntime,fprecip,factor,missing)
ELSE IF (funit_sfc(i) == 'in' .AND. precipunit == 2) THEN
fprecip = fvar_sfc(:,:,:,i)
ELSE IF (funit_sfc(i) == 'in' .AND. precipunit == 1) THEN
fprecip = fvar_sfc(:,:,:,i)
factor = 1000./39.37
CALL cvtpcpunit(fnx,fny,fntime,fprecip,factor,missing)
ELSE
WRITE(6,*)'ERROR: funit_sfc = ',funit_sfc
STOP
ENDIF
ENDIF
END DO
DO i = 1,vnvar_sfc
IF (vid_sfc(i) == 'APCP') THEN
IF (vunit_sfc(i) == 'mm' .AND. precipunit == 1) THEN
vprecip = vvar_sfc(:,:,:,i)
ELSE IF (vunit_sfc(i) == 'mm' .AND. precipunit == 2) THEN
vprecip = vvar_sfc(:,:,:,i)
factor = 39.37/1000.
CALL cvtpcpunit(vnx,vny,vntime,vprecip,factor,missing)
ELSE IF (vunit_sfc(i) == 'in' .AND. precipunit == 2) THEN
vprecip = vvar_sfc(:,:,:,i)
ELSE IF (vunit_sfc(i) == 'in' .AND. precipunit == 1) THEN
vprecip = vvar_sfc(:,:,:,i)
factor = 1000./39.37
CALL cvtpcpunit(vnx,vny,vntime,vprecip,factor,missing)
ELSE
WRITE(6,*)'ERROR: vunit_sfc = ',vunit_sfc
STOP
ENDIF
ENDIF
END DO
!-----------------------------------------------------------------------
!
! Dump HDF bitmap data into mask array.
!
!-----------------------------------------------------------------------
ALLOCATE (mask(mnx,mny),STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
mask = missing
found_mask = 0
DO i = 1,mnvar_sfc
IF (mid_sfc(i) == 'MASK') THEN
mask = INT(mvar_sfc(:,:,1,i))
found_mask = 1
END IF
END DO
IF (found_mask == 0) THEN
WRITE(6,*)'ERROR: Could not find bitmap. Aborting...'
STOP
END IF
!-----------------------------------------------------------------------
!
! Calculate equitable threat score and bias
!
!-----------------------------------------------------------------------
ALLOCATE (bias(vntime,numthresh),ets(vntime,numthresh), &
CTA(vntime,numthresh),CTB(vntime,numthresh), &
CTC(vntime,numthresh),CTD(vntime,numthresh), &
STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
CALL calcetsb(vnx,vny,vntime,vprecip,fprecip,numthresh,thresh, &
mask,bias,ets,CTA,CTB,CTC,CTD)
100 FORMAT(1x,a,f4.2,a)
110 FORMAT(1x,A,i5,A,i5,A,i5)
!-----------------------------------------------------------------------
!
! Dump statistical data to a new HDF file.
!
!-----------------------------------------------------------------------
ALLOCATE(ctstats(vntime,numthresh,nstat),ctstatsid(nstat), &
ctstatsname(nstat),STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
ctstatsid(1) = 'BIAS'
ctstatsid(2) = 'ETS '
ctstatsname(1) = 'Bias'
ctstatsname(2) = 'Equitable Threat Score'
IF (precipunit == 1) THEN
threshunit = 'mm'
ELSE IF (precipunit == 2) THEN
threshunit = 'in'
ELSE
WRITE(6,*)'ERROR: Invalid precipunit value'
WRITE(6,*)'precipunit = ',precipunit
WRITE(6,*)'Aborting...'
STOP
END IF
DO i = 1,nstat
IF (i == 1) THEN
ctstats(:,:,i) = bias
ELSE IF (i == 2) THEN
ctstats(:,:,i) = ets
ENDIF
END DO
CALL wrtqpfstats(vnx,vny,vntime,numthresh,nstat,missing, &
hdfoutfilename,vmodel,vgrid, &
vdate,vtimesec,vmapproj,vscale,vtrulat1,vtrulat2, &
vtrulon,vdx,vdy,vscswlat,vscswlon,&
vmapswlat,vmapswlon,vmapnwlat,vmapnwlon, &
vmapselat,vmapselon,vmapnelat,vmapnelon, &
thresh,threshunit, &
ctstatsid,ctstatsname,ctstats, &
CTA,CTB,CTC,CTD)
!-----------------------------------------------------------------------
!
! Output ASCII data.
!
!-----------------------------------------------------------------------
ALLOCATE(vxsc(vnx),vysc(vny),STAT=status)
IF (status /= 0) CALL alloc_fail (status, 'f1')
vtrulat(1) = vtrulat1
vtrulat(2) = vtrulat2
CALL setmapr(vmapproj,vscale,vtrulat,vtrulon)
CALL lltoxy(1,1,vscswlat,vscswlon,vx0,vy0)
DO i=1,vnx
vxsc(i)=vx0+(i-1)*vdx
END DO
DO j=1,vny
vysc(j)=vy0+(j-1)*vdy
END DO
CALL xytoll(vnx,vny,vxsc,vysc,vlat2d,vlon2d)
WRITE(6,*)'Writing text output to ',TRIM(ascoutfilename)
OPEN(UNIT=13,FILE=ascoutfilename,STATUS="REPLACE")
WRITE(13,600) 'Forecast initialized ',fdate(1),'-',fdate(2), &
'-',fdate(3),'.',fdate(4),':',fdate(5),':',fdate(6)
WRITE(13,610) 'Forecast precipitation data:', TRIM(fmodel)
WRITE(13,610) 'Analysis precipitation data:', TRIM(vmodel)
WRITE(13,*)
WRITE(13,*) 'Verification Region:'
WRITE(13,625) NINT((vnx)*vdx/1000.0), NINT((vny)*vdx/1000.0), &
vdx/1000.0, vnx,vny
WRITE(13,650) 'NW/NE', &
vlat2d(1,vny),vlon2d(1,vny),vlat2d(vnx,vny),vlon2d(vnx,vny)
WRITE(13,650) 'SW/SE', &
vlat2d(1,1),vlon2d(1,1), vlat2d(vnx,1),vlon2d(vnx,1)
WRITE(13,*)'Used bit map from file ',TRIM(minfilename)
WRITE(13,*)
DO l = 1, numthresh
WRITE(13,*) line72
WRITE(13,810) 'Precipitation Threshold: ',thresh(l), TRIM(threshunit)
WRITE(13,*) 'Time Grid False '
WRITE(13,*) '(hr) Points Hits Misses Alarms Bias', &
' ETS'
DO k = 1,vntime
WRITE(13,700) vtimesec(k)/3600, &
CTA(k,l)+CTB(k,l)+CTC(k,l)+CTD(k,l), &
CTA(k,l),CTC(k,l),CTB(k,l), &
bias(k,l),ETS(k,l)
END DO
END DO
600 FORMAT (/1x,a21,i4.4,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2)
610 FORMAT (1X,A,1X,A)
625 FORMAT (' Domain:',I5, ' x', I5, ' km, dx=', F5.1, &
' km (',I4,' x',I4,' )' )
650 FORMAT (1X,A,' corner lat/lon:', 2F10.3, 2X, 2F10.3)
700 FORMAT (3x,i2,4x,i4,2x,i4,4x,i4,4x,i4,2x,f8.2,3x,f8.2)
810 FORMAT ( 1X,a,f4.2," (",A,")" )
!-----------------------------------------------------------------------
!
! The end.
!
!-----------------------------------------------------------------------
WRITE(6,*)'Program QPFSTATS successfully completed.'
END PROGRAM QPFSTATS
!########################################################################
!########################################################################
!###### ######
!###### SUBROUTINE CALCETSB ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!########################################################################
!########################################################################
SUBROUTINE calcetsb(nx,ny,ntime,vprecip,fprecip,numthresh,thresh, & 1
mask,bias,ets,cta,ctb,ctc,ctd)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculates equitable threat score and bias scores for precipitation.
!
! AUTHOR: Eric Kemp, March 2000
!
!-----------------------------------------------------------------------
!
! Variable declarations
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx,ny,ntime,numthresh
REAL, INTENT(IN) :: vprecip(nx,ny,ntime), fprecip(nx,ny,ntime)
REAL, INTENT(OUT) :: ETS(ntime,numthresh),Bias(ntime,numthresh)
REAL, INTENT(IN) :: thresh(numthresh)
INTEGER, INTENT(IN) :: mask(nx,ny)
INTEGER, INTENT(INOUT), DIMENSION(ntime,numthresh) :: CTA,CTB,CTC,CTD
!-----------------------------------------------------------------------
!
! Variables for statistical calculations
!
! LA = Number of "yes" forecast, "yes" observation grid points.
! LB = Number of "yes" forecast, "no" observation grid points.
! LC = Number of "no" forecast, "yes" observation grid points.
! LD = Number of "no" forecast, "no" observation grid points.
! R = Number of random correct forecasts expected due to chance within
! the total number of model and observation pairs.
! ETS = Equitable Threat Score.
!
! NOTE: R = (LA + LB)*(LA + LC)/(LA + LB + LC + LD)
! ETS = (LA - R)/(LA + LB + LC - R)
! Bias = (LA + LB)/(LA + LC)
!
!-----------------------------------------------------------------------
REAL :: R
INTEGER :: LA, LB, LC, LD
!-----------------------------------------------------------------------
!
! Other variables
!
!-----------------------------------------------------------------------
INTEGER :: i,j,k,l
REAL, PARAMETER :: missing = -9999.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Bias = missing
ETS = missing
CTA = 0
CTB = 0
CTC = 0
CTD = 0
DO l = 1,numthresh
DO k = 1,ntime
LA = 0
LB = 0
LC = 0
LD = 0
DO j = 1,ny
DO i = 1,nx
IF (mask(i,j) == 1 .AND. fprecip(i,j,k) /= missing .AND. &
vprecip(i,j,k) /= missing ) THEN
IF (fprecip(i,j,k) >= thresh(l) .AND. &
vprecip(i,j,k) >= thresh(l)) THEN
LA = LA + 1
ELSE IF (fprecip(i,j,k) >= thresh(l) .AND. &
vprecip(i,j,k) < thresh(l)) THEN
LB = LB + 1
ELSE IF (fprecip(i,j,k) < thresh(l) .AND. &
vprecip(i,j,k) >= thresh(l)) THEN
LC = LC + 1
ELSE
LD = LD + 1
ENDIF
ENDIF
END DO
END DO
IF ((LA + LB + LC + LD) == 0) THEN
R = missing
ELSE
R = REAL(LA + LB)*REAL(LA + LC)/REAL(LA + LB + LC + LD)
ENDIF
IF ((R == missing) .OR. (LA + LB + LC - R) == 0 ) THEN
ETS(k,l) = missing
ELSE
ETS(k,l) = REAL(LA - R)/REAL(LA + LB + LC - R)
ENDIF
IF ((LA+LC) == 0) THEN
Bias(k,l) = missing
ELSE
Bias(k,l) = REAL(LA+LB)/REAL(LA+LC)
ENDIF
CTA(k,l) = LA
CTB(k,l) = LB
CTC(k,l) = LC
CTD(k,l) = LD
END DO
END DO
END SUBROUTINE calcetsb
!########################################################################
!########################################################################
!###### ######
!###### SUBROUTINE CVTPCPUNIT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!########################################################################
!########################################################################
SUBROUTINE cvtpcpunit(nx,ny,ntime,precip,factor,missing) 4
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Converts units of precipitation array.
!
! AUTHOR: Eric Kemp, March 2000
!
!-----------------------------------------------------------------------
!
! Variable declarations
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER,INTENT(IN) :: nx,ny,ntime
REAL,INTENT(INOUT) :: precip(nx,ny,ntime)
REAL,INTENT(IN) :: factor,missing
!-----------------------------------------------------------------------
!
! Other variables
!
!-----------------------------------------------------------------------
INTEGER :: i,j,k
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
DO k = 1,ntime
DO j = 1,ny
DO i = 1,nx
IF (precip(i,j,k).ne.missing) THEN
precip(i,j,k) = precip(i,j,k)*factor
END IF
END DO
END DO
END DO
END SUBROUTINE cvtpcpunit