!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRDTOSNG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE grdtosng(nx,ny,nz,nz_tab,mxsng,nvar,nobsng, & 1,2
xs,ys,zp,icatg,anx,qback,hqback,nlvqback, &
shght,su,sv,spres,stheta,sqv, &
stnsng,isrcsng,icatsng,elevsng,xsng,ysng, &
obsng,qobsng,qualsng,odifsng,oanxsng,thesng)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Gridded data are interpolated to single-level locations to
! produce observation differences. For now, bilinear interpolation
! is used to create a vertical sounding at the obs location
! and that is interpolated in the vertical to obtain the
! bacground value.
!
!-----------------------------------------------------------------------
!
! AUTHOR:
!
! Keith Brewster, CAPS, July, 1995
!
! MODIFICATION HISTORY:
! Sept, 1996 (KB)
! Added count-reporting diagnostic output.
!
! Dec., 1998 (KB)
! Added correlation categories.
!
! Feb., 1999 (KB)
! Moved calculation of oanxsng variables to outside the IF (good)
! block, so that it will always be available, for future
! calculations. Needed for determining qv thresholding from RH
! thresholding. Added sychronization of wind error flags, if one
! component is judged to be in error, the other is also.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz,nz_tab
INTEGER :: mxsng,nvar
INTEGER :: nobsng
!
!-----------------------------------------------------------------------
!
! Arrays defining model grid
!
!-----------------------------------------------------------------------
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
!
INTEGER :: icatg(nx,ny)
!
!-----------------------------------------------------------------------
!
! Background fields
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
REAL :: qback(nvar,nz_tab)
REAL :: hqback(nz_tab)
INTEGER :: nlvqback
!
!-----------------------------------------------------------------------
!
! Column-interpolated values (used internally only)
!
!-----------------------------------------------------------------------
!
REAL :: shght(nz)
REAL :: su(nz)
REAL :: sv(nz)
REAL :: spres(nz)
REAL :: stheta(nz)
REAL :: sqv(nz)
!
!-----------------------------------------------------------------------
!
! Station variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnsng(mxsng)
INTEGER :: isrcsng(mxsng)
INTEGER :: icatsng(mxsng)
REAL :: elevsng(mxsng)
REAL :: xsng(mxsng)
REAL :: ysng(mxsng)
REAL :: obsng(nvar,mxsng)
REAL :: qobsng(nvar,mxsng)
INTEGER :: qualsng(nvar,mxsng)
REAL :: odifsng(nvar,mxsng)
REAL :: oanxsng(nvar,mxsng)
REAL :: thesng(mxsng)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,kgrid,ktab,kntdom
INTEGER :: ipt,jpt,ireturn,nlevs
REAL :: selev,hgtgnd,hgtagl
REAL :: whigh,wlow,wqhigh,wqlow,qbthe
!
!-----------------------------------------------------------------------
!
! Beginning of executable code
!
!-----------------------------------------------------------------------
!
kntdom=0
DO ista=1,nobsng
!
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL findlc
(nx,ny,xs,ys,xsng(ista),ysng(ista), &
ipt,jpt,ireturn)
!
!-----------------------------------------------------------------------
!
! Interpolated gridded data if extrapolation is not indicated.
!
!-----------------------------------------------------------------------
!
IF(ireturn == 0) THEN
! write(6,'(a,a)') stnsng(ista),
! : ' inside ARPS domain, processing'
kntdom=kntdom+1
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6)')
! : ' location x= ',(0.001*xsng(ista)),
! : ' y= ',(0.001*ysng(ista)),
! : ' found near pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2,/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
!
!-----------------------------------------------------------------------
!
! Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
CALL colint(nx,ny,nz,nvar, &
xs,ys,zp,xsng(ista),ysng(ista),ipt,jpt,anx, &
su,sv,stheta,spres,shght,sqv,selev, &
nlevs)
!
!-----------------------------------------------------------------------
!
! Assign category
!
!-----------------------------------------------------------------------
!
icatsng(ista)=icatg(ipt,jpt)
!
!-----------------------------------------------------------------------
!
! Interpolate in the vertical at each obs level and compute
! obs differences.
!
!-----------------------------------------------------------------------
!
hgtgnd=0.5*(shght(2)+shght(1))
DO kgrid=2,nlevs-2
IF(shght(kgrid) > elevsng(ista)) EXIT
END DO
! 251 CONTINUE
whigh=(elevsng(ista)-shght(kgrid-1))/ &
(shght(kgrid)-shght(kgrid-1))
wlow=1.-whigh
!
hgtagl=AMIN1(hqback(nlvqback), &
AMAX1(hqback(1),(elevsng(ista)-hgtgnd)))
DO ktab=2,nlvqback-1
IF(hqback(ktab) > hgtagl) EXIT
END DO
! 281 CONTINUE
wqhigh=(hgtagl-hqback(ktab-1))/ &
(hqback(ktab)-hqback(ktab-1))
wqlow=1.-wqhigh
!
oanxsng(1,ista)=(wlow*su(kgrid-1)+whigh*su(kgrid))
IF(qualsng(1,ista) > 0) THEN
odifsng(1,ista)=obsng(1,ista)-oanxsng(1,ista)
qobsng(1,ista)=qobsng(1,ista)/ &
(wqlow*qback(1,ktab-1)+wqhigh*qback(1,ktab))
END IF
oanxsng(2,ista)=(wlow*sv(kgrid-1)+whigh*sv(kgrid))
IF(qualsng(2,ista) > 0) THEN
odifsng(2,ista)=obsng(2,ista)-oanxsng(2,ista)
qobsng(2,ista)=qobsng(2,ista)/ &
(wqlow*qback(2,ktab-1)+wqhigh*qback(2,ktab))
END IF
oanxsng(3,ista)=EXP( wlow*LOG(spres(kgrid-1))+ &
whigh*LOG(spres(kgrid)) )
IF(qualsng(3,ista) > 0) THEN
odifsng(3,ista)=obsng(3,ista)-oanxsng(3,ista)
qobsng(3,ista)=qobsng(3,ista)/ &
(wqlow*qback(3,ktab-1)+wqhigh*qback(3,ktab))
END IF
oanxsng(4,ista)=(wlow*stheta(kgrid-1)+whigh*stheta(kgrid))
IF(qualsng(4,ista) > 0) THEN
odifsng(4,ista)=obsng(4,ista)-oanxsng(4,ista)
qbthe=wqlow*qback(4,ktab-1)+wqhigh*qback(4,ktab)
thesng(ista)=( (obsng(4,ista)/qobsng(4,ista))+ &
(oanxsng(4,ista)/qbthe) ) / &
((1./qobsng(4,ista)) + (1./qbthe))
! print *, ' Theta:anx,obs,qsrc,qback: ',
! : oanxsng(4,ista),obsng(4,ista),qobsng(4,ista),qbthe
qobsng(4,ista)=qobsng(4,ista)/qbthe
! print *, ' theta-sng qobs normalized: ',
! : thesng(ista),qobsng(4,ista)
ELSE
thesng(ista)=(wlow*stheta(kgrid-1)+whigh*stheta(kgrid))
END IF
oanxsng(5,ista)=(wlow*sqv(kgrid-1)+whigh*sqv(kgrid))
IF(qualsng(5,ista) > 0) THEN
odifsng(5,ista)=obsng(5,ista)-oanxsng(5,ista)
qobsng(5,ista)=qobsng(5,ista)/ &
(wqlow*qback(5,ktab-1)+wqhigh*qback(5,ktab))
END IF
!
ELSE
!
! write(6,'(a,a)') stnsng(ista),
! : ' outside ARPS domain, skipping'
!
isrcsng(ista)=0
icatsng(ista)=1
! write(6,'(a,a,i4)') ' Extrapolation status returned from',
! : ' subroutine findlc: ',ireturn
! IF(ireturn.eq.-1) write(6,'(a)')
! : ' Extrapolation in x detected'
! IF(ireturn.eq.-2) write(6,'(a)')
! : ' Extrapolation in y detected'
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6/)')
! : ' location x= ',(0.001*xsng(ista)),
! : ' y= ',(0.001*ysng(ista)),
! : ' nearest pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
DO ivar=1,nvar
qualsng(ivar,ista)=-9
END DO
END IF
END DO
!
WRITE(6,'(//,1x,i5,a,i5,a)') kntdom,' of ',nobsng, &
' single-level stations are inside the ARPS domain'
RETURN
END SUBROUTINE grdtosng
!
SUBROUTINE grdtoua(nx,ny,nz,nz_tab,nzua,mxua,nvar,nobsua, & 1,2
xs,ys,zp,anx,qback,hqback,nlvqback, &
shght,su,sv,spres,stheta,sqv, &
stnua,isrcua,elevua,xua,yua,hgtua, &
obsua,qobsua,qualua,nlevsua, &
odifua,oanxua,theua)
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRDTOUA ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
! PURPOSE:
! Gridded data are interpolated to upper air locations to
! produce observation differences. For now, bilinear interpolation
! is used to create a vertical sounding at the obs location
! and that is interpolated in the vertical to obtain the
! bacground value.
!
! AUTHOR:
!
! Keith Brewster, CAPS, July, 1995
!
! MODIFICATION HISTORY:
! Sept, 1996 (KB)
! Added count-reporting diagnostic output.
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz,nz_tab
INTEGER :: nzua,mxua,nvar
INTEGER :: nobsua
!
!-----------------------------------------------------------------------
!
! Arrays defining model grid
!
!-----------------------------------------------------------------------
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
!
!-----------------------------------------------------------------------
!
! Background fields
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
REAL :: qback(nvar,nz_tab)
REAL :: hqback(nz_tab)
INTEGER :: nlvqback
!
REAL :: shght(nz)
REAL :: su(nz)
REAL :: sv(nz)
REAL :: spres(nz)
REAL :: stheta(nz)
REAL :: sqv(nz)
!
!-----------------------------------------------------------------------
!
! Station variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnua(mxua)
INTEGER :: isrcua(mxua)
REAL :: elevua(mxua)
REAL :: xua(mxua),yua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
REAL :: qobsua(nvar,nzua,mxua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: nlevsua(mxua)
REAL :: odifua(nvar,nzua,mxua)
REAL :: oanxua(nvar,nzua,mxua)
REAL :: theua(nzua,mxua)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,k,kgrid,ktab,kntdom
INTEGER :: ipt,jpt,ireturn,nlevs
REAL :: selev,hgtgnd,hgtagl
REAL :: whigh,wlow,wqhigh,wqlow,qbthe
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
kntdom=0
DO ista=1,nobsua
!
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL findlc(nx,ny,xs,ys,xua(ista),yua(ista), &
ipt,jpt,ireturn)
!
!-----------------------------------------------------------------------
!
! Interpolate if extrapoltion is not indicated.
!
!-----------------------------------------------------------------------
!
IF(ireturn == 0) THEN
! write(6,'(a,a)') stnua(ista),
! : ' inside ARPS domain, processing'
kntdom=kntdom+1
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6)')
! : ' location x= ',(0.001*xua(ista)),
! : ' y= ',(0.001*yua(ista)),
! : ' found near pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2,/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
!
!-----------------------------------------------------------------------
!
! Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
CALL colint(nx,ny,nz,nvar, &
xs,ys,zp,xua(ista),yua(ista),ipt,jpt,anx, &
su,sv,stheta,spres,shght,sqv,selev, &
nlevs)
!
!-----------------------------------------------------------------------
!
! Convert pressure to ln(p) for interpolation
!
!-----------------------------------------------------------------------
!
DO k=1,nlevs
spres(k)=LOG(spres(k))
END DO
!
!-----------------------------------------------------------------------
!
! Interpolate in the vertical at each obs level and compute
! obs differences.
!
!-----------------------------------------------------------------------
!
hgtgnd=0.5*(shght(2)+shght(1))
DO k=1,nlevsua(ista)
DO kgrid=2,nlevs-1
IF(shght(kgrid) > hgtua(k,ista)) GO TO 260
END DO
DO ivar=1,nvar
qualua(ivar,k,ista)=-9
END DO
CYCLE
260 CONTINUE
whigh=(hgtua(k,ista)-shght(kgrid-1))/ &
(shght(kgrid)-shght(kgrid-1))
wlow=1.-whigh
!
hgtagl=AMIN1(hqback(nlvqback), &
AMAX1(hqback(1),(hgtua(k,ista)-hgtgnd)))
DO ktab=2,nlvqback-1
IF(hqback(ktab) > hgtagl) EXIT
END DO
! 281 CONTINUE
wqhigh=(hgtagl-hqback(ktab-1))/ &
(hqback(ktab)-hqback(ktab-1))
wqlow=1.-wqhigh
!
IF(qualua(1,k,ista) > 0) THEN
oanxua(1,k,ista)=wlow*su(kgrid-1)+whigh*su(kgrid)
odifua(1,k,ista)=obsua(1,k,ista)-oanxua(1,k,ista)
qobsua(1,k,ista)=qobsua(1,k,ista)/ &
(wqlow*qback(1,ktab-1)+wqhigh*qback(1,ktab))
END IF
IF(qualua(2,k,ista) > 0) THEN
oanxua(2,k,ista)=wlow*sv(kgrid-1)+whigh*sv(kgrid)
odifua(2,k,ista)=obsua(2,k,ista)-oanxua(2,k,ista)
qobsua(2,k,ista)=qobsua(2,k,ista)/ &
(wqlow*qback(2,ktab-1)+wqhigh*qback(2,ktab))
END IF
IF(qualua(3,k,ista) > 0) THEN
oanxua(3,k,ista)=EXP( wlow*spres(kgrid-1)+ &
whigh*spres(kgrid) )
odifua(3,k,ista)=obsua(3,k,ista)-oanxua(3,k,ista)
qobsua(3,k,ista)=qobsua(3,k,ista)/ &
(wqlow*qback(3,ktab-1)+wqhigh*qback(3,ktab))
END IF
IF(qualua(4,k,ista) > 0) THEN
oanxua(4,k,ista)=wlow*stheta(kgrid-1)+whigh*stheta(kgrid)
odifua(4,k,ista)=obsua(4,k,ista)-oanxua(4,k,ista)
qbthe=wqlow*qback(4,ktab-1)+wqhigh*qback(4,ktab)
theua(k,ista)=( (obsua(4,k,ista)/qobsua(4,k,ista))+ &
(oanxua(4,k,ista)/qbthe) ) / &
((1./qobsua(4,k,ista)) + (1./qbthe))
qobsua(4,k,ista)=qobsua(4,k,ista)/qbthe
ELSE
theua(k,ista)=wlow*stheta(kgrid-1)+whigh*stheta(kgrid)
END IF
IF(qualua(5,k,ista) > 0) THEN
oanxua(5,k,ista)=wlow*sqv(kgrid-1)+whigh*sqv(kgrid)
odifua(5,k,ista)=obsua(5,k,ista)-oanxua(5,k,ista)
qobsua(5,k,ista)=qobsua(5,k,ista)/ &
(wqlow*qback(5,ktab-1)+wqhigh*qback(5,ktab))
END IF
END DO
!
ELSE
!
! write(6,'(a,a)') stnua(ista),
! : ' outside ARPS domain, skipping'
isrcua(ista)=0
! write(6,'(a,a,i4)') ' Extrapolation status returned from',
! : ' subroutine findlc: ',ireturn
! IF(ireturn.eq.-1) write(6,'(a)')
! : ' Extrapolation in x detected'
! IF(ireturn.eq.-2) write(6,'(a)')
! : ' Extrapolation in y detected'
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6/)')
! : ' location x= ',(0.001*xua(ista)),
! : ' y= ',(0.001*yua(ista)),
! : ' nearest pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
!
DO k=1,nzua
DO ivar=1,nvar
qualua(ivar,k,ista)=-9
END DO
END DO
END IF ! inside domain?
END DO
! wdt
!WRITE(6,'(//,1x,i5,a,i5,a)') kntdom,' of ',nobsua, &
! ' upper air profiles are inside the ARPS domain'
PRINT *, kntdom,' of ',nobsua, &
' upper air profiles are inside the ARPS domain'
RETURN
END SUBROUTINE grdtoua
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRDTORET ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE grdtoret(nx,ny,nz,nztab, & 1,2
nzret,mxret,mxcolret,nvar,ncolret, &
xs,ys,zp,anx,qback,hqback,nlvqback, &
shght,su,sv,spres,stheta,sqv, &
stnret,iret,xretc,yretc,hgtretc, &
obsret,qobsret,qualret,nlevret, &
odifret,oanxret,theretc)
!-----------------------------------------------------------------------
!
! PURPOSE:
! Gridded data are interpolated to retrieval column locations to
! produce observation differences. For now, bilinear interpolation
! is used to create a vertical sounding at the obs location
! and that is interpolated in the vertical to obtain the
! bacground value.
!
! AUTHOR:
!
! Keith Brewster, CAPS, July, 1995
!
! MODIFICATION HISTORY:
! Sept, 1996 (KB)
! Added count-reporting diagnostic output.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz,nztab
INTEGER :: nzret,mxret,mxcolret,nvar
INTEGER :: ncolret
!
!-----------------------------------------------------------------------
!
! Arrays defining model grid
!
!-----------------------------------------------------------------------
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
!
!-----------------------------------------------------------------------
!
! Background fields
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
REAL :: qback(nvar,nztab)
REAL :: hqback(nztab)
INTEGER :: nlvqback
!
REAL :: shght(nz)
REAL :: su(nz)
REAL :: sv(nz)
REAL :: spres(nz)
REAL :: stheta(nz)
REAL :: sqv(nz)
!
!-----------------------------------------------------------------------
!
! Station variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnret(mxret)
INTEGER :: iret(mxcolret)
REAL :: xretc(mxcolret),yretc(mxcolret)
REAL :: hgtretc(nzret,mxcolret)
REAL :: obsret(nvar,nzret,mxcolret)
REAL :: qobsret(nvar,nzret,mxcolret)
INTEGER :: qualret(nvar,nzret,mxcolret)
INTEGER :: nlevret(mxcolret)
REAL :: odifret(nvar,nzret,mxcolret)
REAL :: oanxret(nvar,nzret,mxcolret)
REAL :: theretc(nzret,mxcolret)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,k,kgrid,ktab,kntdom
INTEGER :: ipt,jpt,ireturn,nlevs
REAL :: selev,hgtgnd,hgtagl
REAL :: whigh,wlow,wqhigh,wqlow,qbthe
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
kntdom=0
DO ista=1,ncolret
!
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL findlc(nx,ny,xs,ys,xretc(ista),yretc(ista), &
ipt,jpt,ireturn)
!
!-----------------------------------------------------------------------
!
! Interpolate if extrapoltion is not indicated.
!
!-----------------------------------------------------------------------
!
IF(ireturn == 0) THEN
! write(6,'(a,i4,a)') 'Column ',ista,
! : ' inside ARPS domain, processing'
kntdom=kntdom+1
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6)')
! : ' location x= ',(0.001*xua(ista)),
! : ' y= ',(0.001*yua(ista)),
! : ' found near pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2,/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
!
!-----------------------------------------------------------------------
!
! Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
CALL colint(nx,ny,nz,nvar, &
xs,ys,zp,xretc(ista),yretc(ista),ipt,jpt,anx, &
su,sv,stheta,spres,shght,sqv,selev, &
nlevs)
!
!-----------------------------------------------------------------------
!
! Convert pressure to ln(p) for interpolation
!
!-----------------------------------------------------------------------
!
DO k=1,nlevs
spres(k)=LOG(spres(k))
END DO
!
!-----------------------------------------------------------------------
!
! Interpolate in the vertical at each obs level and compute
! obs differences.
!
!-----------------------------------------------------------------------
!
hgtgnd=0.5*(shght(2)+shght(1))
DO k=1,nlevret(ista)
DO kgrid=2,nlevs-1
IF(shght(kgrid) > hgtretc(k,ista)) GO TO 260
END DO
DO ivar=1,nvar
qualret(ivar,k,ista)=-9
END DO
CYCLE
260 CONTINUE
whigh=(hgtretc(k,ista)-shght(kgrid-1))/ &
(shght(kgrid)-shght(kgrid-1))
wlow=1.-whigh
!
hgtagl=AMIN1(hqback(nlvqback), &
AMAX1(hqback(1),(hgtretc(k,ista)-hgtgnd)))
DO ktab=2,nlvqback-1
IF(hqback(ktab) > hgtagl) EXIT
END DO
! 281 CONTINUE
wqhigh=(hgtagl-hqback(ktab-1))/ &
(hqback(ktab)-hqback(ktab-1))
wqlow=1.-wqhigh
!
IF(qualret(1,k,ista) > 0) THEN
oanxret(1,k,ista)=wlow*su(kgrid-1)+whigh*su(kgrid)
odifret(1,k,ista)=obsret(1,k,ista)-oanxret(1,k,ista)
qobsret(1,k,ista)=qobsret(1,k,ista)/ &
(wqlow*qback(1,ktab-1)+wqhigh*qback(1,ktab))
END IF
IF(qualret(2,k,ista) > 0) THEN
oanxret(2,k,ista)=wlow*sv(kgrid-1)+whigh*sv(kgrid)
odifret(2,k,ista)=obsret(2,k,ista)-oanxret(2,k,ista)
qobsret(2,k,ista)=qobsret(2,k,ista)/ &
(wqlow*qback(2,ktab-1)+wqhigh*qback(2,ktab))
END IF
IF(qualret(3,k,ista) > 0) THEN
oanxret(3,k,ista)=EXP( wlow*spres(kgrid-1)+ &
whigh*spres(kgrid) )
odifret(3,k,ista)=obsret(3,k,ista)-oanxret(3,k,ista)
qobsret(3,k,ista)=qobsret(3,k,ista)/ &
(wqlow*qback(3,ktab-1)+wqhigh*qback(3,ktab))
END IF
IF(qualret(4,k,ista) > 0) THEN
oanxret(4,k,ista)=wlow*stheta(kgrid-1)+whigh*stheta(kgrid)
odifret(4,k,ista)=obsret(4,k,ista)-oanxret(4,k,ista)
qbthe=wqlow*qback(4,ktab-1)+wqhigh*qback(4,ktab)
theretc(k,ista)=( (obsret(4,k,ista)/qobsret(4,k,ista))+ &
(oanxret(4,k,ista)/qbthe) ) / &
((1./qobsret(4,k,ista)) + (1./qbthe))
qobsret(4,k,ista)=qobsret(4,k,ista)/qbthe
ELSE
theretc(k,ista)=wlow*stheta(kgrid-1)+whigh*stheta(kgrid)
END IF
IF(qualret(5,k,ista) > 0) THEN
oanxret(5,k,ista)=wlow*sqv(kgrid-1)+whigh*sqv(kgrid)
odifret(5,k,ista)=obsret(5,k,ista)-oanxret(5,k,ista)
qobsret(5,k,ista)=qobsret(5,k,ista)/ &
(wqlow*qback(5,ktab-1)+wqhigh*qback(5,ktab))
END IF
END DO
!
ELSE
!
! write(6,'(a,a)') stnret(ista),
! : ' outside ARPS domain, skipping'
iret(ista)=0
! write(6,'(a,a,i4)') ' Extrapolation status returned from',
! : ' subroutine findlc: ',ireturn
! IF(ireturn.eq.-1) write(6,'(a)')
! : ' Extrapolation in x detected'
! IF(ireturn.eq.-2) write(6,'(a)')
! : ' Extrapolation in y detected'
! write (6,'(2x,a,f12.2,a,f12.2,/
! : 2x,a,i6,a,i6/)')
! : ' location x= ',(0.001*xretc(ista)),
! : ' y= ',(0.001*yretc(ista)),
! : ' nearest pt i= ',ipt,' j= ',jpt
! write (6,'(12x,a,f12.2,a,f12.2/)')
! : 'x= ',(0.001*xs(ipt)),' y= ',(0.001*ys(jpt))
!
DO k=1,nzret
DO ivar=1,nvar
qualret(ivar,k,ista)=-9
END DO
END DO
END IF ! inside domain?
END DO
WRITE(6,'(//,1x,i5,a,i5,a)') kntdom,' of ',ncolret, &
' retrieval columns are inside the ARPS domain'
RETURN
END SUBROUTINE grdtoret
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COLINT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE colint(nx,ny,nz,nvar, & 5,8
xs,ys,zp,xpt,ypt,ipt,jpt,anx, &
su,sv,stheta,spres,shght,sqv,selev, &
nlevs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolates ARPS history data in the horizontal to create
! a column of data located at point xpt, ypt.
!
! Bilinear interpolation is used.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! April 1992.
!
! MODIFICATION HISTORY:
!
! October, 1992 (K. Brewster)
! Conversion to ARPS 3.0.
!
! October, 1994 (K. Brewster)
! Conversion to ARPS 4.0.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx,ny,nz Dimensions of ARPS grids.
!
! xs x coordinate of scalar points in physical/comp. space (m)
! ys y coordinate of scalar points in physical/comp. space (m)
! zp z coordinate of scalar grid points in physical space (m)
!
! xpt x coordinate of desired sounding (m)
! ypt y coordinate of desired sounding (m)
!
! ipt i index of grid point just west of xpt,ypt
! jpt j index of grid point just south of xpt,ypt
!
! anx Background field
!
! OUTPUT:
!
! su Interpolated u wind component. (m/s)
! sv Interpolated v wind component. (m/s)
! stheta Interpolated potential temperature (K).
! spres Interpolated pressure. (Pascals)
! shght Interpolated height (meters)
! sqv Interpolated specific humidity (kg/kg)
! selev Interpolated surface elevation (m)
! nlevs Number of above-ground sounding levels.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
! Arguments -- location data
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nvar ! Dimensions of ARPS grids.
REAL :: xs(nx) ! x coordinate of grid points in
! physical/comp. space (m)
REAL :: ys(ny) ! y coordinate of grid points in
! physical/comp. space (m)
REAL :: zp(nx,ny,nz) ! z coordinate of grid points in
! physical space (m)
REAL :: xpt ! location to find in x coordinate (m)
REAL :: ypt ! location to find in y coordinate (m)
INTEGER :: ipt ! i index to the west of desired
! location
INTEGER :: jpt ! j index to the south of desired
! location
!
!-----------------------------------------------------------------------
!
! Arguments -- background field
!
!-----------------------------------------------------------------------
!
REAL :: anx(nx,ny,nz,nvar)
!
!-----------------------------------------------------------------------
!
! Arguments -- Extracted sounding variables
!
!-----------------------------------------------------------------------
!
REAL :: su(nz),sv(nz),stheta(nz),sqv(nz)
REAL :: spres(nz),shght(nz)
REAL :: selev
INTEGER :: nlevs
!
!-----------------------------------------------------------------------
!
! Functions called
!
!-----------------------------------------------------------------------
!
REAL :: aint2d
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k,in,jn
REAL :: delx,ddx,dely,ddy,w1,w2,w3,w4
REAL :: t2,t3,hmid,tmid
!
!-----------------------------------------------------------------------
!
! Find corner weights
!
!-----------------------------------------------------------------------
!
in=ipt+1
delx=xs(in)-xs(ipt)
IF(ABS(delx) > 0.) THEN
ddx=(xpt-xs(ipt))/delx
ELSE
ddx=0.
END IF
jn=jpt+1
dely=ys(jn)-ys(jpt)
IF(ABS(dely) > 0.) THEN
ddy=(ypt-ys(jpt))/dely
ELSE
ddy=0.
END IF
w1=(1.-ddx)*(1.-ddy)
w2=ddx*(1.-ddy)
w3=ddx*ddy
w4=(1.-ddx)*ddy
!
!-----------------------------------------------------------------------
!
! Interpolate all variables at all levels.
!
!-----------------------------------------------------------------------
!
nlevs=nz-1
DO k=1,nz-1
shght(k)= &
aint2d(nx,ny,nz, zp,ipt,jpt,k,in,jn,w1,w2,w3,w4)
su(k)= &
aint2d(nx,ny,nz, anx(1,1,1,1),ipt,jpt,k,in,jn,w1,w2,w3,w4)
sv(k)= &
aint2d(nx,ny,nz, anx(1,1,1,2),ipt,jpt,k,in,jn,w1,w2,w3,w4)
spres(k)= &
aint2d(nx,ny,nz, anx(1,1,1,3),ipt,jpt,k,in,jn,w1,w2,w3,w4)
stheta(k)= &
aint2d(nx,ny,nz, anx(1,1,1,4),ipt,jpt,k,in,jn,w1,w2,w3,w4)
sqv(k)= &
aint2d(nx,ny,nz, anx(1,1,1,5),ipt,jpt,k,in,jn,w1,w2,w3,w4)
END DO
!
!-----------------------------------------------------------------------
!
! Get height at scalar points, since zp was defined at w points.
!
!-----------------------------------------------------------------------
!
selev=shght(2)
DO k=1,nz-1
shght(k)=0.5*(shght(k+1)+shght(k))
END DO
!
!-----------------------------------------------------------------------
!
! Get a value at the surface, by linearly interpolating
! between the 1st and second levels.
!
!-----------------------------------------------------------------------
!
w2=(selev-shght(1))/(shght(2)-shght(1))
w1=1.-w2
su(1)=w1* su(1) + w2* su(2)
sv(1)=w1* sv(1) + w2* sv(2)
stheta(1)=w1*stheta(1) + w2*stheta(2)
sqv(1)=w1* sqv(1) + w2* sqv(2)
shght(1)=selev
!
!-----------------------------------------------------------------------
!
! Integrate downward to get the pressure at level 1.
!
!-----------------------------------------------------------------------
!
t3=stheta(3)*(spres(3)/100000.)**rddcp
t2=stheta(2)*(spres(2)/100000.)**rddcp
hmid=0.5*(shght(2)+shght(1))
tmid=t3+(((shght(3)-hmid)/(shght(3)-shght(2)))*(t2-t3))
spres(1)=spres(2)*EXP(g*(shght(2)-shght(1))/(rd*tmid))
RETURN
END SUBROUTINE colint
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE QCDIFF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE qcdiff(nvar,nvarrad,nvarradin,mxsng,nsrcsng, & 1,3
nzua,mxua,nsrcua, &
nzrad,mxrad,mxcolrad,nsrcrad, &
nzret,mxret,mxcolret,nsrcret, &
nobsng,nobsua,ncolrad,ncolret,nam_var, &
stnsng,isrcsng,hgtsng,obsng,oanxsng,odifsng, &
qcthrsng,qualsng, &
stnua,isrcua,hgtua,obsua,oanxua,odifua, &
qcthrua,qualua,nlevsua, &
stnrad,irad,isrcrad,hgtradc,obsrad,odifrad, &
qcthrrad,qualrad,nlevrad, &
stnret,iret,isrcret,hgtretc, &
obsret,oanxret,odifret, &
qcthrret,qualret,nlevret, &
obsknt,rejknt)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Observation differences are compared to expected errors
! and rejected if they exceed a threshold.
!
! AUTHOR:
!
! Keith Brewster, CAPS, July, 1995
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! Sizing variables
!
INTEGER :: nvar,nvarrad,nvarradin,mxsng,nsrcsng
INTEGER :: nzua,mxua,nsrcua
INTEGER :: nzrad,mxrad,mxcolrad,nsrcrad
INTEGER :: nzret,mxret,mxcolret,nsrcret
!
INTEGER :: nobsng,nobsua,ncolrad,ncolret
CHARACTER (LEN=6) :: nam_var(nvar)
!
!-----------------------------------------------------------------------
!
! Station variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=5) :: stnsng(mxsng)
INTEGER :: isrcsng(mxsng)
REAL :: obsng(nvar,mxsng)
REAL :: hgtsng(mxsng)
REAL :: qcthrsng(nvar,nsrcsng)
INTEGER :: qualsng(nvar,mxsng)
REAL :: oanxsng(nvar,mxsng)
REAL :: odifsng(nvar,mxsng)
!
CHARACTER (LEN=5) :: stnua(mxua)
INTEGER :: isrcua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
REAL :: oanxua(nvar,nzua,mxua)
REAL :: odifua(nvar,nzua,mxua)
REAL :: qcthrua(nvar,nsrcua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: nlevsua(mxua)
!
CHARACTER (LEN=5) :: stnrad(mxrad)
INTEGER :: irad(mxcolrad)
INTEGER :: isrcrad(0:mxrad)
REAL :: hgtradc(nzrad,mxcolrad)
REAL :: obsrad(nvarradin,nzrad,mxcolrad)
REAL :: odifrad(nvarrad,nzrad,mxcolrad)
REAL :: qcthrrad(nvarradin,nsrcrad)
INTEGER :: qualrad(nvarrad,nzrad,mxcolrad)
INTEGER :: nlevrad(mxcolrad)
!
CHARACTER (LEN=5) :: stnret(mxret)
INTEGER :: iret(mxcolret)
INTEGER :: isrcret(0:mxret)
REAL :: hgtretc(nzret,mxcolret)
REAL :: obsret(nvar,nzret,mxcolret)
REAL :: oanxret(nvar,nzret,mxcolret)
REAL :: odifret(nvar,nzret,mxcolret)
REAL :: qcthrret(nvar,nsrcret)
INTEGER :: qualret(nvar,nzret,mxcolret)
INTEGER :: nlevret(mxcolret)
REAL :: obsknt(nvar)
REAL :: rejknt(nvar)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,klev,jsrc
REAL :: pctrej,vrdiff
REAL :: tk,qvsat,qcthrqv
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Check single-level observation differences
!
!-----------------------------------------------------------------------
!
DO ivar=1,nvar
obsknt(ivar)=0.
rejknt(ivar)=0.
END DO
DO ista=1,nobsng
IF(isrcsng(ista) /= 0) THEN
DO ivar=1,nvar-1
IF(qualsng(ivar,ista) > 0) THEN
obsknt(ivar)=obsknt(ivar)+1.
IF(ABS(odifsng(ivar,ista)) > qcthrsng(ivar,isrcsng(ista)) ) THEN
WRITE(6,'(a,a,a,a)') &
' Discarding ',nam_var(ivar), &
' at ',stnsng(ista)
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtsng(ista),' Ob= ',obsng(ivar,ista), &
' Diff= ',odifsng(ivar,ista), &
' Thresh=',qcthrsng(ivar,isrcsng(ista))
qualsng(ivar,ista)=-199
rejknt(ivar)=rejknt(ivar)+1.
END IF
END IF
END DO
IF(qualsng(5,ista) > 0) THEN
obsknt(5)=obsknt(5)+1.
tk=oanxsng(4,ista)*((oanxsng(3,ista)/p0)**rddcp)
qvsat=f_qvsat(oanxsng(3,ista),tk)
qcthrqv=qcthrsng(5,isrcsng(ista))*qvsat
IF(ABS(odifsng(5,ista)) > qcthrqv) THEN
WRITE(6,'(a,a,a,a)') &
' Discarding ',nam_var(5), &
' at ',stnsng(ista)
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtsng(ista),' Ob= ',obsng(5,ista), &
' Diff= ',odifsng(5,ista), &
' Thresh=',qcthrqv
qualsng(5,ista)=-199
rejknt(5)=rejknt(5)+1.
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Synchronize wind component flags
!
!-----------------------------------------------------------------------
!
IF(qualsng(1,ista) == -199) qualsng(2,ista)=-199
IF(qualsng(2,ista) == -199) qualsng(1,ista)=-199
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Write statistics
!
!-----------------------------------------------------------------------
!
WRITE(6,'(//a/a)') ' QC Statistics for single-level data', &
' Variable Obs Rejected Percent'
DO ivar=1,nvar
pctrej=0.
IF(obsknt(ivar) > 0.) pctrej=100.*rejknt(ivar)/obsknt(ivar)
WRITE(6,'(2x,a,i10,i10,f8.1)') &
nam_var(ivar),nint(obsknt(ivar)),nint(rejknt(ivar)),pctrej
END DO
WRITE(6,'(/a)') ' '
!
!-----------------------------------------------------------------------
!
! Check upper-air observation differences
!
!-----------------------------------------------------------------------
!
DO ivar=1,nvar
obsknt(ivar)=0.
rejknt(ivar)=0.
END DO
DO ista=1,nobsua
IF(isrcua(ista) /= 0) THEN
DO klev=1,nlevsua(ista)
DO ivar=1,nvar-1
IF(qualua(ivar,klev,ista) > 0) THEN
obsknt(ivar)=obsknt(ivar)+1.
IF(ABS(odifua(ivar,klev,ista)) > qcthrua(ivar,isrcua(ista)) ) THEN
WRITE(6,'(a,a,a,a,a,i4)') ' Discarding UA ', &
nam_var(ivar),' at sta ',stnua(ista), &
' klev =',klev
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtua(klev,ista), &
' Ob= ',obsua(ivar,klev,ista), &
' Diff= ',odifua(ivar,klev,ista), &
' Thresh=',qcthrua(ivar,isrcua(ista))
qualua(ivar,klev,ista)=-199
rejknt(ivar)=rejknt(ivar)+1.
END IF
END IF
END DO
IF(qualua(5,klev,ista) > 0) THEN
obsknt(5)=obsknt(5)+1.
tk=oanxua(4,klev,ista)* &
((oanxua(3,klev,ista)/p0)**rddcp)
qvsat=f_qvsat(oanxua(3,klev,ista),tk)
qcthrqv=qcthrua(5,isrcua(ista))*qvsat
IF(ABS(odifua(5,klev,ista)) > qcthrqv) THEN
WRITE(6,'(a,a,a,a)') &
' Discarding ',nam_var(5), &
' at ',stnua(ista)
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtua(klev,ista), &
' Ob= ',obsua(5,klev,ista), &
' Diff= ',odifua(5,klev,ista), &
' Thresh=',qcthrqv
qualua(5,klev,ista)=-199
rejknt(5)=rejknt(5)+1.
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Synchronize wind component flags
!
!-----------------------------------------------------------------------
!
IF(qualua(1,klev,ista) == -199) qualua(2,klev,ista)=-199
IF(qualua(2,klev,ista) == -199) qualua(1,klev,ista)=-199
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Write statistics
!
!-----------------------------------------------------------------------
!
WRITE(6,'(//a/a)') ' QC Statistics for UA data', &
' Variable Obs Rejected Percent'
DO ivar=1,nvar
pctrej=0.
IF(obsknt(ivar) > 0.) pctrej=100.*rejknt(ivar)/obsknt(ivar)
WRITE(6,'(2x,a,i10,i10,f8.1)') &
nam_var(ivar),nint(obsknt(ivar)),nint(rejknt(ivar)),pctrej
END DO
WRITE(6,'(/a)') ' '
!
!-----------------------------------------------------------------------
!
! Check radar data observation differences
! Here radial wind difference is computed and compared only.
!
!-----------------------------------------------------------------------
!
obsknt(1)=0.
rejknt(1)=0.
DO ista=1,ncolrad
IF(irad(ista) /= 0) THEN
jsrc=isrcrad(irad(ista))
DO klev=1,nlevrad(ista)
IF(qualrad(1,klev,ista) > 0 .AND. qualrad(2,klev,ista) > 0 ) THEN
obsknt(1)=obsknt(1)+1.
vrdiff=SQRT( &
odifrad(1,klev,ista)*odifrad(1,klev,ista)+ &
odifrad(2,klev,ista)*odifrad(2,klev,ista))
IF(vrdiff > qcthrrad(2,jsrc) ) THEN
WRITE(6,'(a,a,i5,a,i4)') ' Discarding radar Vr', &
' at col ',ista,' klev =',klev
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtradc(klev,ista), &
' Ob= ',obsrad(2,klev,ista), &
' Diff= ',vrdiff, &
' Thresh=',qcthrrad(2,jsrc)
qualrad(1,klev,ista)=-199
qualrad(2,klev,ista)=-199
rejknt(1)=rejknt(1)+1.
END IF
END IF
END DO
END IF
END DO
WRITE(6,'(//a/a)') ' QC Statistics for radar data', &
' Variable Obs Rejected Percent'
pctrej=0.
IF(obsknt(1) > 0.) pctrej=100.*rejknt(1)/obsknt(1)
WRITE(6,'(2x,a,i10,i10,f8.1)') &
'Vr',nint(obsknt(1)),nint(rejknt(1)),pctrej
WRITE(6,'(/a)') ' '
!
!-----------------------------------------------------------------------
!
! Check retrieval observation differences
!
!-----------------------------------------------------------------------
!
DO ivar=1,nvar
obsknt(ivar)=0.
rejknt(ivar)=0.
END DO
DO ista=1,ncolret
IF(iret(ista) /= 0) THEN
jsrc=isrcret(iret(ista))
DO klev=1,nlevret(ista)
DO ivar=1,nvar-1
IF(qualret(ivar,klev,ista) > 0) THEN
obsknt(ivar)=obsknt(ivar)+1.
IF(ABS(odifret(ivar,klev,ista)) > qcthrret(ivar,jsrc) ) THEN
WRITE(6,'(a,a,a,i5,a,i4)') ' Discarding retrieval ', &
nam_var(ivar),' at col ',ista,' klev =',klev
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtretc(klev,ista), &
' Ob= ',obsret(ivar,klev,ista), &
' Diff= ',odifret(ivar,klev,ista), &
' Thresh=',qcthrret(ivar,jsrc)
qualret(ivar,klev,ista)=-199
rejknt(ivar)=rejknt(ivar)+1.
END IF
END IF
END DO
IF(qualret(5,klev,ista) > 0) THEN
obsknt(5)=obsknt(5)+1.
tk=oanxret(4,klev,ista)* &
((oanxret(3,klev,ista)/p0)**rddcp)
qvsat=f_qvsat(oanxret(3,klev,ista),tk)
qcthrqv=qcthrret(5,jsrc)*qvsat
IF(ABS(odifret(5,klev,ista)) > qcthrqv) THEN
WRITE(6,'(a,a,a,i5,a,i4)') ' Discarding retrieval ', &
nam_var(5),' at col ',ista,' klev =',klev
WRITE(6,'(a,f8.0,a,f16.6,a,f16.6,a,f16.6)') &
' Hgt= ',hgtretc(klev,ista), &
' Ob= ',obsret(5,klev,ista), &
' Diff= ',odifret(5,klev,ista), &
' Thresh=',qcthrqv
qualret(5,klev,ista)=-199
rejknt(5)=rejknt(5)+1.
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Synchronize wind component flags
!
!-----------------------------------------------------------------------
!
IF(qualret(1,klev,ista) == -199) qualret(2,klev,ista)=-199
IF(qualret(2,klev,ista) == -199) qualret(1,klev,ista)=-199
END DO
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Write statistics
!
!-----------------------------------------------------------------------
!
WRITE(6,'(//a/a)') ' QC Statistics for retrieval data', &
' Variable Obs Rejected Percent'
DO ivar=1,nvar
pctrej=0.
IF(obsknt(ivar) > 0.) pctrej=100.*rejknt(ivar)/obsknt(ivar)
WRITE(6,'(2x,a,i10,i10,f8.1)') &
nam_var(ivar),nint(obsknt(ivar)),nint(rejknt(ivar)),pctrej
END DO
WRITE(6,'(/a)') ' '
RETURN
END SUBROUTINE qcdiff