peakfqSA-jfe.f [docs/17c/peakfqSA-jfe-distrib/src] Revision: Date:
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****| PeakfqSA
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c this routine requests the name of a specification file which
c contains input information for implementing the EMA flood
c frequency algorithm. The specification file is loosely based
c on the input file for USGS program PEAKFQ
c
c this was prepared by tim cohn, us geological survey, to support
c development of an operational version of peakfq (bulletin 17b
c implementation) which would include the expected moments
c algorithm (ema; cohn et al. 1997; 2001)
c
C Minor modifications and additions
C by John England, USACE john.f.england@usace.army.mil
C
C modifications by John England denoted with C JFE
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c development history
c
c timothy a. cohn 15 aug 2007 (v. 0.910)
c modified 21 sep 2007 (v. 0.920)
c modified 01 feb 2008 (v. 0.933)
c modified 28 feb 2008 (v. 0.934)
c modified 26 jun 2008 (v. 0.935)
c modified 27 jun 2008 (v. 0.936)
c modified 03 jul 2008 (v. 0.937)
c major mods 08 mar 2011 (v. 0.950)
c modified 15 jun 2011 (v. 0.960)
c modified 05 oct 2011 (v. 0.970)
c modified 19 oct 2011 (v. 0.971)
c changes to emafit 19 oct 2011 (v. 0.972)
c added var_est 08 feb 2012 (v. 0.973)
c added var_est 09 feb 2012 (v. 0.974)
c Alt. MGB (JFE) 09 feb 2012 (v. 0.980)
c Redone MGB (TAC) 26 sep 2012 (v. 0.981)
c New Var Comp (TAC) 15 Oct 2012 (v. 0.990)
c Tweak to MGB (TAC) 15 Oct 2012 (v. 0.991)
c as_G_ERL (TAC) 16 Feb 2013 (v. 0.992)
c Full Inv Quad. (TAC) 05 Mar 2013 (v. 0.993)
c Minor correct. (TAC) 31 Mar 2013 (v. 0.994)
c Minor correct. (TAC) 03 Apr 2013 (v. 0.995)
c Minor correct. (TAC) 28 May 2013 (v. 0.996) ; plotting pos.;MSE[G]
c Add MaxFracLO(TAC) 20 Aug 2013 (v. 0.997)
c Minor cor prt (TAC) 06 Feb 2014 (v. 0.997)
c emafit change (TAC) 06 Feb 2014 (v. 0.998) ; nu_min=5.d0
C
C expand quantile output (JFE) 26 Feb 2016 (v. 0.998)
C
C includes modifications to routines in emafit-jfe.f
C
C expand arrays for paleoflood data (JFE) 07 Sep 2017 (v. 0.999)
C
C includes modifications to: H-S plotting routine in probfun.f,
C array lengths in peakfqSA, emafit-jfe.f, imslfake.f
C
C 28-JUN-2018 - JFE
C Testing reveals errors in reporting of EMA Thresholds, data, plotting positions
C diagnose and fix array lengths
C
C Minor corrections paleoflood data (JFE) 28 Jun 2018 (v. 0.9992)
C Turn off MGBT 3-sweep reporting and max frac
C
C 07-MAY-2019 - JFE
C Output format fixes - paleofloods
C
C
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c compiler options:
c gfortran -o PeakfqSQ0996.out -fbounds-check -static -Wampersand -g
c peakfqSA.f
c emafit.f
c probfun.f
c imslfake.f
c utilityfuncs.f
c dcdflib1.f90
c
c for OSX (good to try -Wall, too)
c gfortran -o PeakfqSA0996.app -fbounds-check -Wampersand -g peakfqSA.f emafit.f probfun.f utilityfuncs.f imslfake.f *.f90
c for Windows
c gfortran -o PeakfqSA0996.exe -fbounds-check -Wampersand -g -static peakfqSA.f emafit.f probfun.f utilityfuncs.f imslfake.f *.f90
c
c note that "-static" is needed to create dll-free code; otherwise users
c will get run-time errors suggesting the dll is minQ.
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c Notes:
c 1) Observations are classified as "Systematic" iff the threshold
c corresponding to the observation does not exceed the gage-base
c discharge at the site.
c 2) Low-outlier identification is done using one of three options:
c "NONE"; "GBT" (Grubbs-Beck Test); "MGBT" (Multiple Grubbs-BeckTest)
c The option is stored in gbtype
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
program PeakfqSA
implicit none
real version
! parameter(version=0.9992)
parameter(version=0.9993)
integer p_nq,p_npp,p_nt
C JFE expand quantiles; npp corresponds to nqt in emafit-jfe.f
C JFE see the expanded quantiles in /blkdata within emafit-jfe.f
c parameter (p_nq=20000,p_npp=32,p_nt=100)
parameter (p_nq=20000,p_npp=41,p_nt=100)
double precision minQ,lmissing,lthresh_inf
parameter (minQ=1.d-99,lmissing=-80.d0,lthresh_inf=99.d0)
double precision qP3,mseg
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
integer
1 i,ind(p_nq),indx,j,n_val,n_true
integer
1 n,n_peaks,n_not_used,n_recs,nt,nlow,nzero,nlow_V,nuncens,
2 nint,ncensl,ncensr,npp
integer
1 tyr_beg(p_nq),tyr_end(p_nq)
integer
1 beg_year,end_year,hist_period,n_historic,n_sys
integer
1 wy(p_nq),wy_in(p_nq)
integer
1 nCLT,Nspc,Nspcmax
double precision
1 mean_gen,mean_mse,mean_sd
double precision
1 sd_gen,sd_mse,sd_sd
double precision
1 skew_mse,skew_gen,skew_sd
double precision
1 ql_in(p_nq),qu_in(p_nq),lq_sys(p_nq),lqCLT,pctCLT,w
double precision
1 pvaluew
character*4 q_type(p_nq),gbtype,signif,filext
double precision
1 lql(p_nq),lqu(p_nq)
double precision
1 ltl(p_nq),ltu(p_nq),tl_in(p_nt),tu_in(p_nt),qs
double precision
1 gage_base_q,min_gbq,hi_sys,hi_thresh,latitude,longitude,lo_hist,
2 urb_reg
double precision
1 pq(p_npp),yp(p_npp),ci_low(p_npp),ci_high(p_npp),var_est(p_npp),
2 gbthresh0,lgbthresh0,gbcrit,gbthresh,gbcrit_V,gbthresh_V
double precision
1 cmomsM(3,3),cmoms(3),cmomsl(3),cmomsB17B(3)
equivalence
1 (cmoms(1),cmomsM(1,1)),
2 (cmomsl(1),cmomsM(1,2)),
2 (cmomsB17B(1),cmomsM(1,3))
character*80
1 data_filename,outfile,specfile,station,output_filename,
2 CSV_outfile,cmdfile
double precision
1 mu,sig,ppall(p_nq),ppallx(p_nq),xsynth(p_nq)
double precision
1 PPalpha,PPthr(p_nq),PPpeT(p_nq)
integer
1 PPnt,PPnb(p_nq),ns
double precision fp_z_icdf
double precision
1 as_M_mse,as_S2_mse,as_G_mse,as_G_mse_Syst,as_G_ERL
double precision Alphaout,Alphain,Alphazeroin,MaxFracLO
character*132
1 line,alpha(20),alpha1
character*132
1 mean_opt
character*132
1 sd_opt
character*132
1 skew_opt
logical lprint(p_nq),lprnt2(p_nq),ex,CSV_file,CLT,Lcmdfile
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c N.B. lql2,lqu2,ltl2,ltu2 are the logs of adjusted values coming from gbtest
c They are not the original data, but rather the modified data
c based on the Grubbs-Beck test
c
integer nx
C JFE
C parameter (nx=25000)
C JFE 28Jun2018 revert back to 25000, fixes errors in
C parameter (nx=20000)
parameter (nx=25000)
double precision lql2,lqu2,ltl2,ltu2
integer nGBiter
character*4 at_site_option,at_site_default,at_site_std
C JFE common /tacg01/gbcrit,gbthresh,pvaluew(10000),qs(10000),
common /tacg01/gbcrit,gbthresh,pvaluew(20000),qs(20000),
1 ns,nlow,nzero,gbtype
common /tacg02/lql2(nx),lqu2(nx),ltl2(nx),ltu2(nx)
common /tacg03/gbcrit_V(10),gbthresh_V(10),nlow_V(10),nGBiter
common /tacg04/at_site_option,at_site_default,at_site_std
common /tac005/as_M_mse,as_S2_mse,as_G_mse,as_G_mse_Syst,as_G_ERL
common /MGB001/Alphaout,Alphain,Alphazeroin
common /MGB002/MaxFracLO
character*4 VarS2opt
common /tacR01/VarS2opt
double precision eps ! nominal confidence interval coverage
common /tacci1/eps
double precision qdum1(nx)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c beginning of executable code
c
write(*,'(///,80(''*''),//)')
write(*,'(t40,a12)') 'EMA ANALYSIS'
C JFE write(*,'(t34,a20,f6.3)') 'Computed by PeakfqSA',version
write(*,'(t34,a20,f8.4)') 'Computed by PeakfqSA',version
C JFE write(*,'(t39,a14)') 'Tim Cohn, USGS'
C JFE write(*,'(t39,a15)') 'tacohn@usgs.gov'
C JFE write(*,'(t39,a14)') ' 703/648-5711 '
C JFE write(*,'(///,80(''*''),//)')
C
C JFE
C
C write(*,'(///,80(''*''),//)')
C JFE write(*,'(t34,a30)') 'with PeakfqSA modifications by'
write(*,'(t39,a19)') 'John England, USACE'
write(*,'(t39,a29)') 'john.f.england@usace.army.mil'
write(*,'(t39,a14)') ' 303/963-4524 '
write(*,'(///,80(''*''),//)')
Lcmdfile = .false.
Nspcmax = 1
call getarg(1,specfile)
if(filext(specfile) .eq. '.cmd') then
cmdfile = specfile
inquire(file=cmdfile,exist=ex)
if(ex) then
open(unit=10,file=cmdfile,status='old')
Lcmdfile = .true.
Nspcmax = 999
else
write(*,*) 'File does not exist: ',cmdfile
stop
endif
endif
do 999 Nspc=1,Nspcmax
if(Lcmdfile) read(10,'(a80)',end=9999) specfile
write(*,'(1x,''Input File: '',a80)') specfile
if(specfile .ne. '') goto 8
5 write(*,*) 'Enter SPECIFICATION filename (a80)'
read(*,*) specfile
8 continue
inquire(file=specfile,exist=ex)
if(ex) then
open(unit=11,file=specfile,status='old')
else
write(*,*) 'File does not exist: ',specfile
goto 5
endif
call makefn_out(specfile,'.out',12,outfile)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c set default parameter values
c
station = ' '
mean_opt = 'STATION'
mean_gen = 1.d0
mean_mse = 1.d99
mean_sd = 1.d99
C
sd_opt = 'STATION'
sd_gen = 1.d0
sd_mse = 1.d99
sd_sd = 1.d99
C
skew_opt = 'STATION'
skew_gen = 0.d0
skew_mse = 1.d99
C
C JFE
C beg_year = -10000
C end_year = +10000
beg_year = -20000
end_year = +20000
hist_period = 0
lo_hist = minQ
gbthresh0 = minQ
gbtype = 'MGBT'
hi_sys = minQ
hi_thresh = minQ
gage_base_q = minQ
latitude = lmissing
longitude = lmissing
CSV_file = .FALSE.
CLT = .FALSE.
min_gbq = minQ
at_site_std = at_site_default
PPalpha = 0.d0 ! Weibull plotting positions
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c read in specifications
c
nt = 0
do 10 i=1,1000
read(11,'(a80)',end=15) line
call readio1(line,n_val,alpha)
call cvrtupcase(alpha(1),alpha1)
if(alpha1 .eq. 'I') then
read(alpha(2),*) data_filename
else if(alpha1 .eq. 'O') then
read(alpha(2),*) output_filename
else if(alpha1 .eq. 'STATION') then
j = (80-len(trim(line(9:88))))/2
station(1+j:80-j) = line(9:88)
c regional mean options
else if(alpha1 .eq. 'MEANOPT') then
call cvrtupcase(alpha(2),mean_opt)
else if(alpha1 .eq. 'GENMEAN') then
read(alpha(2),*) mean_gen
else if(alpha1 .eq. 'MEANMSE') then
read(alpha(2),*) mean_mse
mean_sd = sqrt(mean_mse)
else if(alpha1 .eq. 'MEANSD') then
read(alpha(2),*) mean_sd !
mean_mse = mean_sd**2
c regional standard deviation options
else if(alpha1 .eq. 'SDOPT') then
call cvrtupcase(alpha(2),sd_opt)
else if(alpha1 .eq. 'SDGEN') then
read(alpha(2),*) sd_gen
else if(alpha1 .eq. 'SDMSE') then
read(alpha(2),*) sd_mse
sd_sd = sqrt(sd_mse)
else if(alpha1 .eq. 'SDSD') then
read(alpha(2),*) sd_sd !
sd_mse = sd_sd**2
c regional skew options
else if(alpha1 .eq. 'SKEWOPT') then
call cvrtupcase(alpha(2),skew_opt)
else if(alpha1 .eq. 'GENSKEW') then
read(alpha(2),*) skew_gen
else if(alpha1 .eq. 'A_S_SKEW_OPT') then
at_site_std = alpha(2)(1:4)
else if(alpha1 .eq. 'SKEWMSE') then
read(alpha(2),*) skew_mse
skew_sd = sqrt(skew_mse)
else if(alpha1 .eq. 'SKEWSD') then
read(alpha(2),*) skew_sd !
skew_mse = skew_sd**2
else if(alpha1 .eq. 'BEGYEAR') then
read(alpha(2),*) beg_year
else if(alpha1 .eq. 'ENDYEAR') then
read(alpha(2),*) end_year
else if(alpha1 .eq. 'HISTPERIOD') then ! total length of historical period
read(alpha(2),*) hist_period
else if(alpha1 .eq. 'URB/REG') then
read(alpha(2),*) urb_reg
else if(alpha1 .eq. 'LOHIST') then ! smallest historic peak
read(alpha(2),*) lo_hist
else if(alpha1 .eq. 'LOMETHOD') then ! low-outlier threshold
call cvrtupcase(alpha(2),gbtype) ! 'NONE', 'GBT ', 'MBGT', 'FIXE'
C WRITE(*,'(''GB1'',2A40,A4)') ALPHA(1),ALPHA(2),GBTYPE ! CTAC
else if(alpha1 .eq. 'LOTHRESH') then ! low-outlier threshold
read(alpha(2),*) gbthresh0
else if(alpha1 .eq. 'Alphaout') then ! inward sweep alpha level
read(alpha(2),*) Alphaout
else if(alpha1 .eq. 'Alphain') then ! outward sweep alpha level
read(alpha(2),*) Alphain
else if(alpha1 .eq. 'Alphazeroin') then ! alpha level sweep from X[1]
read(alpha(2),*) Alphazeroin
else if(alpha1 .eq. 'MAXFRACLO') then ! maximum fraction low outliers
read(alpha(2),*) MaxFracLO
else if(alpha1 .eq. 'HISYS') then ! highest systematic peak
read(alpha(2),*) hi_sys
else if(alpha1 .eq. 'HITHRESH') then
read(alpha(2),*) hi_thresh
else if(alpha1 .eq. 'GAGEBASE') then
read(alpha(2),*) gage_base_q
gage_base_q = max(min_gbq,gage_base_q)
else if(alpha1 .eq. 'LATITUDE') then
read(alpha(2),*) latitude
else if(alpha1 .eq. 'LONGITUDE') then
read(alpha(2),*) longitude
else if(alpha1 .eq. 'PP_ALPHA') then ! This determines Weibull/Hazen/etc.
read(alpha(2),*) PPalpha ! plotting positions
else if(alpha1 .eq. 'EPS') then
read(alpha(2),*) eps
else if(alpha1 .eq. 'CSV') then
if(alpha(2) .eq. 'YES') then
CSV_file = .TRUE.
call makefn_out(specfile,'.csv',14,CSV_outfile)
endif
else if(alpha1 .eq. 'CLT') then ! Censor the left-hand tail
read(alpha(2),*) pctCLT
if(pctCLT .gt. 0.d0) then
CLT = .TRUE.
pctCLT = max(50.d0,pctCLT) ! limit the level of deliberate censoring
endif
else if( (alpha1 .eq. 'PCPT_THRESH') .or.
1 (alpha1.eq. 'THRESHOLD') ) then
nt = nt+1
read(alpha(2),*) tyr_beg(nt)
read(alpha(3),*) tyr_end(nt)
read(alpha(4),*) tl_in(nt)
read(alpha(5),*) tu_in(nt)
endif
10 continue
15 continue
write(*,*) 'Specification lines read: ',i-1
close(11)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c set up low-outlier test (note: < -6 implies standard test (see notes))
c
lgbthresh0 = log10(max(minQ,gbthresh0))
n = 1 + (end_year - beg_year)
hist_period = max(n,hist_period) ! default value of hist_period is 0
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c set up regional mean info
c
if(mean_opt .eq. 'WEIGHTED') then ! weighted mean
mean_mse = mean_mse
else if(mean_opt .eq. 'STATION') then ! at-site mean w/uncertainty
mean_mse = 1.d10
else if(mean_opt .eq. 'GENERALIZED') then ! reg. mean w/o uncertainty
mean_mse = 1.d-20
c else if(mean_opt .eq. 'GENERAL/VAR') then ! reg. mean w/uncertainty
c mean_mse = -mean_mse
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c set up regional standard deviation info
c
if(sd_opt .eq. 'WEIGHTED') then ! weighted sd (default)
sd_mse = sd_mse
else if(sd_opt .eq. 'STATION') then ! at-site sd w/uncertainty
sd_mse = 1.d10
else if(sd_opt .eq. 'GENERALIZED') then ! reg. sd w/o uncertainty
sd_mse = 1.d-20
c else if(sd_opt .eq. 'GENERAL/VAR') then ! reg. sd w/uncertainty
c sd_mse = -sd_mse
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c set up regional skew info
c
if(skew_opt .eq. 'WEIGHTED') then ! weighted skew (default)
skew_mse = skew_mse
else if(skew_opt .eq. 'STATION') then ! at-site skew w/uncertainty
skew_mse = 1.d10
else if(skew_opt .eq. 'GENERALIZED') then ! reg. skew w/o uncertainty
skew_mse = 1.d-20
else if(skew_opt .eq. 'GENERAL/VAR') then ! reg. skew w/uncertainty
skew_mse = -skew_mse
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c read in flow data (note format)
c N.B. There are only two types of flow data: Q and QINT
c
inquire(file=specfile,exist=ex)
if(ex) then
open(unit=13,file=data_filename,status='old')
else
write(*,*) 'Data file does not exist: ',data_filename
stop
endif
n_recs = 0
do 20 i=1,p_nq
read(13,'(a80)',end=25) line
call readio1(line,n_val,alpha)
call cvrtupcase(alpha(1),alpha1)
q_type(i) = 'Othr'
if(alpha1 .eq. 'QINT') then
n_recs = n_recs+1
read(alpha(2),*) wy_in(n_recs)
read(alpha(3),*) ql_in(n_recs)
read(alpha(4),*) qu_in(n_recs)
else if(alpha1 .eq. 'Q') then
n_recs = n_recs+1
read(alpha(2),*) wy_in(n_recs)
read(alpha(3),*) ql_in(n_recs)
qu_in(n_recs) = ql_in(n_recs)
endif
20 continue
25 continue
write(*,*) 'Data lines read: ',i-1
close(13)
do 30 i = 1,n
wy(i) = (i + beg_year-1)
lql(i) = lmissing
lqu(i) = lmissing ! lthresh_inf*0.99 Make obs. range huge XXXX????
ltl(i) = lthresh_inf ! Assume no information from this year
ltu(i) = lthresh_inf ! unless user provides valid threshold
q_type(i) = "Cens" ! Assume we don't know anything about q-type
do 32 j=1,nt
if(tyr_beg(j) .le. wy(i) .and. tyr_end(j) .ge. wy(i)) then
if(tl_in(j) .le. gage_base_q) q_type(i) = "Syst" ! Note Def'n
ltl(i) = log10(max(minQ,tl_in(j)))
ltu(i) = log10(max(minQ,tu_in(j)))
endif
32 continue
if(nt .eq. 0) then ! Problem: No threshold
write(*,*) ' ** Error in Specification File: **'
write(*,*) ' ** Observation thresholds minQ **'
stop
endif
30 continue
n_peaks = 0
n_not_used = 0
n_sys = 0
n_historic = 0
do 40 i = 1,n_recs
if(wy_in(i) .ge. beg_year .and. wy_in(i) .le. end_year) then
indx = 1 + wy_in(i) - beg_year
lql(indx) = log10(max(minQ,ql_in(i)))
lqu(indx) = log10(max(minQ,qu_in(i)))
n_peaks = n_peaks + 1
if(ltl(indx) .le. log10(gage_base_q)) then
n_sys = n_sys + 1
lq_sys(n_sys)= lqu(indx)
ltl(indx) = log10(gage_base_q) ! Cannot measure q < gage-base
else
n_historic = n_historic + 1
q_type(indx) = 'Hist'
endif
else
n_not_used = n_not_used + 1
endif
40 continue
if(CLT) then
call dsvrgn(n_sys,lq_sys,lq_sys)
nCLT = (n_sys+1)*pctCLT/100.d0
if(nCLT .lt. 1 .or. nCLT .ge. (n_sys-3)) then
lqCLT = -99.d0
else
w = (n_sys+1.d0)*pctCLT/100.d0 - nCLT
lqCLT = (1.d0-w)*lq_sys(nCLT) + w*lq_sys(nCLT+1)
endif
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c For years without peaks, we encounter two cases
c 1) if CLT (censored left tail) fix tl
c 2) peak is less than threshold tl, if threshold reasonable
c 3) peak is completely unknown, if tl is not reasonable
c
do 50 i=1,n
if(q_type(i) .eq. 'Syst' .and. lqu(i) .le. ltl(i)) then
write(*,1001) wy(i),10**lql(i),10**lqu(i),
1 10**ltl(i),10**ltu(i)
write(12,1001) wy(i),10**lql(i),10**lqu(i),
1 10**ltl(i),10**ltu(i)
1001 format(78('*'),//,' *** Check for error in input data ***',/,
2 ' Zero flow or missing data in systematic record',/,
3 '==> (WY,Ql,Qu,Tl,Tu)',(i8,2f15.2,2d15.5))
if(lqu(i) .eq. lmissing) then
write(*,*) ' There is is likely an error in the'
write(*,*) ' input data value for WY, Q or Threshold'
write(*,*) ' PeakfqSA (line 485) '
endif
endif
if(CLT) then
ltl(i) = max(ltl(i),lqCLT)
endif
if(lqu(i) .lt. ltl(i)) then ! censored peak for this year
lqu(i) = ltl(i) ! recode qu to lower threshold
endif
50 continue
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c compress data (eliminate years for which we have no information)
c
n_true = 0
nint = 0
ncensl = 0
ncensr = 0
nuncens = 0
do 60 i=1,n
if(ltl(i) .lt. lthresh_inf) then
n_true = n_true + 1
wy(n_true) = wy(i)
lql(n_true) = lql(i)
lqu(n_true) = lqu(i)
ltl(n_true) = ltl(i)
ltu(n_true) = ltu(i)
q_type(n_true) = q_type(i)
if(lql(n_true) .ne. lqu(n_true)) then
nint = nint + 1
else
nuncens = nuncens+1
endif
if(lql(n_true) .lt. ltl(n_true)) ncensl = ncensl + 1
if(lqu(n_true) .gt. ltu(n_true)) ncensr = ncensr + 1
endif
60 continue
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c perform calculations
c
call emafitpr(n_true,
1 lql,lqu,ltl,ltu,q_type,
2 mean_gen,mean_mse,
3 sd_gen,sd_mse,
4 skew_gen,skew_mse,
5 gbtype,lgbthresh0,
6 cmomsM,pq,npp,yp,ci_low,ci_high,var_est)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c begin output of results
c
write(12,'(///,80(''*''),//)')
write(12,'(t40,a12)') 'EMA ANALYSIS'
C JFE write(12,'(t33,a20,f6.3)') 'Computed by PeakfqSA',version
write(12,'(t33,a20,f8.4)') 'Computed by PeakfqSA',version
C JFE write(12,'(t39,a14)') 'Tim Cohn, USGS'
C JFE write(12,'(t39,a15)') 'tacohn@usgs.gov'
C JFE write(12,'(t39,a14)') ' 703/648-5711 '
c
C
C JFE
C
C write(12,'(//,80(''*''),/)')
C write(12,'(t34,a30)') 'with PeakfqSA modifications by'
write(12,'(t39,a19)') 'John England, USACE'
write(12,'(t39,a29)') 'john.f.england@usace.army.mil'
write(12,'(t39,a14)') ' 303/963-4524 '
write(12,'(//,80(''*''),/)')
C JFE
write(12,'(t34,a18)') 'Run Date and Time'
call GET_DATE_TIME(12)
write(12,'(/,80(''*''),//)')
write(12,'(5x,a80)') station
write(12,*)
write(12,'(t30,'' Spec. file'',t50,a80)') specfile
write(12,'(t30,'' Data file '',t50,a80)') data_filename
write(12,*)
write(12,'(t30,'' Historical Period Begins'',t59,i9)') wy(1)
write(12,'(t30,'' Historical Period Ends'',t59,i9)') wy(n)
write(12,'(t30,'' Length of Historical Period'',t59,i9)') n
write(12,'(t30,'' Total No. of Observations'',t59,i9)') n_true
write(12,'(t30,'' Years Without Information'',t59,i9)') n-n_true
write(12,'(t30,'' Number of Peaks'',t59,i9)') n_peaks
write(12,'(t30,'' Peaks Not Used'',t59,i9)') n_not_used
write(12,'(t30,'' Point or Finite Int. Peaks'',t59,i9)') n_sys
write(12,'(t30,'' Historic Peaks'',t59,i9)') n_historic
write(12,'(t30,'' Years of Historic Record'',t59,i9)')
1 hist_period
write(12,'(t30,'' Number of Uncensored Peaks'',t59,i9)') nuncens
write(12,'(t30,'' Number of Interval Peaks'',t59,i9)') nint
write(12,'(t30,'' Number Left Censored Peaks'',t59,i9)') ncensl
write(12,'(t30,'' Number Right Censored '',t59,i9)') ncensr
write(12,'(t30,'' Number Zero Flows '',t59,i9)') nzero
c regional mean is employed
write(12,'(t30,'' Mean Option'',t61,a11)') mean_opt
if(mean_opt .ne. 'STATION') then
write(12,'(t30,'' Generalized Mean (mean_gen)'',t60,g12.4)')
1 mean_gen
write(12,'(t30,'' MSE[M_gen]'',t60,g12.4)') mean_mse
endif
c regional sd is employed
write(12,'(t30,'' Std. Dev. Option'',t61,a11)') sd_opt
if(sd_opt .ne. 'STATION') then
write(12,'(t30,'' -S2 Weighting Option'',t57,a11)') VarS2opt
write(12,'(t30,'' Generalized Std Dev (Sd_gen)'',t60,g12.4)')
1 sd_gen
write(12,'(t30,'' Standard Error SD[Sd_gen]'',t60,g12.4)')
1 sd_sd
write(12,'(t30,'' MSE[Sd_gen]'',t60,g12.4)') sd_mse
write(12,'(t30,'' df[Sd_gen]'',t60,g12.4)')
1 sd_gen**2/(2.d0*sd_mse)
endif
c
write(12,'(t30,'' At-Site Skew Method '',t59,5x,a4)') at_site_std
write(12,'(t30,'' Skew Option'',t61,a11)') skew_opt
write(12,'(t30,'' Generalized Skew'',t59,f9.2)') skew_gen
write(12,'(t30,'' Standard Error'',t59,f9.4)') skew_sd
write(12,'(t30,'' (Pseudo)-MSE Employed'',t60,g12.4)') skew_mse
c
write(12,'(t30,'' Gage Base Discharge'',t59,f9.1)') gage_base_q
if(CLT) then
write(12,'(t30,'' Cens. Left Tail Pct'',t59,f9.1)') pctCLT
write(12,'(t30,'' Cens. Threshold (Q)'',t59,f9.1)') 10**lqCLT
endif
if(CSV_file) then
write(14,'(a80)') station
write(14,'(''Spec. file'','','',a80)') specfile
write(14,'(''Data file '','','',a80)') data_filename
write(14,'(''Historical Period Begins'','','',i9)') wy(1)
write(14,'(''Historical Period Ends'','','',i9)') wy(n)
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c low outlier information
c
write(12,'(t30,'' Low-Outlier Type'',t64,a4)') gbtype
if(gbtype .eq. 'GBT ') then
if(nGBiter .gt. 1) then
write(12,'(t30,'' Low-Outlier Bound'',t57,''Iterated GB'')')
else
write(12,'(t30,'' Low-Outlier Bound'',t57,
1 ''Grubbs-Beck Test'')')
endif
nlow = 0
do 90 j=1,nGBiter
write(12,'(t30,'' Iteration'',t59,i9)') j
write(12,'(t30,'' - No. LOs in Iter.'',i2,t59,i9)')
1 j,nlow_V(j)
write(12,'(t30,'' - Crit. Val. Iter.'',i2,t59,f9.2)')
1 j,10**gbcrit_V(j)
write(12,'(t30,'' - Thresh. Val. Iter.'',i2,t59,f9.2)')
1 j,10**gbthresh_V(j)
nlow = nlow + nlow_V(j)
90 continue
else if (gbtype .eq. 'FIXE') then
write(12,'(t30,'' Low-Outlier Bound'',t55,''User Supplied'')')
write(12,'(t30,'' User supplied LO Bound'',t59,f9.1)')
1 10**gbthresh
nlow = nlow_V(1)
else if (gbtype .eq. 'MGBT') then
C JFE
C JFE write(12,'(t30,'' Three-Sweep Low-Outlier Detection'')')
write(12,'(t30,'' Two-Sweep Low-Outlier Detection'')')
C JFE write(12,'(t32,"Sweep Out: Alpha Out[N2]:",t59,f9.6)')
write(12,'(t32,"Sweep Out: Alpha Out:",t59,f9.6)')
1 Alphaout
C JFE write(12,'(t32,"Sweep 2: Alpha In[J1]:",t59,f9.6)')
C JFE 1 Alphain
C JFE write(12,'(t32,"Sweep 3: Alpha In[1]:",t59,f9.6)')
write(12,'(t32,"Sweep In: Alpha In:",t59,f9.6)')
1 Alphazeroin
C JFE write(12,'(t32,"Max Frac LO:",t59,f9.6)')
C JFE 1 MaxFracLO
write(12,'(t30,'' Low-Outlier Bound'',t56,''Multiple G-B'')')
write(12,'(t30,'' MGB L-O Threshold'',t59,f9.1)')
1 10**gbthresh
nlow = nlow_V(1)
write(12,'(t30,'' Number of Low Outliers'',t59,i9)') nlow
do 91 j=1,nlow+1
write(12,'(t32,'' P-value obs:'',i3,'' (Q='',
1 e11.4,'') p='',f6.4,a4)')
2 j,qs(j),pvaluew(j),signif(pvaluew(j))
91 continue
endif
c write(12,'(t30,'' Total Number Interval Obs.'',t59,i9)') nint
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c print out annual frequency curve quantiles and ci_s
c
write(12,2002)
2002 format(//,t20,'ANNUAL FREQUENCY CURVE PARAMETERS ',
1 ' -- LOG-PEARSON TYPE III',//)
write(12,'(t10,''Fitted Log10 Moments'',t50,''M'',
1 t60,''S'',t70,''G'')')
write(12,'(t10,65(''-''))')
write(12,'(t10,''EMA, At-Site Data, w/o Reg. Info'',t45,3f10.6)')
1 cmomsl(1),sqrt(cmomsl(2)),cmomsl(3)
write(12,'(t10,'' EMA w/ Reg. Info & B17B MSE(G)'',t45,3f10.6)')
1 cmomsB17B(1),sqrt(cmomsB17B(2)),cmomsB17B(3)
write(12,'(t10,''EMA w/Reg. Info & Spec. MSE(G)'',t45,3f10.6)')
1 cmoms(1),sqrt(cmoms(2)),cmoms(3)
c MSEs and weighting schemes
if( (mean_sd .ge. 0.d0) .and. (mean_sd .lt. 1.d9) ) then
write(12,'(//,t15,''Estimates & MSEs Corr. to Mean'')')
write(12,'(t17,''Option'',t47,a12)') mean_opt
write(12,'(t17,''EMA Est. of M (M_atsite)'',t45,f12.6)')
1 cmomsl(1)
write(12,'(t17,''MSE[M_atsite]'',t45,f12.6)')
1 as_M_mse
write(12,'(t17,''Regional M (M_reg)'',t45,f12.6)') mean_gen
write(12,'(t17,''MSE[M_reg]'',t45,f12.6)') mean_mse
write(12,'(t17,''Weighted M'',t45,f12.6)') cmoms(1)
endif
if( (sd_sd .ge. 0.d0) .and. (sd_sd .lt. 1.d9) ) then
write(12,'(/,t15,''Estimates & MSEs Corr. to S (std dev)'')')
write(12,'(t17,''Option'',t47,a12)') sd_opt
write(12,'(t17,''EMA Est. of S (S_atsite)'',t45,f12.6)')
1 sqrt(cmomsl(2))
write(12,'(t17,''MSE(S^2_atsite)'',t45,f12.6)')
1 as_S2_mse
write(12,'(t17,''MSE(S_atsite)'',t45,f12.6)')
1 as_S2_mse/(4.d0*cmomsl(2))
write(12,'(t17,''df(S_atsite)'',t45,f12.6)')
1 2.d0*(cmomsl(2))**2/as_S2_mse
write(12,'(t17,''Regional S (S_reg)'',t45,f12.6)') sd_gen
write(12,'(t17,''MSE[S_reg]'',t45,f12.6)') sd_mse
write(12,'(t17,''df(S_reg)'',t45,f12.6)')
1 sd_gen**2/(2.d0*sd_mse)
write(12,'(t17,''Weighted S'',t45,f12.6)') sqrt(cmoms(2))
endif
C JFE adjusted width of output in this table from 45 to 49
write(12,'(/,t15,''Estimates & MSEs Corr. to Skew '')')
C JFE write(12,'(t17,''Option'',t47,a12)') skew_opt
write(12,'(t17,''Option'',t53,a12)') skew_opt
C JFE write(12,'(t17,''EMA Est. of G_atsite'',t45,f12.6)')
write(12,'(t17,''EMA Est. of G_atsite'',t49,f12.6)')
1 cmomsl(3)
C JFE write(12,'(t17,''B17B Est. of MSE[G_atsite]'',t45,f12.6)')
write(12,'(t17,''B17B Est. of MSE[G_atsite]'',t49,f12.6)')
1 mseg(n_true,cmomsl(3))
write(12,
1 '(t17,''EMA ('',A4,'') Est. of MSE[G_atsite]'',t49,f12.6)')
1 at_site_option,as_G_mse
write(12,'(t17,''Regional G (G_reg)'',t49,f12.6)') skew_gen
C JFE write(12,'(t17,''MSE[G_reg]'',t45,f12.6)') skew_mse
write(12,'(t17,''MSE[G_reg]'',t49,g12.6)') skew_mse
write(12,'(t17,''Weighted G'',t49,f12.6)') cmoms(3)
write(12,'(t17,''MSE[G_atsite_systematic]'',t49,f12.6)')
1 as_G_mse_Syst
write(12,'(t17,''ERL[G_atsite]'',t49,f12.6)')
1 as_G_ERL
c Print the frequency curve
write(12,2003) 100.*eps
C JFE
C 2003 format(//,t35,'EMA FREQUENCY ESTIMATES',//,77('-'),/,
2003 format(//,t35,'EMA FREQUENCY ESTIMATES',//,80('-'),/,
1 t56,f5.1,'% CI on EMA Est.',/,
1 'Ann. Exc. Prob.',
2 t18,'EMA Est.',
3 t29,'V[log(EMA)]',
4 t42,'At-Site Est.',
5 t57,'CI-Low',
6 t69,'CI-High',/,77('-'))
if(CSV_file) then
C JFE write(14,'(a12,3('','',g10.5),//)') 'Moments',
write(14,'(a12,3('','',f11.6),//)') 'Moments',
1 cmoms(1),sqrt(cmoms(2)),cmoms(3)
write(14,'(a12,5('','',a12))') 'PP','RetPeriod','Z-Score',
1 'Q_Est','CI_L','CI_U'
endif
do 410 i=1,npp
write(12,2004) 1.d0-pq(i),
1 10**yp(i),
2 var_est(i),
3 10**qp3(pq(i),cmomsl),
4 10**ci_low(i),
5 10**ci_high(i)
C JFE
2004 format('EMA-Q',1x,f9.7,f13.2,f13.5,3f13.2)
C JFE 2004 format('EMA-Q',1x,f6.4,f13.2,f13.5,3f13.2)
if(CSV_file) then
C JFE write(14,'(g12.5,5('',''g12.5))')
C JFE write(14,'(g20.7,5('',''g20.7))')
write(14,2014)
1 1.d0-pq(i),1.d0/(1.d0-pq(i)),fp_z_icdf(pq(i)),
2 10**yp(i),10**ci_low(i),10**ci_high(i)
2014 format(f9.7,',',f15.3,',',f10.7,',',f13.2,',',f13.2,','f13.2)
endif
410 continue
write(12,'(///,80(''*''),//,''1'')')
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c generate plotting position data
c
CALL PLOTPOSHS(N_TRUE,LQL2,LQU2,LTL2,LTU2,PPALPHA,
1 QDUM1,PPALLX,ppnt,PPthr,pppet,PPnb)
DO 987 I=1,N_TRUE
PPALL(I) = 1.D0-PPALLX(I)
987 CONTINUE
call porder(lql2,n_true,ind)
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c print out ema representation of annual peak floods and corresponding
c thresholds
c
write(12,'(//,t33,''PERCEPTION THRESHOLDS FOR EMA'',/,
1 t33,29(''-''),/)')
write(12,1002)
1002 format(t11,61('-'),/,t11,' T',t16,'Beg. Year',t26,'End Year',
1 t48,'T_lower',t63,'T_upper',/,t11,61('-'))
do 70 j=1,nt
Cjfe write(12,'(t11,i2,t21,i4,t30,i4,t45,1p,g12.5,3x,g12.5)')
write(12,'(t11,i2,t21,i6,t30,i6,t45,1p,g12.5,3x,g12.5)')
1 j,tyr_beg(j),tyr_end(j),tl_in(j),min(tu_in(j),1.d50)
70 continue
if(CLT) then
write(12,'(t11,a3,t21,i4,t30,i4,t45,1p,g12.5,3x,g12.5)')
1 'CLT',beg_year,end_year,10**lqCLT
endif
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c print out ema thresholds, corresponding plotting positions,
c and number of observations censored beneath each threshold
c --This info is used to correctly locate glyphs on graph (TAC 21Jun11)
c
write(12,'(///,t28,''THRESHOLD EXCEED. PROBS & CENSORED OBS'',/,
1 t28,38(''-''),/)')
write(12,1012)
1012 format(t11,61('-'),/,t11,'Sorted(T)',t23,'T_lower',t36,
1 'Non-Exc. Prob.',
1 t54,'Cens. Obs. below T',/,t11,61('-'))
do 80 j=1,PPnt
write(12,'(t11,i2,t21,1p,g12.4,t40,0p,f8.5,t60,i6)')
1 j,10**PPthr(j),1.d0-PPpeT(j),PPnb(j)
80 continue
write(12,'(////,t35,''EMA REPRESENTATION OF DATA'',
1 /,t35,26(''-''))')
write(12,'(/,5x,a80,/)') station
do 300 i=1,n_true
if(lql2(i) .eq. lqu2(i)) then
lprint(i) = .TRUE.
else if(i .eq. 1 .or. i .eq. n) then
lprint(i) = .TRUE.
else if(lql2(i) .ne. lql2(i-1)) then
lprint(i) = .TRUE.
else if(lqu2(i) .ne. lqu2(i-1)) then
lprint(i) = .TRUE.
else if(ltl2(i) .ne. ltl2(i-1)) then
lprint(i) = .TRUE.
else if(ltu2(i) .ne. ltu2(i-1)) then
lprint(i) = .TRUE.
else if(wy(i) .ne. wy(i-1)+1) then
lprint(i) = .TRUE.
else
lprint(i) = .FALSE.
endif
C JFE Need to check/modify this
if(lql2(i) .gt. lqu2(i) - 5.d0) then ! if it is a fixed interval
lprnt2(i) = .TRUE. ! print plotting position
else
lprnt2(i) = .FALSE. ! o.w. dont
endif
300 continue
c
if(CSV_file) then
write(14,'(//,a12,5('','',a12))') 'WY','PP','Z-Score','Ql',
1 'Qu','Q-fit'
endif
write(12,1004)
1004 format(80('-'),/,t3,'Year',t13,'Q_obs.',t24,'Q_lower',
1 t36,'Q_upper',t47,'T_lower',t63,'T_upper',t75,'Q_type',
2 /,80('-'))
do 310 i=1,n_true
c
do 311 j=1,n_recs ! find the original datapoint for this wateryear
if(wy_in(j) .eq. wy(i)) goto 312
311 continue
312 continue
c
if(lprint(i) .or. lprint(i+1)) then
C JFE print lqu2 in g13.3 if infinity
if(10**lqu2(i).lt.1000000000.d0) then
write(12,'(1x,i5,t9,f10.1,2(1x,f12.1),2(1x,1p,g13.3),0p,a8)')
1 wy(i),ql_in(j),10**lql2(i),10**lqu2(i),10**ltl2(i),
2 min(10**ltu2(i),1.d10),q_type(i)
else
write(12,'(1x,i5,t9,f10.1,1x,f12.1,1x,g13.3,
3 2(1x,1p,g13.3),0p,a8)')
1 wy(i),ql_in(j),10**lql2(i),10**lqu2(i),10**ltl2(i),
2 min(10**ltu2(i),1.d10),q_type(i)
endif
C JFE
else if(lprint(i-1)) then
write(12,'(t35,'' . . .'')')
endif
310 continue
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c display above-threshold obserations and respective plotting positions
c
write(12,'(///,t20,'' SORTED OBSERVATIONS '')')
write(12,' (t20,''(~Hirsch/Stedinger Plotting Positions)'')')
write(12,' (t22,''Alpha for Plotting Postions:'',f5.2)')
1 PPAlpha
write(12,'(T20,37(''-''),//)')
write(12,1005)
1005 format(t13,53('-'),/,t13,'Year',t18,'Plot Pos',t31,'Obs. Q',
1 t43,'Fit Value Q',t60,'% Diff',
2 /,t13,53('-'))
do 320 i=n_true,1,-1
if(lprnt2(ind(i))) then
write(12,'(1x,i15,f9.4,2(g14.5,1x),f10.2)')
1 wy(ind(i)),1.d0-ppall(ind(i)),10**lql2(ind(i)),
2 10**qP3(ppall(ind(i)),cmoms),
3 100*(1.d0-(10**lql2(ind(i))/10**qP3(ppall(ind(i)),cmoms)))
if(CSV_file)
1 write(14,'(i12,5('','',g12.5))')
2 wy(ind(i)),1.d0-ppall(ind(i)),fp_z_icdf(ppall(ind(i))),
3 10**lql2(ind(i)),10**lqu2(ind(i)),
4 10**qP3(ppall(ind(i)),cmoms)
endif
320 continue
write(*,*) 'Finished Processing: ',specfile
write(*,*) 'Normal Termination'
write(*,*)
close(12)
999 continue
9999 continue
write(*,*)
write(*,'(''Processed '',i3,'' specification files'',/)') Nspc-1
stop
end