subroutine statsco(stats_co,bwork,awork,ndata)
!$$$ subprogram documentation block
! . . . .
! subprogram: statsco prints statistics for carbon monoxide
! prgmmr: derber org: np23 date: 2003-05-22
!
! abstract: The routine computes and prints statistics regarding the
! use of carbon monoxide observations. Printed information includes
! that about data counts, quality control decisions,
! statistics based on the innovations, and penalties - all
! as a function of 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-12-23 treadon - use module jfunc to pass jiter,first
! 2005-11-29 derber - move ozmz array from obsmod to guess_grids
! 2006-02-03 derber - modify for new obs control and clean up output
! 2007-02-21 sienkiewicz - bring in stats for ozone level data
! 2007-05-30 h.liu - remove ozmz
! 2010-04023 tangborn - create statsco
!
! input argument list:
! stats_co - array holding sums from 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
! bwork - array containing information for statistics (level o3)
! awork - array containing information for data counts and gross checks (level o3)
!
! output argument list:
! stats_co - array holding sums from various statistical output
! bwork - array containing information for statistics (level o3)
! awork - array containing information for data counts and gross checks (level o3)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: zero,one
use obsmod, only: ndat,iout_co,dtype,dsis,dplat,ditype
use coinfo, only: error_co,nusis_co,nulev,iuse_co,jpch_co
use jfunc, only: jiter
use qcmod, only: npres_print,pboto3, ptopo3
use convinfo, only: nconvtype,ioctype
use gridmod, only: nsig
implicit none
! Declare passed variables
integer(i_kind),dimension(ndat,3) ,intent(in ) :: ndata
real(r_kind),dimension(9,jpch_co) ,intent(inout) :: stats_co
real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout) :: bwork
real(r_kind),dimension(7*nsig+100) ,intent(inout) :: awork
! Declare local variables
logical,dimension(ndat):: idisplay
integer(i_kind) j,iasim,i,icerr,ii
real(r_kind) svar,rsum,stdev,cpen,penalty_all,qcpenalty_all
real(r_kind),dimension(ndat):: rpenal,qcpenal
integer(i_kind),dimension(ndat):: icount_asim,iqccount_asim
real(r_kind) o3plty, o3qcplty, o3t, qco3t, rat, rat3
integer(i_kind) ntot, nread, nkeep, numgross, numfailqc, k, numlow, numhgh
integer(i_kind),dimension(nsig)::num
logical,dimension(nconvtype):: pflag
character(100) mesage
!******************************************************************************
! Compute and print statistics for carbon monoxide data
open(iout_co,position='append')
write(iout_co,850) jiter
850 format('OUTER ITERATION jiter=',i5)
! Print total penalty over all satellites
penalty_all=zero
qcpenalty_all=zero
do i=1,jpch_co
iasim = nint(stats_co(7,i))
if (iasim>0) then
penalty_all=penalty_all+stats_co(5,i)
qcpenalty_all=qcpenalty_all+stats_co(8,i)
endif
end do
write(iout_co,*)'carbon monoxide total penalty_all=',penalty_all
write(iout_co,*)'carbon monoxide total qcpenalty_all=',qcpenalty_all
! Print counts, bias, rms, stndev as a function of level
icount_asim=0
iqccount_asim=0
rpenal=zero
qcpenal=zero
idisplay=.false.
do ii=1,ndat
if(ditype(ii) == 'co')then
do i = 1,jpch_co
iasim = nint(stats_co(1,i))
if (iasim > 0 .and. nusis_co(i) == dsis(ii) ) then
if (iuse_co(i)==1) then
icount_asim(ii) = icount_asim(ii) + iasim
rpenal(ii) = rpenal(ii) + stats_co(5,i)
iqccount_asim(ii) = iqccount_asim(ii) + nint(stats_co(9,i))
qcpenal(ii) = qcpenal(ii) + stats_co(8,i)
end if
idisplay(ii)=.true.
end if
end do
end if
end do
do i = 1,jpch_co
iasim = nint(stats_co(1,i))
if (iasim > 0) then
svar = error_co(i)
if (iuse_co(i)/=1) svar = -svar
rsum = one/float(iasim)
icerr = nint(stats_co(2,i))
do j=3,6 ! j=3=obs-mod(w_biascor)
! j=4=(obs-mod(w_biascor))**2
! j=5=penalty contribution
! j=6=obs
stats_co(j,i) = stats_co(j,i)*rsum
end do
stats_co(4,i) = sqrt(stats_co(4,i))
if (iasim > 1) then
stdev = sqrt(stats_co(4,i)*stats_co(4,i)-stats_co(3,i)*stats_co(3,i))
else
stdev = zero
end if
write(iout_co,1102) i,nulev(i),nusis_co(i),iasim,icerr,svar,&
stats_co(6,i),stats_co(6,i)-stats_co(3,i),stats_co(3,i),&
stats_co(5,i),stats_co(4,i),stdev
endif
end do
! Write obs count to runtime output file
write(iout_co,1109)
do i=1,ndat
if (idisplay(i)) then
cpen=zero
if (icount_asim(i)>0) cpen=rpenal(i)/float(icount_asim(i))
write(iout_co,1115) jiter,dplat(i),dtype(i),ndata(i,2), &
ndata(i,3),icount_asim(i),rpenal(i),cpen,qcpenal(i),iqccount_asim(i)
endif
end do
2000 format(a7,2x,A4,6x,8(a7,1x))
2010 format(i7,1x,A10,1x,8(i7,1x))
2011 format(3x,f16.8,7(i7,1x))
2012 format(7x,A7,5x,7(a7,1x))
2999 format(' Illegal satellite type ')
1102 format(1x,i4,i4,1x,a20,2i7,1x,f8.3,1x,6(f11.7,1x))
1109 format(t5,'it',t11,'sat',t22,'inst',t36,'# read',t46,'# keep',t55,'# assim',&
t63,'penalty',t78,'cpen',t88,'qcpen',t101,'qcfail')
1115 format('o-g',1x,i2.2,1x,'co ',a10,1x,a10,1x,3(i9,1x),3(g12.5,1x),i8)
! End of carbon monoxide diagnostic print block.
! carbon monoxide level data diagnostics
o3plty=zero; o3qcplty=zero ; ntot=0
o3t=zero ; qco3t = zero;
nread = 0
nkeep = 0
do i=1,ndat
if (dtype(i)== 'o3lev') then
nread=nread+ndata(i,2)
nkeep=nkeep+ndata(i,3)
end if
end do
if (nkeep > 0) then
mesage='current fit of carbon monoxide level data, ranges in ppmv $'
do j = 1,nconvtype
pflag(j)=trim(ioctype(j)) == 'o3lev'
end do
call dtast(bwork,npres_print,pboto3,ptopo3,mesage,jiter,iout_co,pflag)
do k=1,nsig
num(k)=nint(awork(5*nsig+k+100))
rat=zero ; rat3=zero
if(num(k) > 0) then
rat=awork(6*nsig+k+100)/float(num(k))
rat3=awork(3*nsig+k+100)/float(num(k))
end if
ntot=ntot+num(k); o3plty=o3plty+awork(6*nsig+k+100)
o3qcplty=o3qcplty+awork(3*nsig+k+100)
write(iout_co,240) 'o3l',num(k),k,awork(6*nsig+k+100), &
awork(3*nsig+k+100),rat,rat3
end do
numgross=nint(awork(4))
numfailqc=nint(awork(21))
write(iout_co,925) 'o3l',numgross,numfailqc
numlow = nint(awork(2))
numhgh = nint(awork(3))
write(iout_co,900) 'o3l',numhgh,numlow
o3t=o3plty/ntot
qco3t=o3qcplty/ntot
end if
write(iout_co,950) 'o3l',jiter,nread,nkeep,ntot
write(iout_co,951) 'o3l',o3plty,o3qcplty,o3t,qco3t
close(iout_co)
240 format(' num(',A,') = ',i6,' at lev ',i4,' pen,qcpen,cpen,cqcpen = ',6(g12.5,1x))
925 format(' number of ',a5,' obs that failed gross test = ',I5,' nonlin qc test = ',I5)
900 format(' number of ',a5,' obs extrapolated above',&
' top sigma layer=',i8,/,10x,' number extrapolated below',&
' bottom sigma layer=',i8)
950 format(' type ',a5,' jiter ',i3,' nread ',i7,' nkeep ',i7,' num ',i7)
951 format(' type ',a5,' pen= ',e25.18,' qcpen= ',g13.6,' r= ',g13.6,' qcr= ',g13.6)
! End of routine
return
end subroutine statsco