!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UNFNQC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE UNFNQC(maxgate,maxazim, ngate,nazim, &,3
lvlprof,zsnd,usnd,vsnd, &
gsp_vel,rfrst_vel, &
rlat_ptr,rlon_ptr,ralt_ptr, &
rvel,spw,elev,azim,vnyq, &
unfvel,tmp1,tmp2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Unfold and Quality Control Doppler radial velocities.
!
! AUTHOR:
! Keith Brewster, Center for Analysis and Prediction of Storms
!
! MODIFICATION HISTORY:
! Nov. 2000, First Version
! Oct. 2001, Updated and added spectrum width and std dev thresholding
!
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum number of gates in a radial
! maxazim Maximum number of radials in a tilt
! ngate Number of gates in radial
! nazim Number of radials
! rvel Doppler radial velocity
! spw Doppler spectrum width
! elev Elevation angle
! azim Azimuth
! vnyq Nyquist velocity
! unfvel Unfolded Doppler radial velocity
!
! WORK ARRAYS:
!
! tmp1
! tmp2
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
INTEGER, INTENT(IN) :: lvlprof
INTEGER, INTENT(IN) :: gsp_vel
INTEGER, INTENT(IN) :: rfrst_vel
REAL, INTENT(IN) :: rvel(maxgate,maxazim)
REAL, INTENT(IN) :: spw(maxgate,maxazim)
REAL, INTENT(IN) :: elev(maxazim)
REAL, INTENT(IN) :: azim(maxazim)
REAL, INTENT(IN) :: vnyq
REAL, INTENT(IN) :: rlat_ptr
REAL, INTENT(IN) :: rlon_ptr
REAL, INTENT(IN) :: ralt_ptr
REAL, INTENT(IN) :: zsnd(lvlprof)
REAL, INTENT(IN) :: usnd(lvlprof)
REAL, INTENT(IN) :: vsnd(lvlprof)
REAL, INTENT(OUT) :: unfvel(maxgate,maxazim)
REAL, INTENT(OUT) :: tmp1(maxgate,maxazim)
REAL, INTENT(OUT) :: tmp2(maxgate,maxazim)
INTEGER :: iray,kray
INTEGER :: bgate,igate,kgate,kkgate,bkgate,ekgate,lgate,kgate1,kgate2
INTEGER :: egate,k,kclose
INTEGER :: kntall,kntgood,kntfold,kntrej
INTEGER :: istatwrt,rem_points,found
REAL :: bknt,fknt,tknt,tpoint
REAL :: velok,twonyq,inv2nyq,thrpri
REAL :: sum,sum2,vavg,sdev,thresh,tstdev,tstvel,sum1,vavg1
REAL :: rlen,valid_vel,height,diff,refvel
!
!-----------------------------------------------------------------------
!
! Constants
!
!-----------------------------------------------------------------------
!
REAL :: rvchek = -90.
REAL :: spwflag = -931.
CHARACTER (LEN=6) :: varid
CHARACTER (LEN=20):: varname
CHARACTER (LEN=20):: varunits = 'm/s'
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
kntall=0
kntgood=0
kntfold=0
kntrej=0
velok = -(vnyq+1.0E-05)
twonyq = 2.0*vnyq
inv2nyq = 1./twonyq
thrpri = min(max((0.5*vnyq),10.),vnyq)
!
!-----------------------------------------------------------------------
!
! Apply Spectrum width thresholding
!
!-----------------------------------------------------------------------
!
CALL SPWTHR
(maxgate,maxazim, &
ngate,nazim,rvchek,spwflag,vnyq,rvel,spw)
!
!-----------------------------------------------------------------------
!
! Set up tmp1 and tmp2 to be quality arrays. 1=good 0=bad/missing
!
!-----------------------------------------------------------------------
!
tmp1=0.
tmp2=0.
DO iray=1,nazim
DO igate=1,ngate
unfvel(igate,iray)=rvel(igate,iray)
IF( rvel(igate,iray) > velok) THEN
tmp1(igate,iray)=1.
tmp2(igate,iray)=1.
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Apply gate-to-gate shear checking for velocity folding.
! Follows Eilts and Smith, 1989, JTech, 118.
!
!-----------------------------------------------------------------------
!
DO iray=1,nazim
bgate=ngate
DO igate=1,ngate-1
IF( tmp1(igate,iray) > 0.) THEN
bgate=igate
EXIT
END IF
END DO
IF( bgate < ngate ) THEN
! DO igate=bgate+1,ngate
DO igate=bgate,ngate
IF( tmp1(igate,iray) > 0. ) THEN
kntall=kntall+1
lgate = 0
ekgate = max(1,igate-5)
bkgate = max(1,igate-1)
DO kgate = igate-1,ekgate,-1
IF(tmp2(kgate,iray) >0. ) THEN
lgate = kgate
EXIT
ENDIF
ENDDO
IF(lgate >0) THEN
refvel = unfvel(lgate,iray)
ELSE
refvel = 999.
ENDIF
IF( abs(rvel(igate,iray) - refvel) > thrpri) THEN
sum=0.
sum2=0.
fknt=0.
bknt=0.
bkgate=max((igate-10),bgate)
DO kgate=igate-1,bkgate,-1
bknt=bknt+tmp2(kgate,iray)
sum=sum+tmp2(kgate,iray)*unfvel(kgate,iray)
sum2=sum2+tmp2(kgate,iray)*(unfvel(kgate,iray)*unfvel(kgate,iray))
IF (bknt > 3. ) EXIT
END DO
IF( iray > 1 ) THEN
kray=iray-1
ekgate=min((igate+10),ngate)
DO kgate=igate,ekgate
fknt=fknt+tmp2(kgate,kray)
sum=sum+tmp2(kgate,kray)*unfvel(kgate,kray)
sum2=sum2+tmp2(kgate,kray)*(unfvel(kgate,kray)*unfvel(kgate,kray))
IF ( fknt > 4. ) EXIT
END DO
IF( fknt < 5. .AND. kray > 1 ) THEN
kray=kray-1
DO kgate=igate,ekgate
fknt=fknt+tmp2(kgate,kray)
sum=sum+tmp2(kgate,kray)*unfvel(kgate,kray)
sum2=sum2+tmp2(kgate,kray)*(unfvel(kgate,kray)*unfvel(kgate,kray))
IF ( fknt > 4. ) EXIT
END DO
END IF
END IF
tknt=bknt+fknt
IF(tknt > 0. ) THEN
vavg=sum/tknt
thresh=max((0.4*abs(vavg)),thrpri)
IF ( tknt > 2. ) THEN
sdev=sqrt((sum2-(sum*sum/tknt))/(tknt-1.))
thresh=max((2.0*sdev),thresh)
END IF
IF( abs(rvel(igate,iray)-vavg) < thresh ) THEN
kntgood=kntgood+1
ELSE
tstdev=twonyq*NINT((vavg-rvel(igate,iray))*inv2nyq)
tstvel=rvel(igate,iray)+tstdev
IF( abs(tstvel-vavg) < thresh ) THEN
IF(abs(tstdev) > 0.) THEN
kntfold=kntfold+1
! write(6,'(a,f10.1,a,f10.1,a,f10.1)') &
! ' Unfold: Meas=',rvel(igate,iray),' Avgvel=', &
! vavg,' New=',tstvel
unfvel(igate,iray)=tstvel
ELSE
kntgood=kntgood+1
END IF
ELSE
kntrej=kntrej+1
! write(8,'(a,i5,a,i5,a,f10.1,a,f10.1,a,f10.1,a,f10.1)') &
! 'igate=',igate,'iray=',iray, &
! ' Reject: Meas=',rvel(igate,iray),' Avgvel=', &
! vavg,' Try=',tstvel,'vnyq=',vnyq
unfvel(igate,iray)=-777.
tmp2(igate,iray)=0.
ENDIF ! tstvel-vavg >thresh
END IF ! rvel-vavg >thresh
ELSE ! tknt =0.
!
!-----------------------------------------------------------------------
!
! If no points are found to comparison,use environmental wind as a reference
!
!-----------------------------------------------------------------------
!
IF(initopt == 2 .or. initopt ==3)THEN !no valid vel,use snd
rlen = rfrst_vel + (igate-1)*gsp_vel
height = rlen*sin(elev(iray)*3.1415926/180.)
diff = 500.0
DO k=1,lvlprof
IF( abs(zsnd(k)-height)<diff ) THEN
diff = abs(zsnd(k)-height)
kclose = k
ENDIF
ENDDO
CALL uv2vr(rlen,elev(iray),azim(iray), &
rlat_ptr,rlon_ptr,ralt_ptr, &
usnd(kclose),vsnd(kclose),valid_vel)
! IF(abs(rvel(igate,iray)-valid_vel) < 1.5*thrpri) THEN
IF(abs(rvel(igate,iray)-valid_vel) < thrpri) THEN
kntgood = kntgood+1
ELSE
IF(valid_vel*rvel(igate,iray)<0)THEN
tstdev=twonyq*NINT(valid_vel/abs(valid_vel))
ELSE
tstdev=twonyq*NINT((valid_vel-rvel(igate,iray))*inv2nyq)
ENDIF
tstvel=rvel(igate,iray)+tstdev
! IF( abs(tstvel-valid_vel) < 1.5*thrpri ) THEN
IF( abs(tstvel-valid_vel) < 3.0*thrpri ) THEN
IF(abs(tstdev) > 0.) THEN
kntfold=kntfold+1
ELSE
kntgood = kntgood+1
END IF
unfvel(igate,iray)=tstvel
tmp2(igate,iray) = 1.
ELSE
kntrej=kntrej+1
unfvel(igate,iray)=-777.
tmp2(igate,iray)=0.
ENDIF !tstvel-valid_vel
ENDIF !rval-valid_vel<1.5
ELSE
kntrej=kntrej+1
unfvel(igate,iray)=-777.
tmp2(igate,iray)=0.
ENDIF !initopt=1 analytical
ENDIF ! tknt >0.
ELSE
kntgood=kntgood+1
END IF !(rvel(igate,iray) - refvel) > thrpri
END IF ! tmp1(igate,iray) > 0.
END DO
END IF ! bgate < ngate
END DO
!
!-----------------------------------------------------------------------
!
! Second pass to deal with those failed points which have been removed
! during 1st pass. Use average of 15 points in 3 radials
!
!-----------------------------------------------------------------------
!
DO igate = 2,ngate
DO iray = 2,nazim-1
IF(unfvel(igate,iray) == -777. .AND. tmp2(igate,iray) == 0.) THEN
sum=0.
tpoint=0.
sum2=0.
bkgate=max((igate-10),1)
ekgate=min((igate+10),ngate)
DO kray = iray-1,iray+1
fknt = 0.
tknt=0.
DO kgate=igate,ekgate
tknt=tknt+tmp2(kgate,kray)
sum=sum+tmp2(kgate,kray)*unfvel(kgate,kray)
sum2=sum2+tmp2(kgate,kray)*unfvel(kgate,kray)*unfvel(kgate,kray)
IF(tknt > 3.) EXIT
ENDDO
DO kgate=igate-1,bkgate,-1
fknt=fknt+tmp2(kgate,kray)
sum=sum+tmp2(kgate,kray)*unfvel(kgate,kray)
sum2=sum2+tmp2(kgate,kray)*unfvel(kgate,kray)*unfvel(kgate,kray)
IF(fknt > 2.) EXIT
ENDDO
tpoint = tpoint + tknt + fknt
ENDDO
IF(tpoint > 0.) THEN
vavg = sum/tpoint
thresh=max((0.4*abs(vavg)),thrpri)
IF ( tpoint > 2. ) THEN
sdev=sqrt((sum2-(sum*sum/tpoint))/(tpoint-1.))
thresh=max((2.0*sdev),thresh)
END IF
IF( abs(rvel(igate,iray)-vavg) < thresh ) THEN
kntgood=kntgood+1
unfvel(igate,iray)=rvel(igate,iray)
tmp2(igate,iray)=1.
ELSE
tstdev=twonyq*NINT((vavg-rvel(igate,iray))*inv2nyq)
tstvel=rvel(igate,iray)+tstdev
IF( abs(tstvel-vavg) < thresh ) THEN
IF(abs(tstdev) > 0.) THEN
kntfold=kntfold+1
unfvel(igate,iray)=tstvel
tmp2(igate,iray)=1.
ELSE
kntgood=kntgood+1
unfvel(igate,iray)=rvel(igate,iray)
tmp2(igate,iray)=1.
ENDIF
ENDIF
ENDIF
ELSE !no valid vel,use snd
!-----------------------------------------------------------------------
!
! If no points are found to comparison,use environmental wind as a reference
!
!-----------------------------------------------------------------------
!
IF(initopt ==2 .or. initopt ==3)THEN
rlen = rfrst_vel + (igate-1)*gsp_vel
height = rlen*sin(elev(iray)*3.1415926/180.)
diff = 500.0
DO k=1,lvlprof
IF( abs(zsnd(k)-height)<diff ) THEN
diff = abs(zsnd(k)-height)
kclose = k
ENDIF
ENDDO
CALL uv2vr(rlen,elev(iray),azim(iray), &
rlat_ptr,rlon_ptr,ralt_ptr, &
usnd(kclose),vsnd(kclose),valid_vel)
! IF(abs(rvel(igate,iray)-valid_vel) < 1.5*thrpri) THEN
IF(abs(rvel(igate,iray)-valid_vel) < thrpri) THEN
kntgood = kntgood+1
ELSE
IF(valid_vel*rvel(igate,iray)<0)THEN
tstdev=twonyq*NINT(valid_vel/abs(valid_vel))
ELSE
tstdev=twonyq*NINT((valid_vel-rvel(igate,iray))*inv2nyq)
ENDIF
tstvel=rvel(igate,iray)+tstdev
! IF( abs(tstvel-valid_vel) < 1.5*thrpri ) THEN
IF( abs(tstvel-valid_vel) < 3*thrpri ) THEN
IF(abs(tstdev) > 0.) THEN
kntfold=kntfold+1
ELSE
kntgood = kntgood+1
END IF
unfvel(igate,iray)=tstvel
tmp2(igate,iray) = 1.
ELSE
kntrej=kntrej+1
unfvel(igate,iray)=-777.
tmp2(igate,iray)=0.
ENDIF !tstvel-valid_vel
ENDIF !rval-valid_vel<1.5
ELSE
kntrej=kntrej+1
unfvel(igate,iray)=-777.
tmp2(igate,iray)=0.
ENDIF ! initopt = 1
ENDIF ! tknt >0.
ENDIF ! -777.0
ENDDO
ENDDO
write(6,'(a/,4i12)') &
' Gates Good-Gates Folded Rejected', &
kntall,kntgood,kntfold,kntrej
RETURN
END SUBROUTINE unfnqc
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DESPEKL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE DESPEKL(maxgate,maxazim, & 3
ngate,nazim,varchek,rdata,tmp)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Despeckle radar data.
! Can be used for reflectivity or velocity.
!
!
! AUTHOR: Keith Brewster
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
! maxgate Maximum number of gates in a radial
! maxazim Maximum number of radials in a tilt
! ngate Number of gates in radial
! nazim Number of radials
! varchek Threshold to determine good data vs. flagged
! rdata Doppler radial velocity
!
! OUTPUT:
! rdata
!
! WORK ARRAYS:
!
! tmp
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
REAL, INTENT(IN) :: varchek
REAL, INTENT(INOUT) :: rdata(maxgate,maxazim)
REAL, INTENT(OUT) :: tmp(maxgate,maxazim)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
INTEGER :: kntgd,kntdsp
!
REAL :: sum,pctdsp
REAL, PARAMETER :: dspmiss = -991.
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Create temporary array where tmp=1 if non-missing tmp=0 if missing
!
!-----------------------------------------------------------------------
!
!
tmp=0.
DO j=1,nazim
DO i=1,ngate
IF ( rdata(i,j) > varchek ) tmp(i,j)=1.
END DO
END DO
!
!-----------------------------------------------------------------------
!
! If datum has fewer than 3 non-missing neighbors in a 3x3
! square template set it to missing
!
!-----------------------------------------------------------------------
!
kntgd=0
kntdsp=0
DO j=2,nazim-1
DO i=2,ngate-1
IF(rdata(i,j) > varchek ) THEN
kntgd=kntgd+1
sum=tmp(i-1,j+1)+tmp(i,j+1)+tmp(i+1,j+1) &
+tmp(i-1,j )+ tmp(i+1,j) &
+tmp(i-1,j-1)+tmp(i,j-1)+tmp(i+1,j-1)
IF( sum < 3. ) THEN
kntdsp=kntdsp+1
rdata(i,j) = dspmiss
END IF
END IF
END DO
END DO
IF(kntgd > 0 ) THEN
pctdsp=100.*(float(kntdsp)/float(kntgd))
ELSE
pctdsp=0.
END IF
write(6,'(a,i8,a,i8,a,f6.1,a)') &
' Despeckled ',kntdsp,' of ',kntgd,' data =',pctdsp,' percent.'
!
RETURN
!
END SUBROUTINE despekl
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SPWTHR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE SPWTHR(maxgate,maxazim, & 1
ngate,nazim,rvchek,spwflag,vnyq,rvel,spw)
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Threshold radial velocity data based on spectrum width.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum number of gates in a radial
! maxazim Maximum number of radials in a tilt
! ngate Number of gates in radial
! nazim Number of radials
! rvchek Threshold for good data vs. flagged
! spwflag Flag for spectrum width threshold fail
! vnyq Nyquist velocity
! rvel Doppler radial velocity
! spw Doppler spectrum width
!
! OUTPUT:
! rvel Doppler radial velocity
!
!
! WORK ARRAYS:
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
REAL, INTENT(IN) :: rvchek
REAL, INTENT(IN) :: spwflag
REAL, INTENT(IN) :: vnyq
REAL, INTENT(INOUT) :: rvel(maxgate,maxazim)
REAL, INTENT(IN) :: spw(maxgate,maxazim)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: spthrat = 0.67
REAL :: spthr
INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! If spectrum width is greater than threshold, set to missing
!
!-----------------------------------------------------------------------
!
spthr=spthrat*vnyq
DO j=1,nazim
DO i=1,ngate
IF(rvel(i,j) > rvchek .AND. spw(i,j) > spthr ) rvel(i,j) = spwflag
END DO
END DO
!
RETURN
!
END SUBROUTINE spwthr