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',& 'grde4','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