!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SET_LBCOPT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE set_lbcopt(lbcnew)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! An interface to access bndry.inc commons from c routines.
! Allows program to reset lbcopt to save memory in call
! to initgrdvar when it is used for purposes other than
! initializing the model.
!
! AUTHOR:
! Keith Brewster, CAPS
! Jan 30, 2002
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: lbcnew
INCLUDE 'bndry.inc'
lbcopt=lbcnew
END SUBROUTINE set_lbcopt
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTENVPRF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE extenvprf(nx,ny,nz,lvlprof, &,3
x,y,zp,xs,ys,zps, &
u,v,ptprt,pprt,qv,ptbar,pbar,usc,vsc,rfrct, &
radarx,radary,radarz,radius, &
zsnd,ktsnd,usnd,vsnd,rfrctsnd)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Finds environmental profile around radar location from gridded data.
!
! AUTHOR:
! Keith Brewster, CAPS
! Jan 30, 2002
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx,ny,nz
INTEGER, INTENT(IN) :: lvlprof
REAL, INTENT(IN) :: x(nx)
REAL, INTENT(IN) :: y(ny)
REAL, INTENT(IN) :: zp(nx,ny,nz)
REAL, INTENT(OUT) :: xs(nx)
REAL, INTENT(OUT) :: ys(ny)
REAL, INTENT(OUT) :: zps(nx,ny,nz)
REAL, INTENT(IN) :: u(nx,ny,nz)
REAL, INTENT(IN) :: v(nx,ny,nz)
REAL, INTENT(IN) :: ptprt(nx,ny,nz)
REAL, INTENT(IN) :: pprt(nx,ny,nz)
REAL, INTENT(IN) :: qv(nx,ny,nz)
REAL, INTENT(IN) :: ptbar(nx,ny,nz)
REAL, INTENT(IN) :: pbar(nx,ny,nz)
REAL, INTENT(OUT) :: usc(nx,ny,nz)
REAL, INTENT(OUT) :: vsc(nx,ny,nz)
REAL, INTENT(OUT) :: rfrct(nx,ny,nz)
REAL, INTENT(IN) :: radarx
REAL, INTENT(IN) :: radary
REAL, INTENT(IN) :: radarz
REAL, INTENT(IN) :: radius
REAL, INTENT(IN) :: zsnd(lvlprof)
REAL, INTENT(OUT) :: ktsnd(lvlprof)
REAL, INTENT(OUT) :: usnd(lvlprof)
REAL, INTENT(OUT) :: vsnd(lvlprof)
REAL, INTENT(OUT) :: rfrctsnd(lvlprof)
INTEGER :: i,j,k,kbot,ktop,kext
INTEGER :: iradar,jradar
INTEGER :: ibeg,iend,jbeg,jend
INTEGER :: knt,knttot
REAL :: dx,dy
REAL :: radius2,dist2,whigh,wlow,accept,rfrgrad
dx=x(2)-x(1)
dy=y(2)-y(1)
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
xs(nx)=xs(nx-1)+dx
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
ys(ny)=ys(ny-1)+dy
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zps(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
zps(i,j,nz)=2.0*zp(i,j,nz-1)-zp(i,j,nz-2)
END DO
END DO
!-----------------------------------------------------------------------
!
! Bring u and v to a common grid, the scalar points.
!
!-----------------------------------------------------------------------
CALL avgx
(u, 0, nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
usc)
CALL avgy(v, 0, nx,ny,nz, &
1,nx-1,1,ny-1,1,nz-1, &
vsc)
!
!-----------------------------------------------------------------------
!
! Find the index of scalar grid location that is closest to the radar.
! Set bounds of searching to be within a box of
! radarx-radius < xs < radarx+radius
! radary-radius < ys < radary+radius
!
!-----------------------------------------------------------------------
!
iradar=1+nint((radarx-xs(1))/dx)
jradar=1+nint((radary-ys(1))/dy)
ibeg=iradar-(int(radius/dx)+1)
ibeg=min(ibeg,iradar,nx-2)
ibeg=max(ibeg,2)
iend=iradar+(int(radius/dx)+1)
iend=max(iend,iradar,2)
iend=min(iend,nx-2)
jbeg=jradar-(int(radius/dy)+1)
jbeg=min(jbeg,jradar,ny-2)
jbeg=max(jbeg,2)
jend=jradar+(int(radius/dy)+1)
jend=max(jend,jradar,2)
jend=min(jend,ny-2)
iradar=min(max(iradar,2),nx-1)
jradar=min(max(jradar,2),ny-1)
!
!-----------------------------------------------------------------------
!
! Initialize sums to zero
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
ktsnd(k)=0.
usnd(k)=0.
vsnd(k)=0.
rfrctsnd(k)=0.
END DO
!
CALL refract(nx,ny,nz,ptprt,pprt,qv,ptbar,pbar,rfrct)
!
knt=0
knttot=0
radius2=radius*radius
DO j=jbeg,jend
DO i=ibeg,iend
!
!-----------------------------------------------------------------------
!
! Is this point within the ARPS domain?
! Since the ARPS grid is Cartesian, need only compare
! the external grid coordinates to the x and y limits
! of the ARPS grid.
!
!-----------------------------------------------------------------------
!
knttot=knttot+1
dist2 = (xs(i)-radarx)*(xs(i)-radarx) &
+(ys(j)-radary)*(ys(j)-radary)
IF(dist2 < radius2 .OR. &
(i == iradar .AND. j == jradar)) THEN
knt=knt+1
!
!-----------------------------------------------------------------------
!
!
! Interpolate external data in vertical onto profile
! arrays.
!
!-----------------------------------------------------------------------
!
DO k=1,lvlprof
IF(zps(i,j,2) <= zsnd(k)) THEN
DO kext=3,nz-1
IF(zps(i,j,kext) >= zsnd(k)) EXIT
END DO
IF(kext > nz-1) EXIT
whigh=(zsnd(k)-zps(i,j,kext-1))/ &
(zps(i,j,kext)-zps(i,j,kext-1))
wlow=1.-whigh
ktsnd(k)=ktsnd(k)+1.
usnd(k)=usnd(k)+ &
whigh*usc(i,j,kext)+wlow*usc(i,j,kext-1)
vsnd(k)=vsnd(k)+ &
whigh*vsc(i,j,kext)+wlow*vsc(i,j,kext-1)
rfrctsnd(k)=rfrctsnd(k)+ &
whigh*rfrct(i,j,kext)+wlow*rfrct(i,j,kext-1)
END IF
END DO
END IF
END DO
END DO
WRITE(6,'(/a,i6,a,/a,i6,a/)') &
' extenvprf found ',knt,' points within radius', &
' of ',knttot,' checked.'
accept=0.3*float(knt)
!
!-----------------------------------------------------------------------
!
! Find lowest height with data
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a)') ' Finding range of mean profile data ...'
DO k=1,lvlprof
IF(ktsnd(k) > accept) EXIT
END DO
kbot=k
!
!-----------------------------------------------------------------------
!
! Find highest height with data
!
!-----------------------------------------------------------------------
!
DO k=lvlprof,2,-1
WRITE(6,'(a,f10.2,a,f6.0,a,f10.0)') ' z = ',zsnd(k), &
' knt = ',ktsnd(k),' accept = ',accept
IF(ktsnd(k) > accept) EXIT
END DO
ktop=k
!
WRITE(6,'(a,f10.2,a,f10.2,a)') &
' Height of data for wind profile spans from ', &
zsnd(kbot),' to ',zsnd(ktop),' meters.'
!
!-----------------------------------------------------------------------
!
! Divide through to find average
!
!-----------------------------------------------------------------------
!
DO k=kbot,ktop
usnd(k)=usnd(k)/ktsnd(k)
vsnd(k)=vsnd(k)/ktsnd(k)
rfrctsnd(k)=rfrctsnd(k)/ktsnd(k)
END DO
!
!-----------------------------------------------------------------------
!
! Set variables "below-ground"
! Zero gradiant assumed in usnd, vsnd.
! Constant gradiant assumed in rfrsnd.
!
!-----------------------------------------------------------------------
!
rfrgrad=(rfrctsnd(kbot+1)-rfrctsnd(kbot))/(zsnd(kbot+1)-zsnd(kbot))
DO k=kbot-1,1,-1
usnd(k)=usnd(kbot)
vsnd(k)=vsnd(kbot)
rfrctsnd(k)=rfrctsnd(kbot)+(zsnd(k)-zsnd(kbot))*rfrgrad
END DO
!
!-----------------------------------------------------------------------
!
! Set variables "above-top"
! Zero gradiant assumed.
!
!-----------------------------------------------------------------------
!
rfrgrad=(rfrctsnd(ktop)-rfrctsnd(ktop-1))/(zsnd(ktop)-zsnd(ktop-1))
DO k=ktop+1,lvlprof
usnd(k)=usnd(ktop)
vsnd(k)=vsnd(ktop)
rfrctsnd(k)=rfrctsnd(ktop)+(zsnd(k)-zsnd(ktop))*rfrgrad
END DO
RETURN
END SUBROUTINE extenvprf
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UV2VR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE uv2vr(rlen,elev, azimu, & 2,6
rlat_rad,rlon_rad,ralt_rad,ugrid,vgrid,vr)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Project u,v to radial direction
!
!-----------------------------------------------------------------------
!
! INPUT:
! rlen distance to radar location
! elev elevation angle
! azimu azimuth angle
! rlat_rad latitude of radar
! rlon_rad latitude of radar
! ralt_rad altitude of radar
! ugrid u
! vgrid v
!
! OUTPUT:
! vr radial wind
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INCLUDE FILES
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
REAL :: rlen, elev, azimu
REAL :: rlat_rad,rlon_rad,ralt_rad
REAL :: ugrid,vgrid
REAL :: vr
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: azim,dist,dz,eleva,range,dhdr,dsdr,ddrot
REAL :: xrad, yrad, r_horiz, z_verti
REAL :: xgrid,ygrid,zgrid
REAL :: rlat_grid,rlon_grid
REAL :: uazmrad,vazmrad
REAL :: DPI
PARAMETER ( DPI=3.1415926535898/180.0 )
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Get x and y locations of each radar ob
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,rlat_rad,rlon_rad, xrad, yrad)
r_horiz = rlen*COS(elev*DPI)
z_verti = rlen*SIN(elev*DPI)
xgrid = r_horiz*SIN( azimu*DPI ) + xrad
ygrid = r_horiz*COS( azimu*DPI ) + yrad
zgrid = z_verti
CALL xytoll(1,1,xgrid, ygrid, rlat_grid,rlon_grid)
!
!-----------------------------------------------------------------------
!
! Find heading and distance from radar along ground
! Store correlation of azimuth with u and v drections
! in uazmrad and vazmrad, respectively.
!
!-----------------------------------------------------------------------
!
CALL disthead(rlat_grid,rlon_grid, &
rlat_rad,rlon_rad, &
azim,dist)
!
CALL ddrotuv(1,rlon_grid, azim,1.,ddrot, &
uazmrad,vazmrad)
!
!-----------------------------------------------------------------------
!
! Loop in height
!
!-----------------------------------------------------------------------
!
dz= zgrid - ralt_rad
CALL beamelv(dz,dist,eleva,range)
CALL dhdrange(eleva,range,dhdr)
dsdr=SQRT(AMAX1(0.,(1.-dhdr*dhdr)))
!
!-----------------------------------------------------------------------
!
! Project u-v to radial direction to get radial velocity
!
!-----------------------------------------------------------------------
!
vr=(uazmrad*ugrid + vazmrad*vgrid) * dsdr
!
RETURN
END SUBROUTINE uv2vr
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE UV2VRN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE uv2vrn(nzsnd,zsnd,rfrsnd,rfropt,rlen,elev,azimu, &,6
rlat_rad,rlon_rad,ralt_rad,ugrid,vgrid,vr)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Project u,v to radial direction
!
!-----------------------------------------------------------------------
!
! INPUT:
! rlen distance to radar location
! elev elevation angle
! azimu azimuth angle
! rlat_rad latitude of radar
! rlon_rad latitude of radar
! ralt_rad altitude of radar
! ugrid u
! vgrid v
!
! OUTPUT:
! vr radial wind
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INCLUDE FILES
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
INTEGER, INTENT(IN) :: nzsnd
REAL, INTENT(IN) :: zsnd(nzsnd)
REAL, INTENT(IN) :: rfrsnd(nzsnd)
INTEGER, INTENT(IN) :: rfropt
REAL, INTENT(IN) :: rlen, elev, azimu
REAL, INTENT(IN) :: rlat_rad,rlon_rad,ralt_rad
REAL, INTENT(IN) :: ugrid,vgrid
REAL, INTENT(OUT) :: vr
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: azim,dist,dz,eleva,range,dhdr,dsdr,ddrot
REAL :: xrad, yrad, r_horiz, z_verti
REAL :: xgrid,ygrid,zgrid
REAL :: rlat_grid,rlon_grid
REAL :: uazmrad,vazmrad
REAL :: DPI
PARAMETER ( DPI=3.1415926535898/180.0 )
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Get x and y locations of each radar ob
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,rlat_rad,rlon_rad, xrad, yrad)
r_horiz = rlen*COS(elev*DPI)
z_verti = rlen*SIN(elev*DPI)
xgrid = r_horiz*SIN( azimu*DPI ) + xrad
ygrid = r_horiz*COS( azimu*DPI ) + yrad
zgrid = z_verti
CALL xytoll(1,1,xgrid, ygrid, rlat_grid,rlon_grid)
!
!-----------------------------------------------------------------------
!
! Find heading and distance from radar along ground
! Store correlation of azimuth with u and v drections
! in uazmrad and vazmrad, respectively.
!
!-----------------------------------------------------------------------
!
CALL disthead(rlat_grid,rlon_grid, &
rlat_rad,rlon_rad, &
azim,dist)
!
CALL ddrotuv(1,rlon_grid, azim,1.,ddrot, &
uazmrad,vazmrad)
!
!-----------------------------------------------------------------------
!
! Loop in height
!
!-----------------------------------------------------------------------
!
dz= zgrid - ralt_rad
CALL beamelvn(nzsnd,zsnd,rfrsnd,ralt_rad,rfropt,dz,dist,eleva,range)
CALL dhdrangn(nzsnd,zsnd,rfrsnd,ralt_rad,rfropt,eleva,range,dhdr)
dsdr=SQRT(AMAX1(0.,(1.-dhdr*dhdr)))
!
!-----------------------------------------------------------------------
!
! Project u-v to radial direction to get radial velocity
!
!-----------------------------------------------------------------------
!
vr=(uazmrad*ugrid + vazmrad*vgrid) * dsdr
!
RETURN
END SUBROUTINE uv2vrn
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE QUADFIT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE quadfit(maxgate,maxazim,maxelev,nzsnd,maxsort, &,3
varchek,vmedlim,dazlim,iorder,rfropt, &
kntgate,kntazim,kntelev, &
radarx,radary,rdralt,dazim, &
rngmin,rngmax,zsnd,rfrsnd, &
rngvol,azmvol,elvvol, &
rxvol,ryvol,rzvol, &
varsort,elvmean,varvol, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Uses a least squares local quadratic fit to refill any data region
! rejected during unfolding process
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
!
! MODIFICATION HISTORY:
!
! 14-Jan-2004 Keith Brewster
! Modified logic to avoid possibility of overflow of varsort.
! Updated azimuth searching to match design of remapvol.
! varsort and elvmean arrays are now allocated outside and
! passed into this routine.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
! nzsnd Number of levels in sounding arrays
! maxsort Number of elements in sorting vector for computing median
!
! varchek Threshold for checking data, good vs. flagged
! vmedlim Threshold limit for median check
! dazlim Maximum value of azimuth difference (grid vs data) to accept
! Generally should be 30 degrees or less for velocity, 360 for refl
! iorder Order of polynomial to fit (1: linear, 2: quadratic)
! rfropt Rafractivity option (1: std atmos lapse, 2: avg actual lapse rate)
! kntgate Number of gates
! kntazim Number of azimuths for each tilt
! kntelev Number of elevation angles (tilts) of data
!
! radarx x-location of radar
! radary y-location of radar
! rdralt Elevation (m MSL) of radar
! dazim Approximate range resolution of data (degrees of azimuth)
! 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
!
! zsnd Heights of levels in refractivity profile
! rfrsnd Refractivity profile
!
! 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
!
! rxvol x coordinate of radar volume elements
! ryvol y coordinate of radar volume elements
! rzvol z coordinate of radar volume elements
!
!
! OUTPUT:
!
! varsort Temporary array used for computing the mean
! elvmean Temporary array average elevation angle for each tilt
!
! varvol Radar data 3-D volume
! istatus Status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: nzsnd
INTEGER, INTENT(IN) :: maxsort
REAL, INTENT(IN) :: varchek
REAL, INTENT(IN) :: vmedlim
REAL, INTENT(IN) :: dazlim
INTEGER, INTENT(IN) :: iorder
INTEGER, INTENT(IN) :: rfropt
INTEGER, INTENT(IN) :: kntgate(maxazim,maxelev)
INTEGER, INTENT(IN) :: kntazim(maxelev)
INTEGER, INTENT(IN) :: kntelev
REAL, INTENT(IN) :: radarx
REAL, INTENT(IN) :: radary
REAL, INTENT(IN) :: rdralt
REAL, INTENT(IN) :: dazim
REAL, INTENT(IN) :: rngmin
REAL, INTENT(IN) :: rngmax
REAL, INTENT(IN) :: zsnd(nzsnd)
REAL, INTENT(IN) :: rfrsnd(nzsnd)
REAL, INTENT(IN) :: rngvol(maxgate,maxelev)
REAL, INTENT(IN) :: azmvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvvol(maxazim,maxelev)
REAL, INTENT(OUT) :: rxvol(maxgate,maxazim,maxelev)
REAL, INTENT(OUT) :: ryvol(maxgate,maxazim,maxelev)
REAL, INTENT(OUT) :: rzvol(maxgate,maxazim,maxelev)
REAL, INTENT(OUT) :: varsort(maxsort)
REAL, INTENT(OUT) :: elvmean(maxelev)
REAL, INTENT(INOUT) :: varvol(maxgate,maxazim,maxelev)
INTEGER, INTENT(OUT) :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. Local Variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: n = 6
REAL, PARAMETER :: eps = 1.0E-25
REAL :: avar(n,n)
REAL :: rhsvar(n)
REAL :: avel(n,n)
REAL :: rhsvel(n)
REAL :: sol(n)
REAL :: work(n,n+1)
REAL :: work1d(n+1)
REAL :: array(4,4)
REAL :: rhsv(4)
REAL :: solv(4)
INTEGER :: ii,jj,kk,i,j,k,knt,kinbox
INTEGER :: kok,isort,jsort,mid,maxkok
INTEGER :: kbgn,kend
INTEGER :: igate,jazim,kelev,jazmin,jmirror,jend,jknt
INTEGER :: iigate,jjazim,kkelev
INTEGER :: istatal,istatwrt
INTEGER :: nbeam,nrang
!
! jsrchmn is the minimum number of radials to search to find data
! that are within the distance limits of the grid center.
!
INTEGER, PARAMETER :: jsrchmn = 2
REAL :: deg2rad,rad2deg
REAL :: elevmin,elevmax
REAL :: delx,dely,delz,dazimr,daz,azdiff
REAL :: ddx,ddxy,ddx2,ddy,ddy2,dxthr,dxthr0
REAL :: cosaz,sinaz,sfcr,zagl
REAL :: sum,sum2,sdev,thresh,slrange,elijk,azimijk,time
REAL :: varmax,varmin,varavg,varmean,varmed
REAL :: xs,ys,zps,fitvalue
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
deg2rad = atan(1.)/45.
rad2deg = 1./deg2rad
dazimr = dazim*deg2rad
dxthr0=5*250.
dxthr=dxthr0
maxkok=0
!
time=0.
elevmax=-999.
elevmin=999.
DO kelev=1,kntelev
sum=0.
DO jazim=1,kntazim(kelev)
sum=sum+elvvol(jazim,kelev)
elevmin=min(elevmin,elvvol(jazim,kelev))
elevmax=max(elevmax,elvvol(jazim,kelev))
END DO
elvmean(kelev)=sum/float(kntazim(kelev))
END DO
!
DO kelev=1,kntelev
DO jazim=1,kntazim(kelev)
cosaz=cos(deg2rad*azmvol(jazim,kelev))
sinaz=sin(deg2rad*azmvol(jazim,kelev))
DO igate=1,kntgate(jazim,kelev)
CALL beamhgtn(nzsnd,zsnd,rfrsnd,rdralt,rfropt, &
elvvol(jazim,kelev),rngvol(igate,kelev),zagl,sfcr)
rxvol(igate,jazim,kelev)=radarx+sinaz*sfcr
ryvol(igate,jazim,kelev)=radary+cosaz*sfcr
rzvol(igate,jazim,kelev)=rdralt+zagl
END DO
END DO
END DO
!
DO kkelev = 1,kntelev
DO jjazim = 1,kntazim(kkelev)
DO iigate= 1,kntgate(jjazim,kkelev)
IF(varvol(iigate,jjazim,kkelev)==-777.0)THEN
xs = rxvol(iigate,jjazim,kkelev)
ys = ryvol(iigate,jjazim,kkelev)
zps = rzvol(iigate,jjazim,kkelev)
kok=0
sum=0.
sum2=0.
sdev=0.
varavg=999999.
varsort=999999.
delx=xs-radarx
dely=ys-radary
delz=zps-rdralt
slrange=sqrt(delx*delx+dely*dely+delz*delz)
!
IF(delz >= 0. .AND. slrange > 0. ) THEN
elijk=rad2deg*asin(delz/slrange)
IF(elijk >= elevmin .AND. elijk <= elevmax) THEN
!
!-----------------------------------------------------------------------
!
! Determine azimuth to this grid cell
!
!-----------------------------------------------------------------------
!
IF(delx == 0.) THEN
IF(dely >= 0.) THEN
azimijk=0.
ELSE
azimijk=180.
END IF
ELSE
azimijk=rad2deg*atan(delx/dely)
IF(dely < 0.) azimijk=azimijk+180.
IF(azimijk < 0.) azimijk=azimijk+360.
END IF
DO kelev=2,kntelev-1
IF(elvmean(kelev) > elijk) EXIT
END DO
!
varmax=-999.
varmin=999.
DO jj=1,n
DO ii=1,n
avar(ii,jj)=0.
END DO
END DO
!
DO ii=1,n
rhsvar(ii)=0.
END DO
!
kbgn=max((kelev-1),1)
kend=min(kelev,kntelev)
DO kk=kend,kend
!
!-----------------------------------------------------------------------
!
! Find nearest azimuth at this level
!
!-----------------------------------------------------------------------
!
azdiff=181.
jazmin=1
DO jazim=1,kntazim(kk)
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
daz=abs(daz)
IF(daz < azdiff) THEN
azdiff=daz
jazmin=jazim
END IF
END DO
jmirror=jazmin+(kntazim(kk)/2)
IF(jmirror > kntazim(kk)) jmirror=jmirror-kntazim(kk)
!
!-----------------------------------------------------------------------
!
! First pass, find median, avg, std dev.
!
!
! Loop forward from jazmin
!
!-----------------------------------------------------------------------
!
jend=kntazim(kk)
IF(jmirror > jazmin) jend=jmirror-1
jknt=0
DO jazim=jazmin,jend
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
IF(varvol(igate,jazim,kk) > varchek ) THEN
IF(kok < maxsort) THEN
IF(kok > 0 ) THEN
DO isort=1,kok
IF(varvol(igate,jazim,kk) < varsort(isort)) EXIT
END DO
DO jsort=kok,isort,-1
varsort(jsort+1)=varsort(jsort)
END DO
ELSE
isort=1
END IF
varsort(isort)=varvol(igate,jazim,kk)
END IF
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
END IF ! data ok
END IF ! inside box
END DO ! igate
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
!
!-----------------------------------------------------------------------
!
! IF kinbox > 0 continue from jazim=1
!
!-----------------------------------------------------------------------
!
IF((kinbox > 0 .OR. jknt <= jsrchmn) .and. jend==kntazim(kk)) THEN
DO jazim=1,jmirror-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
IF(varvol(igate,jazim,kk) > varchek ) THEN
IF(kok < maxsort) THEN
IF(kok > 0) THEN
DO isort=1,kok
IF(varvol(igate,jazim,kk) < varsort(isort)) EXIT
END DO
DO jsort=kok,isort,-1
varsort(jsort+1)=varsort(jsort)
END DO
ELSE
isort=1
END IF
varsort(isort)=varvol(igate,jazim,kk)
END IF
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
END IF
END IF
END DO
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Loop backward from jazmin
!
!-----------------------------------------------------------------------
!
jend= 1
IF(jmirror < jazmin) jend=jmirror
jknt=0
DO jazim=jazmin-1,jend,-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
IF(varvol(igate,jazim,kk) > varchek ) THEN
IF(kok < maxsort) THEN
IF(kok > 0) THEN
DO isort=1,kok
IF(varvol(igate,jazim,kk) < varsort(isort)) EXIT
END DO
DO jsort=kok,isort,-1
varsort(jsort+1)=varsort(jsort)
END DO
ELSE
isort=1
END IF
varsort(isort)=varvol(igate,jazim,kk)
END IF
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
END IF
END IF
END DO
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO
!
!-----------------------------------------------------------------------
!
! If not yet outside box, continue from last radial.
!
!-----------------------------------------------------------------------
!
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' OUT3 jazim,kinbox= ',jazim,kinbox
IF((kinbox > 0 .OR. jknt <= jsrchmn) .and. jend==1 ) THEN
DO jazim=kntazim(kk),jmirror,-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
IF(varvol(igate,jazim,kk) > varchek ) THEN
IF(kok < maxsort ) THEN
IF(kok > 0 ) THEN
DO isort=1,kok
IF(varvol(igate,jazim,kk) < varsort(isort)) EXIT
END DO
DO jsort=kok,isort,-1
varsort(jsort+1)=varsort(jsort)
END DO
ELSE
isort=1
END IF
varsort(isort)=varvol(igate,jazim,kk)
END IF
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
END IF
END IF
END DO ! igate
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
END IF
END DO ! kk
!
IF( kok == 0 ) CYCLE
varavg=sum/float(kok)
IF( kok <= maxsort ) THEN
mid=(kok/2)+1
varmed=varsort(mid)
ELSE
varmed=varavg
END IF
IF ( kok > 1 ) THEN
sdev=sqrt((sum2-(sum*sum/float(kok)))/float(kok-1))
thresh=max((2.*sdev),vmedlim)
ELSE
thresh=vmedlim
END IF
maxkok=max(maxkok,kok)
!
!-----------------------------------------------------------------------
!
! Process data for local quadratic fit
!
!-----------------------------------------------------------------------
!
DO kk=kend,kend
!
!-----------------------------------------------------------------------
!
! Find nearest azimuth at this level
!
!-----------------------------------------------------------------------
!
azdiff=181.
jazmin=1
DO jazim=1,kntazim(kk)
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
daz=abs(daz)
IF(daz < azdiff) THEN
azdiff=daz
jazmin=jazim
END IF
END DO
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' kk,azdiff,jazmin=', &
! kk,azdiff,jazmin,azmvol(jazmin,kk)
jmirror=jazmin+(kntazim(kk)/2)
IF(jmirror > kntazim(kk)) jmirror=jmirror-kntazim(kk)
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' jmirror=',jmirror
!
!-----------------------------------------------------------------------
!
! Loop forward from jazmin
!
!-----------------------------------------------------------------------
!
jend=kntazim(kk)
IF(jmirror > jazmin) jend=jmirror-1
jknt=0
DO jazim=jazmin,jend
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
ddxy=ddx*ddy
ddx2=ddx*ddx
ddy2=ddy*ddy
IF(varvol(igate,jazim,kk) > varchek .AND. &
abs(varvol(igate,jazim,kk)-varmed) < thresh ) THEN
!
varmax=max(varmax,varvol(igate,jazim,kk))
varmin=min(varmin,varvol(igate,jazim,kk))
!
rhsvar(1)=rhsvar(1)+varvol(igate,jazim,kk)
rhsvar(2)=rhsvar(2)+varvol(igate,jazim,kk)*ddx
rhsvar(3)=rhsvar(3)+varvol(igate,jazim,kk)*ddy
rhsvar(4)=rhsvar(4)+varvol(igate,jazim,kk)*ddxy
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddx2
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddy2
!
avar(1,1)=avar(1,1)+1.
avar(1,2)=avar(1,2)+ddx
avar(1,3)=avar(1,3)+ddy
avar(1,4)=avar(1,4)+ddxy
avar(1,5)=avar(1,5)+ddx2
avar(1,6)=avar(1,6)+ddy2
!
avar(2,1)=avar(2,1)+ddx
avar(2,2)=avar(2,2)+ddx2
avar(2,3)=avar(2,3)+ddx*ddy
avar(2,4)=avar(2,4)+ddx*ddxy
avar(2,5)=avar(2,5)+ddx*ddx2
avar(2,6)=avar(2,6)+ddx*ddy2
!
avar(3,1)=avar(3,1)+ddy
avar(3,2)=avar(3,2)+ddy*ddx
avar(3,3)=avar(3,3)+ddy2
avar(3,4)=avar(3,4)+ddy*ddx2
avar(3,5)=avar(3,5)+ddy*ddx2
avar(3,6)=avar(3,6)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddxy
avar(4,2)=avar(4,2)+ddxy*ddx
avar(4,3)=avar(4,3)+ddxy*ddy
avar(4,4)=avar(4,4)+ddxy*ddxy
avar(4,5)=avar(4,5)+ddxy*ddx2
avar(4,6)=avar(4,6)+ddxy*ddy2
!
avar(5,1)=avar(6,1)+ddx2
avar(5,2)=avar(6,2)+ddx2*ddx
avar(5,3)=avar(6,3)+ddx2*ddy
avar(5,4)=avar(6,4)+ddx2*ddxy
avar(5,5)=avar(6,5)+ddx2*ddx2
avar(5,6)=avar(6,6)+ddx2*ddy2
!
avar(6,1)=avar(6,1)+ddy2
avar(6,2)=avar(6,2)+ddy2*ddx
avar(6,3)=avar(6,3)+ddy2*ddy
avar(6,4)=avar(6,4)+ddy2*ddxy
avar(6,5)=avar(6,5)+ddy2*ddx2
avar(6,6)=avar(6,6)+ddy2*ddy2
!
END IF
!
END IF
END DO ! igate
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' jazim,azim,kinbox= ',jazim,azmvol(jazim,kk),kinbox
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
!
!-----------------------------------------------------------------------
!
! IF kinbox > 0 continue from jazim=1
!
!-----------------------------------------------------------------------
!
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' OUT1 jazim,azim,kinbox= ',&
! jazim,azmvol(jazim,kk),kinbox
IF((kinbox > 0 .OR. jknt <= jsrchmn) .and. jend==kntazim(kk)) THEN
DO jazim=1,jmirror-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
ddxy=ddx*ddy
ddx2=ddx*ddx
ddy2=ddy*ddy
IF(varvol(igate,jazim,kk) > varchek .AND. &
abs(varvol(igate,jazim,kk)-varmed) < thresh ) THEN
!
varmax=max(varmax,varvol(igate,jazim,kk))
varmin=min(varmin,varvol(igate,jazim,kk))
!
rhsvar(1)=rhsvar(1)+varvol(igate,jazim,kk)
rhsvar(2)=rhsvar(2)+varvol(igate,jazim,kk)*ddx
rhsvar(3)=rhsvar(3)+varvol(igate,jazim,kk)*ddy
rhsvar(4)=rhsvar(4)+varvol(igate,jazim,kk)*ddxy
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddx2
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddy2
!
avar(1,1)=avar(1,1)+1.
avar(1,2)=avar(1,2)+ddx
avar(1,3)=avar(1,3)+ddy
avar(1,4)=avar(1,4)+ddxy
avar(1,5)=avar(1,5)+ddx2
avar(1,6)=avar(1,6)+ddy2
!
avar(2,1)=avar(2,1)+ddx
avar(2,2)=avar(2,2)+ddx2
avar(2,3)=avar(2,3)+ddx*ddy
avar(2,4)=avar(2,4)+ddx*ddxy
avar(2,5)=avar(2,5)+ddx*ddx2
avar(2,6)=avar(2,6)+ddx*ddy2
!
avar(3,1)=avar(3,1)+ddy
avar(3,2)=avar(3,2)+ddy*ddx
avar(3,3)=avar(3,3)+ddy2
avar(3,4)=avar(3,4)+ddy*ddxy
avar(3,5)=avar(3,5)+ddy*ddx2
avar(3,6)=avar(3,6)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddxy
avar(4,2)=avar(4,2)+ddxy*ddx
avar(4,3)=avar(4,3)+ddxy*ddy
avar(4,4)=avar(4,4)+ddxy*ddxy
avar(4,5)=avar(4,5)+ddxy*ddx2
avar(4,6)=avar(4,6)+ddxy*ddy2
avar(5,1)=avar(5,1)+ddx2
avar(5,2)=avar(5,2)+ddx2*ddx
avar(5,3)=avar(5,3)+ddx2*ddy
avar(5,4)=avar(5,4)+ddx2*ddxy
avar(5,5)=avar(5,5)+ddx2*ddx2
avar(5,6)=avar(5,6)+ddx2*ddy2
!
avar(6,1)=avar(6,1)+ddy2
avar(6,2)=avar(6,2)+ddy2*ddx
avar(6,3)=avar(6,3)+ddy2*ddy
avar(6,4)=avar(6,4)+ddy2*ddxy
avar(6,5)=avar(6,5)+ddy2*ddx2
avar(6,6)=avar(6,6)+ddy2*ddy2
!
END IF
!
END IF
END DO ! igate
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' jazim,azim,kinbox= ',jazim,azmvol(jazim,kk),kinbox
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
END IF
!
!-----------------------------------------------------------------------
!
! Loop backward from jazmin
!
!-----------------------------------------------------------------------
!
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' OUT2 jazim,kinbox= ',jazim,kinbox
jend= 1
IF(jmirror < jazmin) jend=jmirror
jknt=0
DO jazim=jazmin-1,jend,-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
ddxy=ddx*ddy
ddx2=ddx*ddx
ddy2=ddy*ddy
IF(varvol(igate,jazim,kk) > varchek .AND. &
abs(varvol(igate,jazim,kk)-varmed) < thresh ) THEN
!
varmax=max(varmax,varvol(igate,jazim,kk))
varmin=min(varmin,varvol(igate,jazim,kk))
!
rhsvar(1)=rhsvar(1)+varvol(igate,jazim,kk)
rhsvar(2)=rhsvar(2)+varvol(igate,jazim,kk)*ddx
rhsvar(3)=rhsvar(3)+varvol(igate,jazim,kk)*ddy
rhsvar(4)=rhsvar(4)+varvol(igate,jazim,kk)*ddxy
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddx2
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddy2
!
avar(1,1)=avar(1,1)+1.
avar(1,2)=avar(1,2)+ddx
avar(1,3)=avar(1,3)+ddy
avar(1,4)=avar(1,4)+ddxy
avar(1,5)=avar(1,5)+ddx2
avar(1,6)=avar(1,6)+ddy2
!
avar(2,1)=avar(2,1)+ddx
avar(2,2)=avar(2,2)+ddx2
avar(2,3)=avar(2,3)+ddx*ddy
avar(2,4)=avar(2,4)+ddx*ddxy
avar(2,5)=avar(2,5)+ddx*ddx2
avar(2,6)=avar(2,6)+ddx*ddy2
!
avar(3,1)=avar(3,1)+ddy
avar(3,2)=avar(3,2)+ddy*ddx
avar(3,3)=avar(3,3)+ddy2
avar(3,4)=avar(3,4)+ddy*ddxy
avar(3,5)=avar(3,5)+ddy*ddx2
avar(3,6)=avar(3,6)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddxy
avar(4,2)=avar(4,2)+ddxy*ddx
avar(4,3)=avar(4,3)+ddxy*ddy
avar(4,4)=avar(4,4)+ddxy*ddxy
avar(4,5)=avar(4,5)+ddxy*ddx2
avar(4,6)=avar(4,6)+ddxy*ddy2
!
avar(5,1)=avar(5,1)+ddx2
avar(5,2)=avar(5,2)+ddx2*ddx
avar(5,3)=avar(5,3)+ddx2*ddy
avar(5,4)=avar(5,4)+ddx2*ddxy
avar(5,5)=avar(5,5)+ddx2*ddx2
avar(5,6)=avar(5,6)+ddx2*ddy2
!
avar(6,1)=avar(6,1)+ddy2
avar(6,2)=avar(6,2)+ddy2*ddx
avar(6,3)=avar(6,3)+ddy2*ddy
avar(6,4)=avar(6,4)+ddy2*ddxy
avar(6,5)=avar(6,5)+ddy2*ddx2
avar(6,6)=avar(6,6)+ddy2*ddy2
!
END IF
!
END IF
END DO ! igate
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' jazim,azim,kinbox= ',jazim,azmvol(jazim,kk),kinbox
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
!
!-----------------------------------------------------------------------
!
! If not yet outside box, continue from last radial.
!
!-----------------------------------------------------------------------
!
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' OUT3 jazim,kinbox= ',jazim,kinbox
IF((kinbox > 0 .OR. jknt <= jsrchmn) .and. jend==1 ) THEN
DO jazim=kntazim(kk),jmirror,-1
kinbox=0
jknt=jknt+1
daz=azmvol(jazim,kk)-azimijk
IF(daz > 180.) daz=daz-360.
IF(daz < -180.) daz=daz+360.
IF(abs(daz) > dazlim) EXIT
DO igate=1,kntgate(jazim,kk)
!
ddx=rxvol(igate,jazim,kk)-xs
ddy=ryvol(igate,jazim,kk)-ys
!
IF( rngvol(igate,kk) > rngmin .AND. &
rngvol(igate,kk) < rngmax .AND. &
abs(ddx) < dxthr .AND. abs(ddy) < dxthr ) THEN
kinbox=kinbox+1
ddxy=ddx*ddy
ddx2=ddx*ddx
ddy2=ddy*ddy
IF(varvol(igate,jazim,kk) > varchek .AND. &
abs(varvol(igate,jazim,kk)-varmed) < thresh ) THEN
!
varmax=max(varmax,varvol(igate,jazim,kk))
varmin=min(varmin,varvol(igate,jazim,kk))
!
rhsvar(1)=rhsvar(1)+varvol(igate,jazim,kk)
rhsvar(2)=rhsvar(2)+varvol(igate,jazim,kk)*ddx
rhsvar(3)=rhsvar(3)+varvol(igate,jazim,kk)*ddy
rhsvar(4)=rhsvar(4)+varvol(igate,jazim,kk)*ddxy
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddx2
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddy2
!
avar(1,1)=avar(1,1)+1.
avar(1,2)=avar(1,2)+ddx
avar(1,3)=avar(1,3)+ddy
avar(1,4)=avar(1,4)+ddxy
avar(1,5)=avar(1,5)+ddx2
avar(1,6)=avar(1,6)+ddy2
!
avar(2,1)=avar(2,1)+ddx
avar(2,2)=avar(2,2)+ddx2
avar(2,3)=avar(2,3)+ddx*ddy
avar(2,4)=avar(2,4)+ddx*ddxy
avar(2,5)=avar(2,5)+ddx*ddx2
avar(2,6)=avar(2,6)+ddx*ddy2
!
avar(3,1)=avar(3,1)+ddy
avar(3,2)=avar(3,2)+ddy*ddx
avar(3,3)=avar(3,3)+ddy2
avar(3,4)=avar(3,4)+ddy*ddxy
avar(3,5)=avar(3,5)+ddy*ddx2
avar(3,6)=avar(3,6)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddx2
avar(4,2)=avar(4,2)+ddx2*ddx
avar(4,3)=avar(4,3)+ddx2*ddy
avar(4,4)=avar(4,4)+ddx2*ddxy
avar(4,5)=avar(4,5)+ddx2*ddx2
avar(4,6)=avar(4,6)+ddx2*ddy2
!
avar(5,1)=avar(5,1)+ddx2
avar(5,2)=avar(5,2)+ddx2*ddx
avar(5,3)=avar(5,3)+ddx2*ddy
avar(5,4)=avar(5,4)+ddx2*ddxy
avar(5,5)=avar(5,5)+ddx2*ddx2
avar(5,6)=avar(5,6)+ddx2*ddy2
!
avar(6,1)=avar(6,1)+ddy2
avar(6,2)=avar(6,2)+ddy2*ddx
avar(6,3)=avar(6,3)+ddy2*ddy
avar(6,4)=avar(6,4)+ddy2*ddxy
avar(6,5)=avar(6,5)+ddy2*ddx2
avar(6,6)=avar(6,6)+ddy2*ddy2
!
END IF
END IF
END DO ! igate
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' jazim,azim,kinbox=',jazim,azmvol(jazim,kk),kinbox
IF(kinbox == 0 .AND. jknt > jsrchmn) EXIT
END DO ! jazim
END IF
!
END DO ! kk
! IF(i==10 .and. j==10 .and. k==5) &
! print *, ' FINAL OUT jjazim,kinbox= ',jazim,kinbox
!
!-----------------------------------------------------------------------
!
! Solve for variable at grid point
!
!-----------------------------------------------------------------------
!
knt=nint(avar(1,1))
IF ( iorder > 1 .and. knt > 6 ) THEN
varmean=rhsvar(1)/avar(1,1)
CALL GJELIM(n,avar,rhsvar,sol,work,work1d,eps,istatus)
fitvalue=min(varmax,max(varmin,sol(1)))
! write(6,'(3i3,i7,4f7.1)') i,j,k,knt, &
! varmin,varmean,varmax,fitvalue
ELSE IF ( iorder > 0 .and. knt > 5 ) THEN
DO jj=1,4
DO ii=1,4
array(ii,jj)=avar(ii,jj)
END DO
END DO
DO ii=1,4
rhsv(ii)=rhsvar(ii)
END DO
CALL GJELIM(4,array,rhsv,solv,work,work1d,eps,istatus)
fitvalue=min(varmax,max(varmin,solv(1)))
! write(6,'(3i3,a,i6,6f7.1)') i,j,k,'*',knt, &
! varmin,varmean,varmax,gridvar(i,j,k)
ELSE IF ( knt > 0 ) THEN
varmean=rhsvar(1)/avar(1,1)
fitvalue=varmean
! write(6,'(3i3,a,i6,6f7.1)') i,j,k,'*',knt, &
! varmin,varmean,varmax,fitvalue
END IF
varvol(iigate,jjazim,kkelev) = fitvalue
END IF
END IF
ENDIF
END DO
END DO
END DO
IF(maxkok > maxsort) THEN
WRITE(6,'(a,i10,/a,i10,a,/a)') &
' quadfit: Note ... max gate count =',maxkok, &
' exceeded sort space allocated =',maxsort, &
' for median check.', &
' Where maxsort is exceeded avg used in place of median for QC.'
ELSE
WRITE(6,'(a,i10,/a,i10)') &
' quadfit: FYI max gate count in grid vol =',maxkok, &
' Sort space allocated for median check =',maxsort
END IF
RETURN
END SUBROUTINE quadfit
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRTVEL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wrtvel(i_angle,i_tilt,varid, &,2
iyr,imon,iday,ihr,imin,isec, &
gsp_vel,rfrst_vel, &
maxgate,maxazim,ngate,nazim, &
azim,elev,rvel);
!-----------------------------------------------------------------------
!
! PURPOSE:
! Write out a tilt of radar data in radar coordinates for plotting
!
!-----------------------------------------------------------------------
!
! 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
! i_angle elevation angle
! i_tilt tilt number
! varid variable id
! iyr,imon,iday,ihr,imin,isec Data time
! gsp_vel gate spacing
! rfrst_vel distance of first gate to radar
! azim azimuthal angles for each radial
! elev elevations for each radial
! rvel Doppler radial velocity
!
! OUTPUT:
! None
!
!
! WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: maxgate
INTEGER :: maxazim
INTEGER :: ngate
INTEGER :: nazim
INTEGER :: i, j
INTEGER :: i_angle,i_tilt
INTEGER :: iyr,imon,iday,ihr,imin,isec
INTEGER :: gsp_vel,rfrst_vel
REAL :: rvel(maxgate,maxazim)
REAL :: elev(maxazim)
REAL :: azim(maxazim)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=6) :: varid
CHARACTER (LEN=256) :: tmpdirnam
CHARACTER (LEN=256) :: vfnam
INTEGER :: nunit
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL gtlfnkey( runname, lfnkey )
tmpdirnam = dirname
IF ( LEN(tmpdirnam) == 0 .OR. tmpdirnam == ' ' ) THEN
tmpdirnam = '.'
END IF
WRITE (vfnam,'(4A,A6,i6.6)') tmpdirnam(:INDEX(tmpdirnam,' ')-1), &
'/', runname(1:lfnkey), '.', varid, i_angle
CALL getunit( nunit)
write(*,'(a)') vfnam
open(nunit,file=TRIM(vfnam),form='formatted',status='unknown')
write(nunit,1003)iyr,imon,iday,ihr,imin,isec
write(nunit,1004) gsp_vel
write(nunit,1004) rfrst_vel
write(nunit,1004) i_angle
write(nunit,1004) i_tilt
write(nunit,1004) ngate
write(nunit,1004) nazim
write(nunit,1005) (elev(i),i=1,nazim)
write(nunit,1005) (azim(i),i=1,nazim)
DO 315 j=1,nazim
write(nunit,1005) (rvel(i,j),i=1,ngate)
315 CONTINUE
1003 format(6i4)
1004 format(i8)
1005 format(16f8.2)
!
close(nunit)
CALL retunit(nunit)
END SUBROUTINE wrtvel