! !################################################################## !################################################################## !###### ###### !###### 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