subroutine statspcp(aivals,ndata)
!$$$ subprogram documentation block
! . . . .
! subprogram: statspcp prints statistics for precipitation obs
! prgmmr: derber org: np23 date: 2003-05-22
!
! abstract: The routine computes and prints statistics regarding the
! use of precipitation observations. Printed information
! includes that about data counts, quality control decisions,
! statistics based on the innovations, and penalties - all
! as a function of precipitation observation and satellite
! type
!
! program history log:
! 2003-05-22 derber
! 2004-06-15 treadon - update documentation
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-10-06 parrish - increase diagnostic array sizes and add
! output from nonlinear qc
! 2004-12-23 treadon - use module jfunc to pass jiter,first
! 2005-03-09 parrish - add output of count of precip obs rejected
! by nonlinear qc
! 2005-09-28 derber - replace dtype parameter with obstype_all
! 2006-02-03 derber - modify for new obs control and to clean up output
! 2006-07-28 derber - use r1000 from constants
!
! input argument list:
! aivals - array holding sums for various statistical output
! ndata(*,1)- number of observation keep for processing
! ndata(*,2)- number of observation read
! ndata(*,3)- number of observations keep after read
!
! output argument list:
! aivals - array holding sums for various statistical output
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
! Include data type module
use kinds, only: r_kind,i_kind
use pcpinfo, only: dtphys,deltim
use obsmod, only: dtype,ndat,iout_pcp,dplat,ditype
use jfunc, only: jiter,first
use constants, only: zero,half,one,two,r1000,r3600
implicit none
! Declare passed variables
integer(i_kind),dimension(ndat,3) ,intent(in ) :: ndata
real(r_kind) ,dimension(40,ndat),intent(inout) :: aivals
! Declare local variables
integer(i_kind) i,j,is
integer(i_kind) iread,ikeep,iasim
integer(i_kind) nsphys
integer(i_kind) iland
integer(i_kind) isum,icoast,icerr
integer(i_kind) isnoice,itotal,ireduce
integer(i_kind) kdt,iratio,idzmax
integer(i_kind) igradsml,igradbig,idetect
integer(i_kind) inumw,inum,iclip,itoss
integer(i_kind) ntossqc
real(r_kind) sdv,rterm
real(r_kind) pcpobs0,pcpnbc0,pcpbc0,pcpobs1,pcpnbc1,pcpbc1
real(r_kind) pcpsas0,pcpsas1,rsum
real(r_kind) frain,dtp,dtf
real(r_kind) fhour,rtime
real(r_kind) penalty_all,qcpenalty_all
real(r_kind) rmmhr
real(r_kind),dimension(ndat):: rpen,cpen,qcpen,qccpen
logical,dimension(ndat):: display
! Initialize variables
! Open unit to receive printout
open(iout_pcp,position='append')
! On first outer iteration, write constant parameters
if (first ) then
nsphys = max(int(two*deltim/dtphys+0.9999_r_kind),1)
dtp = two*deltim/nsphys
dtf = half*dtp
frain = dtf/dtp
kdt = 1
fhour = kdt*deltim/r3600
rtime = one/(r3600*fhour)
rmmhr = r1000*rtime * r3600
write(iout_pcp,100)'deltim,dtphys=',deltim,dtphys
write(iout_pcp,100)'dtp,dtf =',dtp,dtf
write(iout_pcp,100)'frain,fhour =',frain,fhour
write(iout_pcp,100)'rtime,rmmhr =',rtime,rmmhr
write(iout_pcp,110)'nsphys,kdt =',nsphys,kdt
100 format(a14,1x,2(g25.18,1x))
110 format(a14,1x,2(i6,1x))
endif
! Loop over all observational data types
display=.false.
! Print information regarding qc.
write(iout_pcp,*)' '
write(iout_pcp,2000) 'obtype','num','numw',&
'itotal','ireduce','iland','isnoice','icoast',&
'grd<e-10','grd>e4','ratio>3','smooth',&
'thrhold'
2000 format(a10,1x,15(a8,1x))
do is=1,ndat
inum = nint(aivals(2,is))
display(is) = ditype(is) == 'pcp'
! If precipitation observation has nonzero number of obs, generate output
if (inum > 0 .and. display(is)) then
inumw = nint(aivals(1,is))
itotal = nint(aivals(3,is))
ireduce = nint(aivals(4,is))
iland = nint(aivals(5,is))
isnoice = nint(aivals(6,is))
icoast = nint(aivals(7,is))
iclip = nint(aivals(12,is))
igradsml = nint(aivals(26,is))
igradbig = nint(aivals(27,is))
iratio = nint(aivals(28,is))
idzmax = nint(aivals(29,is))
idetect = nint(aivals(30,is))
write(iout_pcp,2010) dtype(is),inum,inumw,&
itotal,ireduce,iland,isnoice,icoast,&
igradsml,igradbig,iratio,idzmax,idetect
2010 format(a10,1x,12(i8,1x))
end if
end do
write(iout_pcp,2001) 'obtype','qc-flag','satpcp','pcpnbc',&
'pcpsas','pcplrg','pcpbc'
2001 format(a7,1x,a7,1x,5(a18,1x),a6)
penalty_all=zero
qcpenalty_all=zero
ntossqc=0
do is=1,ndat
inum = nint(aivals(2,is))
! If precipitation observation has nonzero number of obs, generate output
if (inum > 0 .and. display(is)) then
rterm=zero
if (aivals(25,is) > zero) rterm=one/aivals(25,is)
pcpobs0 = aivals(21,is)*rterm
pcpnbc0 = aivals(22,is)*rterm
pcpbc0 = aivals(23,is)*rterm
pcpsas0 = aivals(24,is)*rterm
ikeep = nint(aivals(25,is))
if (aivals(35,is) > zero) rterm=one/aivals(35,is)
pcpobs1 = aivals(31,is)*rterm
pcpnbc1 = aivals(32,is)*rterm
pcpbc1 = aivals(33,is)*rterm
pcpsas1 = aivals(34,is)*rterm
itoss = nint(aivals(35,is))
write(iout_pcp,2013)dtype(is),'keep',ikeep,pcpobs0,pcpnbc0,pcpsas0,&
pcpnbc0-pcpsas0,pcpbc0
write(iout_pcp,2013)dtype(is),'toss',itoss,pcpobs1,pcpnbc1,pcpsas1,&
pcpnbc1-pcpsas1,pcpbc1
2013 format(a10,1x,a4,1x,i6,1x,5(g19.12,1x))
! Accumualte total penalty
penalty_all=penalty_all+aivals(15,is)
qcpenalty_all=qcpenalty_all+aivals(39,is)
ntossqc=ntossqc+nint(aivals(40,is))
end if
end do
do is=1,ndat
inum = nint(aivals(2,is))
! If precipitation observation has nonzero number of obs, generate output
if (inum > 0 .and. display(is)) then
write(iout_pcp,2012) dtype(is),aivals(15,is)
2012 format(A10,' penalty = ',g19.12)
end if
end do
! Print total precipitation penalty
write(iout_pcp,*)' '
write(iout_pcp,3000)penalty_all
3000 format('pcp total penalty_all =',g25.18)
write(iout_pcp,3001)qcpenalty_all
3001 format('pcp total qcpenalty_all =',g25.18)
write(iout_pcp,3002)ntossqc
3002 format('pcp total failed nonlinqc=',i8)
! Print counts, bias, rms, stndev as a function of observation type
do is = 1,ndat
isum = nint(aivals(11,is))
if (isum > 0 .and. display(is)) then
rpen(is) = aivals(15,is)
qcpen(is) = aivals(39,is)
rsum = one/float(isum)
icerr = nint(aivals(12,is))
do j=13,16
aivals(j,is)=aivals(j,is)*rsum
end do
aivals(14,is) = sqrt(aivals(14,is))
cpen(is) = aivals(15,is)
qccpen(is) = aivals(39,is)
if (isum > 1) then
sdv = sqrt(aivals(14,is)*aivals(14,is)- &
aivals(13,is)*aivals(13,is))
else
sdv = zero
endif
write(iout_pcp,1102) dtype(is),dplat(is),isum,icerr,&
aivals(16,is),aivals(13,is),aivals(15,is),aivals(14,is),sdv
1102 format(a10,1x,a10,1x,2(i6,1x),6(g17.10,1x))
else
rpen(is) = zero
qcpen(is) = zero
cpen(is) = zero
qccpen(is) = zero
end if
end do
write(iout_pcp,1109)
1109 format(t5,'it',t13,'sat',t28,'# read',t39,'# keep',t49,'# assim',&
t58,'penalty',t74,'cpen',t82,'qcpnlty',t95,'qccpen')
do i=1,ndat
if (aivals(11,i)>zero .and. display(i)) then
iread=ndata(i,2)
ikeep=nint(aivals(1,i))
iasim=nint(aivals(11,i))
write(iout_pcp,1115) jiter,dtype(i),iread,ikeep,iasim,&
rpen(i),cpen(i),qcpen(i),qccpen(i)
endif
end do
1115 format('o-g',1x,i2.2,1x,'pcp',2x,a10,2x,3(i9,2x),4(g12.5,1x))
! Close output unit
close(iout_pcp)
! End of routine
return
end subroutine statspcp