!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE APDETECT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE apdetect(nrefgates,nvelgates,maxazim,maxelev, & 1,1
kntrgate,kntrazim,kntrelev, &
kntvgate,kntvazim,kntvelev, &
refchek,velchek, &
irefgatsp,ivelgatsp, &
winszrad,winszazim,i_vcp, &
rngrvol,azmrvol,elvrvol, &
rngvvol,azmvvol,elvvvol, &
refvol,velvol,tmp, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take data from a radar volume in polar coordinates and map it to
! ARPS Cartesian terrain-following grid. Uses a least squares
! local quadratic fit to remap the data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: LeiLei Wang, CAPS
!
! MODIFICATION HISTORY:
! Keith Brewster,
! Modified to remove additional array allocation
!
! 07/08/2003, Keith Brewster
! Added setting of istatus return values.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
!
! dazlim Maximum value of azimuth difference (grid vs data) to accept
! Generally should be 30 degrees or less for velocity, 360 for refl
! rngmin Minimum range (m) of data to use
! (10 000 m or more to eliminate near field ground targets).
! rngmax Maximum range (m) of data to use
!
! rngvvol Range to gate in velocity 3-D volume
! azmvvol Azimuth angle in velocity 3-D volume
! elvvvol Elevation angle in velocity 3-D volume
! varvol Radar data 3-D volume
!
! xs x coordinate of scalar grid points in physical/comp. space (m)
! ys y coordinate of scalar grid points in physical/comp. space (m)
! zps Vertical coordinate of scalar grid points in physical space(m)
!
! OUTPUT:
! varvol Radar data 3-D volume
! istatus Status indicator
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nrefgates
INTEGER, INTENT(IN) :: nvelgates
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: kntrgate(maxazim,maxelev)
INTEGER, INTENT(IN) :: kntrazim(maxelev)
INTEGER, INTENT(IN) :: kntrelev
INTEGER, INTENT(IN) :: kntvgate(maxazim,maxelev)
INTEGER, INTENT(IN) :: kntvazim(maxelev)
INTEGER, INTENT(IN) :: kntvelev
INTEGER, INTENT(IN) :: i_vcp
REAL, INTENT(IN) :: winszrad
REAL, INTENT(IN) :: winszazim
REAL, INTENT(IN) :: refchek
REAL, INTENT(IN) :: velchek
INTEGER, INTENT(IN) :: irefgatsp
INTEGER, INTENT(IN) :: ivelgatsp
REAL, INTENT(IN) :: rngrvol(nrefgates,maxelev)
REAL, INTENT(IN) :: azmrvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvrvol(maxazim,maxelev)
REAL, INTENT(IN) :: rngvvol(nvelgates,maxelev)
REAL, INTENT(IN) :: azmvvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvvvol(maxazim,maxelev)
REAL, INTENT(INOUT) :: refvol(nrefgates,maxazim,maxelev)
REAL, INTENT(INOUT) :: velvol(nvelgates,maxazim,maxelev)
REAL, INTENT(INOUT) :: tmp(nrefgates,maxazim)
INTEGER, INTENT(OUT) :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. Local Variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ii,jj,kk,kkv,kv,kr2,i,j,k,knt
INTEGER :: mvelok,spin,flags
INTEGER :: igate,ivgate,jazim,kelev,jazmin,jazmax
INTEGER :: iigate,jjazim,jazmref,jazmvel
INTEGER :: irngmin,irngmax
INTEGER :: igatebgn,igateend,jazmbgn,jazmend,indgate,indazm
INTEGER :: igspratio,iwinszrad,iwinszazm
INTEGER :: kntchek,kntapref,kntapvel
REAL :: azmdiff,maxdbz
REAL :: summdbz,sumvel,sumvel2,sumtdbz,sumtvel
REAL :: dbzdiff,sign,refprev,prevvel,veldiff
REAL :: all_counts,spinchange_counts,signcnt
REAL :: delev,avgelv,avgelvv,avgelvr,avgelvr2,appct
REAL :: mdbz,tdbz,deltdbz,spinchange,stdvel,mvel,tvel
LOGICAL :: found
REAL, PARAMETER :: dbzthr = 10.0
REAL, PARAMETER :: apflag = -888.0
REAL, PARAMETER :: spinthr = 15.0
REAL, PARAMETER :: ddbzthr = -12.0
REAL, PARAMETER :: mvelthr = 2.3
REAL, PARAMETER :: apvelthr = 5.0
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Misc initializations
!
!-----------------------------------------------------------------------
!
istatus=0
IF( kntrelev > 1 ) THEN
igspratio=max(int(float(irefgatsp)/float(2*ivelgatsp)),1)
iwinszrad=max(int(winszrad/float(2*irefgatsp)),1)
iwinszazm=max(int(0.5*winszazim),1)
!
!-----------------------------------------------------------------------
!
! Establish the index of the reflectivity tilt above 1.1 degrees.
! For NEXRAD the first two unique tilt angles should be 0.5
! and 1.5 degrees.
!
!-----------------------------------------------------------------------
!
found=.false.
avgelvr2=0.
DO kk=2,kntrelev
avgelvr2=0.
DO jazim=1,kntrazim(kk)
avgelvr2=avgelvr2+elvrvol(jazim,kk)
END DO
avgelvr2=avgelvr2/float(kntrazim(kk))
IF(avgelvr2 > 1.1) THEN
found=.true.
EXIT
END IF
END DO
IF(.NOT. found) THEN
write(6,'(a,f9.2,/a//)') ' Highest elev not above 1.1 deg: ', &
avgelvr2,' Skipping AP Detect'
istatus=-2
RETURN
END IF
kr2=kk
write(6,'(a,f9.2,a,i6)') ' Elev 1st ref tilt above 1.1 deg: ', &
avgelvr2,' kr2=',kr2
!
!-----------------------------------------------------------------------
!
! For now, AP detection is only done on any tilts below 1.1 degree.
!
!-----------------------------------------------------------------------
!
DO kk = 1, (kr2-1)
tmp=0.
kntchek=0
kntapref=0
kntapvel=0
!-----------------------------------------------------------------------
!
! Establish the velocity tilt that most closely matches the
! this reflectivity elevation angle.
!
!-----------------------------------------------------------------------
!
avgelvr=0.
DO jazim=1,kntrazim(kk)
avgelvr=avgelvr+elvrvol(jazim,kk)
END DO
avgelvr=avgelvr/float(kntrazim(kk))
kv=1
delev=999.
avgelvv=-99.
DO kkv=1,kntvelev
avgelv=0.
DO jazim=1,kntvazim(kkv)
avgelv=avgelv+elvvvol(jazim,kkv)
END DO
avgelv=avgelv/float(kntvazim(kkv))
IF(abs(avgelv-avgelvr) < delev ) THEN
delev=abs(avgelv-avgelvr)
avgelvv=avgelv
kv=kkv
END IF
END DO
write(6,'(a,f9.2,a,i6)') ' Elev of nearest velocity tilt: ', &
avgelvv,' kv=',kv
DO jjazim = 1,kntrazim(kk)
!
! Find nearest azimuth in reflectivity data above 1.1
!
jazmref=1
azmdiff = 999.0
DO jazim = 1,kntrazim(kr2)
IF(abs(azmrvol(jazim,kr2)-azmrvol(jjazim,kk))<azmdiff)THEN
azmdiff = abs(azmrvol(jazim,kr2)-azmrvol(jjazim,kk))
jazmref = jazim
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Find nearest azimuth in velocity data
!
!-----------------------------------------------------------------------
!
jazmvel=1
azmdiff = 999.0
DO jazim = 1,kntvazim(kv)
IF(abs(azmvvol(jazim,kv)-azmrvol(jjazim,kk))<azmdiff &
.and. kntvgate(jazim,kv)>0 )THEN
azmdiff = abs(azmvvol(jazim,kv)-azmrvol(jjazim,kk))
jazmvel = jazim
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Calculate mean ref, texture of ref, spin
!
!-----------------------------------------------------------------------
!
DO iigate= 1,kntrgate(jjazim,kk)
IF(refvol(iigate,jjazim,kk) > refchek )THEN
kntchek=kntchek+1
summdbz=0.
sumtdbz=0.
all_counts=0.
spinchange_counts=0.
spin = 1
irngmin = max((iigate - iwinszrad),1)
irngmax = min((iigate + iwinszrad),kntrgate(jjazim,kk))
jazmin = jjazim - iwinszazm
jazmax = jjazim + iwinszazm
IF(jazmin <= 0) jazmin = jazmin + kntrazim(kk)-2
IF(jazmax > kntrazim(kk) )jazmax = jazmax-kntrazim(kk)+2
!
!-----------------------------------------------------------------------
!
! Loop forward from jazmin to jazmax
!
!-----------------------------------------------------------------------
!
IF(jazmax > jazmin)THEN
DO jazim = jazmin,jazmax
DO igate = irngmin,irngmax
IF(refvol(igate,jazim,kk)>refchek)THEN
summdbz=summdbz+refvol(igate,jazim,kk)
sumtdbz=sumtdbz+dbzdiff*dbzdiff
all_counts = all_counts + 1
dbzdiff = refvol(igate,jazim,kk)-refprev
IF(dbzdiff > dbzthr)THEN
IF(spin<0)THEN
spinchange_counts = SPINchange_counts + 1
spin = 1
END IF
END IF
IF(dbzdiff < -dbzthr)THEN
IF(spin>0)THEN
spinchange_counts = SPINchange_counts + 1
spin = -1
END IF
END IF
IF(dbzdiff > 0.)THEN
signcnt = signcnt + 1.
ELSE
signcnt = signcnt - 1.
END IF
refprev = refvol(igate,jazim,kk)
END IF
END DO ! igate
END DO ! jazim
!
!-----------------------------------------------------------------------
!
! if jazmax<kazmin
!
!-----------------------------------------------------------------------
!
ELSE
DO jazim = jazmin,kntrazim(kk)
DO igate = irngmin,irngmax
IF(refvol(igate,jazim,kk)>refchek)THEN
summdbz=summdbz+refvol(igate,jazim,kk)
sumtdbz=sumtdbz+dbzdiff*dbzdiff
all_counts = all_counts + 1
dbzdiff = refvol(igate,jazim,kk)-refprev
IF(dbzdiff > dbzthr)THEN
IF(spin<0)THEN
spinchange_counts = SPINchange_counts + 1
spin = 1
END IF
END IF
IF(dbzdiff < -dbzthr)THEN
IF(spin>0)THEN
spinchange_counts = SPINchange_counts + 1
spin = -1
END IF
END IF
IF(dbzdiff > 0.)THEN
signcnt = signcnt + 1.
ELSE
signcnt = signcnt - 1.
END IF
refprev = refvol(igate,jazim,kk)
END IF
END DO ! igate
END DO ! jazim
DO jazim = 1,jazmax
DO igate = irngmin,irngmax
IF(refvol(igate,jazim,kk)>refchek)THEN
summdbz=summdbz+refvol(igate,jazim,kk)
sumtdbz=sumtdbz+dbzdiff*dbzdiff
all_counts = all_counts + 1
dbzdiff = refvol(igate,jazim,kk)-refprev
IF(dbzdiff > dbzthr)THEN
IF(spin<0)THEN
spinchange_counts = SPINchange_counts + 1
spin = 1
END IF
END IF
IF(dbzdiff < -dbzthr)THEN
IF(spin>0)THEN
spinchange_counts = SPINchange_counts + 1
spin = -1
END IF
END IF
IF(dbzdiff > 0.)THEN
signcnt = signcnt + 1.
ELSE
signcnt = signcnt - 1.
END IF
refprev = refvol(igate,jazim,kk)
END IF
END DO ! igate
END DO ! jazim
END IF
IF(all_counts > 0.)THEN
mdbz= summdbz/all_counts
tdbz= sumtdbz/all_counts
spinchange= (spinchange_counts/all_counts)*100.
END IF
!
! Calculate difference in reflectivity between this gate and the
! first level above 1.1 degrees.
!
IF(iigate <= kntrgate(jazmref,kr2) .and. &
refvol(iigate,jazmref,kr2) > refchek ) THEN
deltdbz= &
refvol(iigate,jazmref,kr2)-refvol(iigate,jjazim,kk)
ELSE
IF(i_vcp==31 .or. i_vcp==32) THEN
deltdbz= -32.0-refvol(iigate,jjazim,kk)
ELSE
deltdbz= -5.0-refvol(iigate,jjazim,kk)
END IF
END IF
!
! If the delta-dBZ check fails when comparing to the nearest gate,
! account for tilt of echoes by finding max of gates in the neighborhood.
!
IF (deltdbz < ddbzthr) THEN
jazmbgn = max((jazmref-1),1)
jazmend = min((jazmref+1),kntrazim(kr2))
maxdbz = refvol(iigate,jazmref,kr2)
DO indazm = jazmbgn,jazmend
igatebgn = max((iigate-1),1)
igateend = min((iigate+1),kntrgate(indazm,kr2))
DO indgate = igatebgn,igateend
maxdbz = max(maxdbz,refvol(indgate,indazm,kr2))
END DO
END DO
IF( maxdbz > refchek ) &
deltdbz = maxdbz - refvol(iigate,jjazim,kk)
END IF
!
!-----------------------------------------------------------------------
!
! Find std of vel, texture of vel, mean vel
!
! First, find nearest velocity gate
!
! Nearest velocity tilt (kv) and velocity azimuth (jazmvel) were
! found earlier
!
!-----------------------------------------------------------------------
!
DO igate=2,kntvgate(jazmvel,kv)-1
IF(rngvvol(igate,kv) >= rngrvol(iigate,kk)) EXIT
END DO
IF(abs(rngvvol(igate,kv)-rngrvol(iigate,kk)) < &
abs(rngvvol(igate-1,kv)-rngrvol(iigate,kk)) ) THEN
ivgate=igate
ELSE
ivgate=igate-1
END IF
IF(abs(rngvvol(ivgate,kv)-rngrvol(iigate,kk))< &
4*max(irefgatsp,ivelgatsp) .and. &
abs(azmvvol(jazmvel,kv)-azmrvol(jjazim,kk))< &
float(iwinszazm))THEN
sumvel2=0.
sumvel =0.
sumtvel=0.
mvelok=0
irngmin = max((ivgate-igspratio),1)
irngmax = min((ivgate+igspratio),kntvgate(jazmvel,kv))
jazmin = jazmvel - iwinszazm
jazmax = jazmvel + iwinszazm
IF( jazmin < 1 ) jazmin = jazmin+kntvazim(kv)
IF( jazmax > kntvazim(kv) ) jazmax = jazmax-kntvazim(kv)
!-----------------------------------------------------------------------
!
! Loop forward from jazmin to jazmax
!
!-----------------------------------------------------------------------
!
IF(jazmax > jazmin)THEN
DO jazim = jazmin,jazmax
prevvel = velvol(irngmin,jazim,kv)
DO igate = irngmin+1,irngmax
IF(velvol(igate,jazim,kv)>velchek)THEN
sumvel=sumvel+velvol(igate,jazim,kv)
sumvel2=sumvel2+ &
velvol(igate,jazim,kv)*velvol(igate,jazim,kv)
mvelok=mvelok+1
veldiff=velvol(igate,jazim,kv) - prevvel
sumtvel=sumtvel+veldiff*veldiff
prevvel=velvol(igate,jazim,kv)
END IF
END DO ! igate
END DO ! jazim
!
!-----------------------------------------------------------------------
!
! if jazmax<kazmin
!
!-----------------------------------------------------------------------
!
ELSE
DO jazim = jazmin,kntvazim(kv)
prevvel = velvol(irngmin,jazim,kv)
DO igate = irngmin+1,irngmax
IF(velvol(igate,jazim,kv)>velchek)THEN
sumvel=sumvel+velvol(igate,jazim,kv)
sumvel2=sumvel2+ &
velvol(igate,jazim,kv)*velvol(igate,jazim,kv)
mvelok=mvelok+1
veldiff=velvol(igate,jazim,kv) - prevvel
sumtvel=sumtvel+veldiff*veldiff
prevvel=velvol(igate,jazim,kv)
END IF
END DO ! igate
END DO ! jazim
DO jazim = 1,jazmax
prevvel = velvol(irngmin,jazim,kv)
DO igate = irngmin,irngmax
IF(velvol(igate,jazim,kv)>velchek)THEN
sumvel=sumvel+velvol(igate,jazim,kv)
sumvel2=sumvel2+ &
velvol(igate,jazim,kv)*velvol(igate,jazim,kv)
mvelok=mvelok+1
veldiff= velvol(igate,jazim,kv) - prevvel
sumtvel=sumtvel+veldiff*veldiff
prevvel = velvol(igate,jazim,kv)
END IF
END DO ! igate
END DO ! jazim
END IF
IF(mvelok>1)THEN
stdvel = &
sqrt((sumvel2-sumvel*sumvel/mvelok)/(mvelok-1))
tvel= sumtvel/mvelok
mvel= abs(sumvel/mvelok)
ELSE
stdvel=999.
tvel=999.
mvel=999.
END IF
ELSE ! if ivgate is too far away from iigate
stdvel=999.
tvel=999.
mvel=999.
END IF
!
! If calculated values match two or more characteristics,
! set the AP marker, tmp=1. The actual resetting of refvol
! is deferred until all points have been checked on this level,
! so that the calculation of spin at subsequent gates is not affected.
!
flags=0
IF(spinchange >= spinthr) flags=flags+1
IF(deltdbz <= ddbzthr) flags=flags+1
IF(abs(mvel) < mvelthr ) flags=flags+1
IF(flags > 1 ) tmp(iigate,jjazim)=1.0
END IF ! a valid reflectivity at iigate,jjazim,kk
END DO !iigate
END DO !jjazim
!
! Where the AP marker has been set, set for reflectivity to apflag,
! and set corresponding velocities to apflag if the absolute
! value of the velocity is less than apvelthr.
!
DO jjazim=1,kntrazim(kk)
!
jazmvel=1
azmdiff = 999.0
DO jazim = 1,kntvazim(kv)
IF(abs(azmvvol(jazim,kv)-azmrvol(jjazim,kk))<azmdiff)THEN
azmdiff = abs(azmvvol(jazim,kv)-azmrvol(jjazim,kk))
jazmvel = jazim
END IF
END DO
!
DO iigate=1,kntrgate(jjazim,kk)
IF(tmp(iigate,jjazim) > 0. ) THEN
kntapref=kntapref+1
refvol(iigate,jjazim,kk) = apflag
DO igate=2,kntvgate(jazmvel,kv)-1
IF(rngvvol(igate,kv) >= rngrvol(iigate,kk)) EXIT
END DO
IF(abs(rngvvol(igate,kv)-rngrvol(iigate,kk)) < &
abs(rngvvol(igate-1,kv)-rngrvol(iigate,kk)) ) THEN
ivgate=igate
ELSE
ivgate=igate-1
END IF
IF(abs(rngvvol(ivgate,kv)-rngrvol(iigate,kk))< &
4*max(irefgatsp,ivelgatsp) .and. &
abs(azmvvol(jazmvel,kv)-azmrvol(jjazim,kk))< &
float(iwinszazm))THEN
irngmin = ivgate - int(irefgatsp/float(2*ivelgatsp))
irngmin = max(irngmin,1)
irngmax = ivgate + int(irefgatsp/float(2*ivelgatsp))
irngmax = min(irngmax,kntvgate(jazmvel,kv))
DO igate=irngmin,irngmax
IF(abs(velvol(igate,jazmvel,kv)) < apvelthr ) THEN
kntapvel=kntapvel+1
velvol(igate,jazmvel,kv)=apflag
END IF
END DO
ENDIF
END IF
END DO ! iigate loop
END DO ! jjazim loop
IF ( kntchek > 0 ) THEN
appct=100.*float(kntapref)/float(kntchek)
ELSE
appct=0.
END IF
write(6,'(a,i6,a,f6.2,/a,i8,/a,i8,f9.2,a,/a,i8)') &
' AP detect completed for level ',kk,' elev=',avgelvr, &
' Reflectivity gates checked:',kntchek, &
' AP flagged ref:',kntapref,appct,' percent', &
' AP flagged vel:',kntapvel
!
! Apply despekl to the edited data
! This should help catch AP residue.
!
CALL despekl
(nrefgates,maxazim,nrefgates,maxazim,refchek, &
refvol(1,1,kk),tmp)
END DO ! kk
! open(21,file='ap.dat',form='unformatted',status='unknown')
! write(21) nrefgates
! write(21) maxazim
! write(21) maxelev
! write(21) kntrelev
! write(21) kntrazim
! write(21) kntrgate
! write(21) rngrvol
! write(21) azmvvol
! write(21) elvrvol
! write(21) refvol
! close(21)
ELSE
write(6,'(a)') ' Need at least two reflectivity levels for AP detection'
write(6,'(a)') ' Skipping apdetect'
istatus=-1
END IF
RETURN
END SUBROUTINE apdetect