!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RADCOORD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE radcoord(nx,ny,nz, & 1,4
x,y,z,zp,xs,ys,zps, &
radar_lat,radar_lon,radarx,radary)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Intialize coordinate fields for the radar remapping routines.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! radar_lat Latitude (degrees) of radar
! radar_lon Longitude (positive East) of radar
!
! OUTPUT:
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space(m)
!
! 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)
!
! radarx x location of radar in grid
! radary y location of radar in grid
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx
INTEGER, INTENT(IN) :: ny
INTEGER, INTENT(IN) :: nz
REAL, INTENT(IN) :: x(nx)
REAL, INTENT(IN) :: y(ny)
REAL, INTENT(IN) :: z(nz)
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) :: radar_lat
REAL, INTENT(IN) :: radar_lon
REAL, INTENT(OUT) :: radarx
REAL, INTENT(OUT) :: radary
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER i,j,k
REAL ctrx,ctry,swx,swy
REAL latnot(2)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Set the coordinates of the scalar grid points
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
xs(nx)=2.0*xs(nx-1)-xs(nx-2)
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
ys(ny)=2.0*ys(ny-1)-ys(ny-2)
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*zps(i,j,nz-1)-zps(i,j,nz-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set up map projection
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
print *, ' latnot,trulon=',latnot(1),latnot(2),trulon
CALL setmapr
(mapproj,sclfct,latnot,trulon)
!
!-----------------------------------------------------------------------
!
! Get lat, lon of grid center and lat, lon of radar
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,ctrlat,ctrlon,ctrx,ctry)
print *, ' ctrlat,ctrlon,ctrx,ctry=',ctrlat,ctrlon,ctrx,ctry
swx=ctrx-0.5*(nx-3)*dx
swy=ctry-0.5*(ny-3)*dy
CALL setorig( 1, swx, swy)
CALL lltoxy(1,1,radar_lat,radar_lon,radarx,radary)
print *, ' radar_lat,radar_lon: ',radar_lat,radar_lon
print *, ' grid radarx = ',radarx,' radary = ',radary
RETURN
END SUBROUTINE radcoord
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RMPINIT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rmpinit(nx,ny,nz,mxrefgat,mxvelgat,maxazim,maxelev, & 1
kntrgat,kntrazm,kntrelv, &
kntvgat,kntvazm,kntvelv, &
nyqvvol,timevol, &
rngrvol,azmrvol,elvrvol,refvol, &
rngvvol,azmvvol,elvvvol,velvol, &
gridvel,gridref,gridnyq,gridtim)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Intialize fields for the radar remapping routines.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
!
! MODIFICATION HISTORY:
!
! 05/17/2002 (Keith Brewster)
! Added initialization of gridvel,gridref,gridnyq,gridtim
! to missing value.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! mxrefgat Maximum number of reflectivity gates
! mxvelgat Maximum number of velocity gates
! maxazim Maximum number of azimuth angles per tilt
! maxelev Maximum number of elevation angles (tilts)
!
! OUTPUT
!
! kntrgat Number of reflectivity gates
! kntrazm Number of reflectivity radials
! kntvelv Number of reflectivity tilts
!
! kntvgat Number of velocity gates
! kntvazm Number of velocity radials
! kntvelv Number of velocity tilts
!
! rngrvol Range to gate in reflectivity 3-D volume
! azmrvol Azimuth angle in reflectivity 3-D volume
! elvrvol Elevation angle in reflectivity 3-D volume
! refvol Reflectivity 3-D volume
!
! 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
! velvol Velocity 3-D volume
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nx
INTEGER, INTENT(IN) :: ny
INTEGER, INTENT(IN) :: nz
INTEGER, INTENT(IN) :: mxrefgat
INTEGER, INTENT(IN) :: mxvelgat
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(OUT) :: kntrgat(maxazim,maxelev)
INTEGER, INTENT(OUT) :: kntrazm(maxelev)
INTEGER, INTENT(OUT) :: kntrelv
INTEGER, INTENT(OUT) :: kntvgat(maxazim,maxelev)
INTEGER, INTENT(OUT) :: kntvazm(maxelev)
INTEGER, INTENT(OUT) :: kntvelv
REAL, INTENT(OUT) :: nyqvvol(maxelev)
INTEGER, INTENT(OUT) :: timevol(maxazim,maxelev)
REAL, INTENT(OUT) :: rngrvol(mxrefgat,maxelev)
REAL, INTENT(OUT) :: azmrvol(maxazim,maxelev)
REAL, INTENT(OUT) :: elvrvol(maxazim,maxelev)
REAL, INTENT(OUT) :: refvol(mxrefgat,maxazim,maxelev)
REAL, INTENT(OUT) :: rngvvol(mxvelgat,maxelev)
REAL, INTENT(OUT) :: azmvvol(maxazim,maxelev)
REAL, INTENT(OUT) :: elvvvol(maxazim,maxelev)
REAL, INTENT(OUT) :: velvol(mxvelgat,maxazim,maxelev)
REAL, INTENT(OUT) :: gridvel(nx,ny,nz)
REAL, INTENT(OUT) :: gridref(nx,ny,nz)
REAL, INTENT(OUT) :: gridnyq(nx,ny,nz)
REAL, INTENT(OUT) :: gridtim(nx,ny,nz)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Initialize counters to zero and fill data arrays with missing value
!
!-----------------------------------------------------------------------
!
kntrelv = 0
kntrazm = 0
kntrgat = 0
kntvelv = 0
kntvazm = 0
kntvgat = 0
nyqvvol = 0.
timevol = 0
elvrvol = -999.
azmrvol = -999.
rngrvol = -999.
refvol = -999.
elvvvol = -999.
azmvvol = -999.
rngvvol = -999.
velvol = -999.
gridvel = -999.
gridref = -999.
gridnyq = -999.
gridtim = -999.
RETURN
END SUBROUTINE rmpinit
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE VOLBUILD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE volbuild(maxgate,maxazim,maxelev,ngate,nazim, & 2
nyqset,timeset, &
kntgate,kntazim,kntelev, &
gatesp,rfirstg,varchek, &
vnyquist,time, &
azim,elev,vartilt, &
nyqvvol,timevol, &
rngvol,azmvol,elvvol,varvol)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Build a 3-D volume of radar data from 2-D tilts.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
!
! ngate Number of gates in each radial
! nazim Number of radials in each tilt
!
!
! gatesp Gate spacing
! rfirstg Range (m) to first gate
! varchek Threshold value to determine good vs. flagged data
! azim Azimuth angles
! elev Elevation angles
! vartilt Radar variable (reflectivity, velocity or spectrum width)
!
! OUTPUT:
!
! kntgate Number of gates in each radial (3-D)
! kntazim Number of radials in each tilt (3-D)
! kntelev Number of elevation angles (tilts)
! rngvol Range to gate
! azmvol Azimuth angle
! elvvol Elevation angle
! varvol Radar variables accumulated in 3-D array
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
INTEGER, INTENT(IN) :: nyqset
INTEGER, INTENT(IN) :: timeset
INTEGER, INTENT(INOUT) :: kntgate(maxazim,maxelev)
INTEGER, INTENT(INOUT) :: kntazim(maxelev)
INTEGER, INTENT(INOUT) :: kntelev
INTEGER, INTENT(IN) :: gatesp
INTEGER, INTENT(IN) :: rfirstg
REAL, INTENT(IN) :: varchek
REAL, INTENT(IN) :: vnyquist
REAL, INTENT(IN) :: time(maxazim)
REAL, INTENT(IN) :: azim(maxazim)
REAL, INTENT(IN) :: elev(maxazim)
REAL, INTENT(IN) :: vartilt(maxgate,maxazim)
REAL, INTENT(OUT) :: nyqvvol(maxelev)
REAL, INTENT(OUT) :: timevol(maxazim,maxelev)
REAL, INTENT(OUT) :: rngvol(maxgate,maxelev)
REAL, INTENT(OUT) :: azmvol(maxazim,maxelev)
REAL, INTENT(OUT) :: elvvol(maxazim,maxelev)
REAL, INTENT(OUT) :: varvol(maxgate,maxazim,maxelev)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: igate,jazim,kelev
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
kntelev = kntelev + 1
kelev = kntelev
kntazim(kelev) = nazim
IF(nyqset > 0 ) nyqvvol(kelev) = vnyquist
DO igate = 1, maxgate
rngvol(igate,kelev)=rfirstg+(igate-1)*gatesp
END DO
IF( timeset > 0 ) THEN
DO jazim = 1, nazim
timevol(jazim,kelev)=time(jazim)
azmvol(jazim,kelev)=azim(jazim)
elvvol(jazim,kelev)=elev(jazim)
END DO
ELSE
DO jazim = 1, nazim
azmvol(jazim,kelev)=azim(jazim)
elvvol(jazim,kelev)=elev(jazim)
END DO
END IF
DO jazim = 1, nazim
DO igate=1, ngate
varvol(igate,jazim,kelev)=vartilt(igate,jazim)
IF (vartilt(igate,jazim) > varchek) &
kntgate(jazim,kelev)=max(kntgate(jazim,kelev),igate)
END DO
END DO
RETURN
END SUBROUTINE volbuild
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE REMAPVOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE remapvol(maxgate,maxazim,maxelev,nx,ny,nz, & 2,4
vardump,nyqset,timeset, &
varid,varname,varunits, &
varchek,varmiss,vmedlim,dazlim,iorder, &
kntgate,kntazim,kntelev, &
radarx,radary,rdralt,dazim, &
rngmin,rngmax,time1st, &
rngvol,azmvol,elvvol, &
varvol,timevol,nyqvvol,rxvol,ryvol,rzvol, &
xs,ys,zps, &
gridvar,gridtim,gridnyq,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: Keith Brewster, CAPS
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! varid Radar variable ID for diagnostic file writing
! varname Name of radar variable for diagnostic file writing
! varunit Units of radar variable for diagnostic file writing
!
! varchek Threshold for checking data, good vs. flagged
! varmiss Value to assign to data for missing
! 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
! 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:
!
! rxvol x-coordinate at radar data location
! ryvol y-coordinate at radar data location
! rzvol z-coordinate at radar data location
!
! gridvar Radar variable remapped to scalar points
!
! istatus Status indicator
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: nx
INTEGER, INTENT(IN) :: ny
INTEGER, INTENT(IN) :: nz
INTEGER, INTENT(IN) :: vardump
INTEGER, INTENT(IN) :: nyqset
INTEGER, INTENT(IN) :: timeset
CHARACTER (LEN=6), INTENT(IN) :: varid
CHARACTER (LEN=20), INTENT(IN) :: varname
CHARACTER (LEN=20), INTENT(IN) :: varunits
REAL, INTENT(IN) :: varchek
REAL, INTENT(IN) :: varmiss
REAL, INTENT(IN) :: vmedlim
REAL, INTENT(IN) :: dazlim
INTEGER, INTENT(IN) :: iorder
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
INTEGER, INTENT(IN) :: time1st
REAL, INTENT(IN) :: nyqvvol(maxelev)
INTEGER, INTENT(IN) :: timevol(maxazim,maxelev)
REAL, INTENT(IN) :: rngvol(maxgate,maxelev)
REAL, INTENT(IN) :: azmvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvvol(maxazim,maxelev)
REAL, INTENT(IN) :: varvol(maxgate,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(IN) :: xs(nx)
REAL, INTENT(IN) :: ys(ny)
REAL, INTENT(IN) :: zps(nx,ny,nz)
REAL, INTENT(OUT) :: gridvar(nx,ny,nz)
REAL, INTENT(OUT) :: gridnyq(nx,ny,nz)
REAL, INTENT(OUT) :: gridtim(nx,ny,nz)
INTEGER, INTENT(OUT) :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. Local Variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: n = 7
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)
REAL, PARAMETER :: eps = 1.0E-25
INTEGER :: ii,jj,kk,i,j,k,knt,kinbox
INTEGER :: kok,isort,jsort,mid
INTEGER :: kbgn,kend
INTEGER :: igate,jazim,kelev,jazmin,jmirror,jend
INTEGER :: istatal,istatwrt
INTEGER, PARAMETER :: maxsort = 5000
REAL :: deg2rad,rad2deg
REAL :: elevmin,elevmax
REAL :: delx,dely,delz,dazimr,daz,azdiff
REAL :: ddx,ddxy,ddx2,ddy,ddy2,ddz,dxthr,dxthr0
REAL :: cosaz,sinaz,sfcr,zagl
REAL :: sum,sum2,sdev,thresh,slrange,elijk,azimijk,time
REAL :: varmax,varmin,varavg,varmean,varmed
REAL, ALLOCATABLE :: varsort(:)
REAL, ALLOCATABLE :: elevmean(:)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
deg2rad = atan(1.)/45.
rad2deg = 1./deg2rad
dazimr = dazim*deg2rad
dxthr0=0.6*max((xs(3)-xs(2)),(ys(3)-ys(2)))
allocate(elevmean(kntelev),stat=istatal)
allocate(varsort(maxsort),stat=istatal)
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
elevmean(kelev)=sum/float(kntazim(kelev))
END DO
!
write(6,'(a)') &
' i j k knt min mean max intrp'
DO k=1,nz
DO j=1,ny
DO i=1,nx
gridvar(i,j,k)=varmiss
gridtim(i,j,k)=0.
gridnyq(i,j,k)=0.
END DO
END DO
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 BEAMHGT(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
print *, ' REMAPVOL: Finished loop 3 '
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
kok=0
sum=0.
sum2=0.
sdev=0.
varavg=999999.
varsort=999999.
delx=xs(i)-radarx
dely=ys(j)-radary
delz=zps(i,j,k)-rdralt
slrange=sqrt(delx*delx+dely*dely+delz*delz)
dxthr=max(dxthr0,((slrange+dxthr0)*dazimr))
!
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(elevmean(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=kbgn,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)
IF(nyqset > 0) gridnyq(i,j,k) = nyqvvol(kk)
IF(timeset > 0) gridtim(i,j,k) = &
float(timevol(jazim,kk)-time1st)
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
!
!-----------------------------------------------------------------------
!
! First pass, find median, avg, std dev.
!
!
! Loop forward from jazmin
!
!-----------------------------------------------------------------------
!
jend=kntazim(kk)
IF(jmirror > jazmin) jend=jmirror-1
DO jazim=jazmin,jend
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
!
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
DO isort=1,kok
IF(varvol(igate,jazim,kk) < varsort(isort)) EXIT
END DO
IF(kok < maxsort) THEN
DO jsort=kok,isort,-1
varsort(jsort+1)=varsort(jsort)
END DO
varsort(isort)=varvol(igate,jazim,kk)
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
ELSE
EXIT
END IF
END IF ! data ok
END IF ! inside box
END DO ! igate
IF(kinbox == 0 .OR. kok == maxsort) 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 .and. jend==kntazim(kk)) THEN
DO jazim=1,jmirror-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
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
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
varsort(isort)=varvol(igate,jazim,kk)
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
ELSE
EXIT
END IF
END IF
END IF
END DO
IF(kinbox == 0 .OR. kok == maxsort) EXIT
END DO
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
DO jazim=jazmin-1,jend,-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
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
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
varsort(isort)=varvol(igate,jazim,kk)
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
ELSE
EXIT
END IF
END IF
END IF
END DO
IF(kinbox == 0 .OR. kok == maxsort) 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 .and. jend==1 ) THEN
DO jazim=kntazim(kk),jmirror,-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
!
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
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
varsort(isort)=varvol(igate,jazim,kk)
sum=sum+varvol(igate,jazim,kk)
sum2=sum2+(varvol(igate,jazim,kk)*varvol(igate,jazim,kk))
kok=kok+1
ELSE
EXIT
END IF
END IF
END IF
END DO ! igate
IF(kinbox == 0 .OR. kok == maxsort) EXIT
END DO ! jazim
END IF
END DO ! kk
!
mid=(kok/2)+1
varmed=varsort(mid)
IF(kok > 0) varavg=sum/float(kok)
IF ( kok > 1 ) &
sdev=sqrt((sum2-(sum*sum/float(kok)))/float(kok-1))
thresh=max((2.*sdev),vmedlim)
!
!-----------------------------------------------------------------------
!
! Process data for local quadratic fit
!
!-----------------------------------------------------------------------
!
DO kk=kbgn,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
DO jazim=jazmin,jend
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
ddz=rzvol(igate,jazim,kk)-zps(i,j,k)
!
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)*ddz
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddxy
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddx2
rhsvar(7)=rhsvar(7)+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)+ddz
avar(1,5)=avar(1,5)+ddxy
avar(1,6)=avar(1,6)+ddx2
avar(1,7)=avar(1,7)+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*ddz
avar(2,5)=avar(2,5)+ddx*ddxy
avar(2,6)=avar(2,6)+ddx*ddx2
avar(2,7)=avar(2,7)+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*ddz
avar(3,5)=avar(3,5)+ddy*ddx2
avar(3,6)=avar(3,6)+ddy*ddx2
avar(3,7)=avar(3,7)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddz
avar(4,2)=avar(4,2)+ddz*ddx
avar(4,3)=avar(4,3)+ddz*ddy
avar(4,4)=avar(4,4)+ddz*ddz
avar(4,5)=avar(4,5)+ddz*ddxy
avar(4,6)=avar(4,6)+ddz*ddx2
avar(4,7)=avar(4,7)+ddz*ddy2
!
avar(5,1)=avar(5,1)+ddxy
avar(5,2)=avar(5,2)+ddxy*ddx
avar(5,3)=avar(5,3)+ddxy*ddy
avar(5,4)=avar(5,4)+ddxy*ddz
avar(5,5)=avar(5,5)+ddxy*ddxy
avar(5,6)=avar(5,6)+ddxy*ddx2
avar(5,7)=avar(5,7)+ddxy*ddy2
!
avar(6,1)=avar(6,1)+ddx2
avar(6,2)=avar(6,2)+ddx2*ddx
avar(6,3)=avar(6,3)+ddx2*ddy
avar(6,4)=avar(6,4)+ddx2*ddz
avar(6,5)=avar(6,5)+ddx2*ddxy
avar(6,6)=avar(6,6)+ddx2*ddx2
avar(6,7)=avar(6,7)+ddx2*ddy2
!
avar(7,1)=avar(7,1)+ddy2
avar(7,2)=avar(7,2)+ddy2*ddx
avar(7,3)=avar(7,3)+ddy2*ddy
avar(7,4)=avar(7,4)+ddy2*ddz
avar(7,5)=avar(7,5)+ddy2*ddxy
avar(7,6)=avar(7,6)+ddy2*ddx2
avar(7,7)=avar(7,7)+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) 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 .and. jend==kntazim(kk)) THEN
DO jazim=1,jmirror-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
ddz=rzvol(igate,jazim,kk)-zps(i,j,k)
!
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)*ddz
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddxy
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddx2
rhsvar(7)=rhsvar(7)+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)+ddz
avar(1,5)=avar(1,5)+ddxy
avar(1,6)=avar(1,6)+ddx2
avar(1,7)=avar(1,7)+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*ddz
avar(2,5)=avar(2,5)+ddx*ddxy
avar(2,6)=avar(2,6)+ddx*ddx2
avar(2,7)=avar(2,7)+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*ddz
avar(3,5)=avar(3,5)+ddy*ddxy
avar(3,6)=avar(3,6)+ddy*ddx2
avar(3,7)=avar(3,7)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddz
avar(4,2)=avar(4,2)+ddz*ddx
avar(4,3)=avar(4,3)+ddz*ddy
avar(4,4)=avar(4,4)+ddz*ddz
avar(4,5)=avar(4,5)+ddz*ddxy
avar(4,6)=avar(4,6)+ddz*ddx2
avar(4,7)=avar(4,7)+ddz*ddy2
!
avar(5,1)=avar(5,1)+ddxy
avar(5,2)=avar(5,2)+ddxy*ddx
avar(5,3)=avar(5,3)+ddxy*ddy
avar(5,4)=avar(5,4)+ddxy*ddz
avar(5,5)=avar(5,5)+ddxy*ddxy
avar(5,6)=avar(5,6)+ddxy*ddx2
avar(5,7)=avar(5,7)+ddxy*ddy2
avar(6,1)=avar(6,1)+ddx2
avar(6,2)=avar(6,2)+ddx2*ddx
avar(6,3)=avar(6,3)+ddx2*ddy
avar(6,4)=avar(6,4)+ddx2*ddz
avar(6,5)=avar(6,5)+ddx2*ddxy
avar(6,6)=avar(6,6)+ddx2*ddx2
avar(6,7)=avar(6,7)+ddx2*ddy2
!
avar(7,1)=avar(7,1)+ddy2
avar(7,2)=avar(7,2)+ddy2*ddx
avar(7,3)=avar(7,3)+ddy2*ddy
avar(7,4)=avar(7,4)+ddy2*ddz
avar(7,5)=avar(7,5)+ddy2*ddxy
avar(7,6)=avar(7,6)+ddy2*ddx2
avar(7,7)=avar(7,7)+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) 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
DO jazim=jazmin-1,jend,-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
ddz=rzvol(igate,jazim,kk)-zps(i,j,k)
!
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)*ddz
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddxy
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddx2
rhsvar(7)=rhsvar(7)+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)+ddz
avar(1,5)=avar(1,5)+ddxy
avar(1,6)=avar(1,6)+ddx2
avar(1,7)=avar(1,7)+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*ddz
avar(2,5)=avar(2,5)+ddx*ddxy
avar(2,6)=avar(2,6)+ddx*ddx2
avar(2,7)=avar(2,7)+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*ddz
avar(3,5)=avar(3,5)+ddy*ddxy
avar(3,6)=avar(3,6)+ddy*ddx2
avar(3,7)=avar(3,7)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddz
avar(4,2)=avar(4,2)+ddz*ddx
avar(4,3)=avar(4,3)+ddz*ddy
avar(4,4)=avar(4,4)+ddz*ddz
avar(4,5)=avar(4,5)+ddz*ddxy
avar(4,6)=avar(4,6)+ddz*ddx2
avar(4,7)=avar(4,7)+ddz*ddy2
!
avar(5,1)=avar(5,1)+ddxy
avar(5,2)=avar(5,2)+ddxy*ddx
avar(5,3)=avar(5,3)+ddxy*ddy
avar(5,4)=avar(5,4)+ddxy*ddz
avar(5,5)=avar(5,5)+ddxy*ddxy
avar(5,6)=avar(5,6)+ddxy*ddx2
avar(5,7)=avar(5,7)+ddxy*ddy2
!
avar(6,1)=avar(6,1)+ddx2
avar(6,2)=avar(6,2)+ddx2*ddx
avar(6,3)=avar(6,3)+ddx2*ddy
avar(6,4)=avar(6,4)+ddx2*ddz
avar(6,5)=avar(6,5)+ddx2*ddxy
avar(6,6)=avar(6,6)+ddx2*ddx2
avar(6,7)=avar(6,7)+ddx2*ddy2
!
avar(7,1)=avar(7,1)+ddy2
avar(7,2)=avar(7,2)+ddy2*ddx
avar(7,3)=avar(7,3)+ddy2*ddy
avar(7,4)=avar(7,4)+ddy2*ddz
avar(7,5)=avar(7,5)+ddy2*ddxy
avar(7,6)=avar(7,6)+ddy2*ddx2
avar(7,7)=avar(7,7)+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) 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 .and. jend==1 ) THEN
DO jazim=kntazim(kk),jmirror,-1
kinbox=0
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(i)
ddy=ryvol(igate,jazim,kk)-ys(j)
ddz=rzvol(igate,jazim,kk)-zps(i,j,k)
!
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)*ddz
rhsvar(5)=rhsvar(5)+varvol(igate,jazim,kk)*ddxy
rhsvar(6)=rhsvar(6)+varvol(igate,jazim,kk)*ddx2
rhsvar(7)=rhsvar(7)+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)+ddz
avar(1,5)=avar(1,5)+ddxy
avar(1,6)=avar(1,6)+ddx2
avar(1,7)=avar(1,7)+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*ddz
avar(2,5)=avar(2,5)+ddx*ddxy
avar(2,6)=avar(2,6)+ddx*ddx2
avar(2,7)=avar(2,7)+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*ddz
avar(3,5)=avar(3,5)+ddy*ddxy
avar(3,6)=avar(3,6)+ddy*ddx2
avar(3,7)=avar(3,7)+ddy*ddy2
!
avar(4,1)=avar(4,1)+ddz
avar(4,2)=avar(4,2)+ddz*ddx
avar(4,3)=avar(4,3)+ddz*ddy
avar(4,4)=avar(4,4)+ddz*ddz
avar(4,5)=avar(4,5)+ddz*ddxy
avar(4,6)=avar(4,6)+ddz*ddx2
avar(4,7)=avar(4,7)+ddz*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*ddz
avar(5,5)=avar(5,5)+ddx2*ddxy
avar(5,6)=avar(5,6)+ddx2*ddx2
avar(5,7)=avar(5,7)+ddx2*ddy2
!
avar(6,1)=avar(6,1)+ddx2
avar(6,2)=avar(6,2)+ddx2*ddx
avar(6,3)=avar(6,3)+ddx2*ddy
avar(6,4)=avar(6,4)+ddx2*ddz
avar(6,5)=avar(6,5)+ddx2*ddxy
avar(6,6)=avar(6,6)+ddx2*ddx2
avar(6,7)=avar(6,7)+ddx2*ddy2
!
avar(7,1)=avar(7,1)+ddy2
avar(7,2)=avar(7,2)+ddy2*ddx
avar(7,3)=avar(7,3)+ddy2*ddy
avar(7,4)=avar(7,4)+ddy2*ddz
avar(7,5)=avar(7,5)+ddy2*ddxy
avar(7,6)=avar(7,6)+ddy2*ddx2
avar(7,7)=avar(7,7)+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) 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 > 7 ) THEN
varmean=rhsvar(1)/avar(1,1)
CALL GJELIM(n,avar,rhsvar,sol,work,work1d,eps,istatus)
gridvar(i,j,k)=min(varmax,max(varmin,sol(1)))
! write(6,'(3i3,i7,4f7.1)') i,j,k,knt, &
! varmin,varmean,varmax,gridvar(i,j,k)
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)
gridvar(i,j,k)=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)
gridvar(i,j,k)=varmean
! write(6,'(3i3,a,i6,6f7.1)') i,j,k,'*',knt, &
! varmin,varmean,varmax,gridvar(i,j,k)
END IF
END IF
END IF
END DO ! i loop
END DO ! j loop
END DO ! k loop
deallocate(elevmean,stat=istatal)
deallocate(varsort,stat=istatal)
IF( vardump == 1 ) &
CALL wrtvar1(nx,ny,nz,gridvar,varid,varname,varunits, &
time,runname,dirname,istatwrt)
RETURN
END SUBROUTINE remapvol
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GJELIM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gjelim(n,a,rhs,sol,work,work1d,eps,istatus) 5
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Solve an N x N system of equations. [a][sol]=[rhs]
! Using Gauss-Jordan with array pivoting.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
! 09/26/00
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! n Matrix order
! a Array
! rhs Right-hand side vector of matrix equation
! eps Error checking threshold
!
! OUTPUT:
! sol Solution vector
! istatus Status indicator
!
! WORK SPACE:
! work Work Array
! work1d Work vector
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: a(n,n)
REAL, INTENT(IN) :: rhs(n)
REAL, INTENT(OUT) :: sol(n)
REAL, INTENT(OUT) :: work(n,n+1)
REAL, INTENT(OUT) :: work1d(n+1)
REAL, INTENT(IN) :: eps
REAL, INTENT(OUT) :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. Local Variables
!
!-----------------------------------------------------------------------
!
REAL :: pivot,const
INTEGER np1
INTEGER :: i,j,k,m
!
!-----------------------------------------------------------------------
!
! Initialize the work array
! First set all elements to zero.
! Fill nxn with elements from input a
! Fill last column with RHS vector.
!
!-----------------------------------------------------------------------
!
np1=n+1
DO j=1, np1
DO i=1, n
work(i,j)=0.0
END DO
END DO
DO j=1, n
DO i=1, n
work(i,j)=a(i,j)
END DO
END DO
DO i=1,n
work(i,np1)=rhs(i)
END DO
DO j=1, n
!
!-----------------------------------------------------------------------
!
! Find largest element in column j
!
!-----------------------------------------------------------------------
!
m=j
pivot=ABS(work(m,j))
DO i=j+1,n
IF(ABS(work(i,j)) > pivot ) THEN
m=i
pivot=ABS(work(m,j))
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Error trapping
!
!-----------------------------------------------------------------------
!
IF( pivot < eps ) THEN
DO i=1, n
sol(i)=0.
END DO
istatus=-1
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Swap rows
!
!-----------------------------------------------------------------------
!
IF(m /= j) THEN
DO k=1,np1
work1d(k)=work(j,k)
END DO
DO k=1,np1
work(j,k)=work(m,k)
work(m,k)=work1d(k)
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Normalize Row
!
!-----------------------------------------------------------------------
!
const=1./work(j,j)
DO k=1,np1
work(j,k)=const*work(j,k)
END DO
work(j,j)=1.0
!
!-----------------------------------------------------------------------
!
! Elimination
!
!-----------------------------------------------------------------------
!
DO i=1,n
IF ( i /= j ) THEN
const=work(i,j)
DO k=1,np1
work(i,k)=work(i,k)-const*work(j,k)
END DO
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Transfer last column to sol vector
!
!-----------------------------------------------------------------------
!
DO i=1,n
sol(i)=work(i,n+1)
END DO
istatus = 1
RETURN
END SUBROUTINE gjelim