! 
!##################################################################
!##################################################################
!######                                                      ######
!######                 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
!
!-----------------------------------------------------------------------
!
!  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 :: dbzthr,azmdiff,maxdbz
  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

  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
!
!-----------------------------------------------------------------------
!
  igspratio=max(int(float(irefgatsp)/float(2*ivelgatsp)),1)
  iwinszrad=max(int(winszrad/float(2*irefgatsp)),1)
  iwinszazm=max(int(0.5*winszazim),1)

  istatus = 0
!
!-----------------------------------------------------------------------
!
!  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.
!
!-----------------------------------------------------------------------
!
  avgelvr2=0.
  DO kk=2,kntrelev-1
    avgelvr2=0.
    DO jazim=1,kntrazim(kk)
      avgelvr2=avgelvr2+elvrvol(jazim,kk)
    END DO
    avgelvr2=avgelvr2/float(kntrazim(kk))
    IF(avgelvr2 > 1.1) EXIT
  END DO
  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)

  RETURN
  END SUBROUTINE APDETECT