!########################################################################
!########################################################################
!######### #########
!######### PROGRAM ncrad2arps #########
!######### #########
!######### Developed by #########
!######### Center for Analysis and Prediction of Storms #########
!######### University of Oklahoma #########
!######### #########
!########################################################################
!########################################################################
PROGRAM ncrad2arps,50
!------------------------------------------------------------------------
!
! PURPOSE:
!
! Reads radar data from NetCDF data files and remaps data on ARPS
! grid to be used for analysis or display.
! CASA radar data, for example, are stored in netCDF files.
!
!------------------------------------------------------------------------
!
! AUTHOR:
!
! Keith Brewster (May 2006)
!
! MODIFICATIONS:
!
! 1 Feb 2007 Keith Brewster
! Added processing for gzipped files.
!
!------------------------------------------------------------------------
IMPLICIT NONE
!------------------------------------------------------------------------
!
! Include files
!
!------------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'netcdf.inc'
!------------------------------------------------------------------------
!
! Dimensions and parameters
!
!------------------------------------------------------------------------
INTEGER, PARAMETER :: maxfile = 90
INTEGER, PARAMETER :: maxvar = 10
!
INTEGER, PARAMETER :: stdout = 6
INTEGER, PARAMETER :: itime1970 = 315619200
REAL, PARAMETER :: refchek=-35.
REAL, PARAMETER :: refmiss=-999.
REAL, PARAMETER :: velchek=-190.
REAL, PARAMETER :: velmiss=-999.
REAL, PARAMETER :: winszazim=10.
REAL, PARAMETER :: winszrad=1000.
!------------------------------------------------------------------------
!
! Variable Declarations.
!
!------------------------------------------------------------------------
INTEGER :: nx, ny, nz, nzsoil, itime, nstyps
REAL :: radelv, radlat, radlon
REAL :: radarx, radary
!------------------------------------------------------------------------
!
! CASA Radar Data Variables
!
!------------------------------------------------------------------------
!
REAL, ALLOCATABLE :: refl(:,:)
REAL, ALLOCATABLE :: attref(:,:)
REAL, ALLOCATABLE :: rhohv(:,:)
REAL, ALLOCATABLE :: radv(:,:)
!
!------------------------------------------------------------------------
!
! Variables for subroutine remapvol.
!
!------------------------------------------------------------------------
!
INTEGER, PARAMETER :: nt=1
INTEGER, PARAMETER :: exbcbufsz=1
INTEGER, PARAMETER :: bkgopt=1
INTEGER, PARAMETER :: shropt=1
INTEGER, PARAMETER :: rfropt = 1 ! 4/3rds earth refraction option
INTEGER, PARAMETER :: maxsort = 10000
INTEGER, PARAMETER :: nzsnd = 201
REAL, PARAMETER :: dzsnd=100.
INTEGER :: ktsnd(nzsnd)
REAL :: zsnd(nzsnd) ! hgt levels for refractivity sounding
REAL :: usnd(nzsnd)
REAL :: vsnd(nzsnd)
REAL :: rfrsnd(nzsnd) ! refractivity sounding
REAL, PARAMETER :: rfrconst = (-1./(4.*6371.E03)) ! -1/(4e)
REAL, PARAMETER :: envavgr=80000.
REAL, ALLOCATABLE :: arfdata(:,:) ! attenuated refl (dBZ) data
REAL, ALLOCATABLE :: refdata(:,:) ! refl (dBZ) data
REAL, ALLOCATABLE :: veldata(:,:) ! velocity data (m/s)
REAL, ALLOCATABLE :: spwdata(:,:)
REAL, ALLOCATABLE :: bkgvel(:,:)
REAL, ALLOCATABLE :: unfvdata(:,:) ! unfolded velocity data (m/s)
REAL, ALLOCATABLE :: rhodata(:,:) ! rho-HV data
REAL, ALLOCATABLE :: snstvty(:) ! radar sensitivity
REAL, ALLOCATABLE :: rtem(:,:) ! temporary array for radar data processing
INTEGER, ALLOCATABLE :: time(:) ! time offset from time1st
REAL, ALLOCATABLE :: azim(:) ! azimuth angle for each radial(degree)
REAL, ALLOCATABLE :: beamw(:) ! beamwidth for each radial(degree)
REAL, ALLOCATABLE :: gtspc(:) ! gate spacing for each radial (mneter)
REAL, ALLOCATABLE :: vnyq(:) ! Nyquist velocity (meters/sec)
REAL, ALLOCATABLE :: elev(:) ! elevation angle for each radial (degree)
INTEGER, ALLOCATABLE :: istrgate(:)
INTEGER, ALLOCATABLE :: bgate(:)
INTEGER, ALLOCATABLE :: egate(:)
INTEGER :: time1st ! time of first radial in volume
INTEGER :: irfirstg ! range to first gate (meters)
INTEGER :: igatesp ! gate spacing (meters)
INTEGER :: irefgatsp ! reflectivity gate spacing (meters)
INTEGER :: ivelgatsp ! velocity gate spacing (meters)
REAL :: elv ! elevation angle for this tilt
REAL :: rfirstg ! range to first gate (m)
REAL :: gatesp ! gate spacing (meters)
REAL :: vnyquist ! Nyquist velocity for this tilt
INTEGER, ALLOCATABLE :: kntrgat(:,:)
INTEGER, ALLOCATABLE :: kntrazm(:)
INTEGER :: kntrelv
INTEGER, ALLOCATABLE :: kntvgat(:,:)
INTEGER, ALLOCATABLE :: kntvazm(:)
INTEGER :: kntvelv
INTEGER, ALLOCATABLE :: timevol(:,:)
REAL, ALLOCATABLE :: nyqvvol(:)
REAL, ALLOCATABLE :: rngrvol(:,:)
REAL, ALLOCATABLE :: azmrvol(:,:)
REAL, ALLOCATABLE :: elvrvol(:,:)
REAL, ALLOCATABLE :: refvol(:,:,:)
REAL, ALLOCATABLE :: rngvvol(:,:)
REAL, ALLOCATABLE :: azmvvol(:,:)
REAL, ALLOCATABLE :: elvvvol(:,:)
REAL, ALLOCATABLE :: velvol(:,:,:)
REAL, ALLOCATABLE :: elvmean(:)
REAL, ALLOCATABLE :: varsort(:)
REAL, ALLOCATABLE :: rxvol(:,:,:)
REAL, ALLOCATABLE :: ryvol(:,:,:)
REAL, ALLOCATABLE :: rzvol(:,:,:)
REAL, ALLOCATABLE :: gridvel(:,:,:)
REAL, ALLOCATABLE :: gridref(:,:,:)
REAL, ALLOCATABLE :: gridnyq(:,:,:)
REAL, ALLOCATABLE :: gridtim(:,:,:)
REAL, ALLOCATABLE :: x(:)
REAL, ALLOCATABLE :: y(:)
REAL, ALLOCATABLE :: z(:)
REAL, ALLOCATABLE :: zp(:,:,:)
REAL, ALLOCATABLE :: xs(:)
REAL, ALLOCATABLE :: ys(:)
REAL, ALLOCATABLE :: zps(:,:,:)
REAL, ALLOCATABLE :: hterain(:,:)
REAL, ALLOCATABLE :: mapfct(:,:,:) ! Map factors at scalar, u and v points
REAL, ALLOCATABLE :: j1 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( x ).
REAL, ALLOCATABLE :: j2 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( y ).
REAL, ALLOCATABLE :: j3 (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ).
REAL, ALLOCATABLE :: tem1d1(:)
REAL, ALLOCATABLE :: tem1d2(:)
REAL, ALLOCATABLE :: tem2d(:,:)
REAL, ALLOCATABLE :: tem2dyz(:,:)
REAL, ALLOCATABLE :: tem2dxz(:,:)
REAL, ALLOCATABLE :: tem2dns(:,:,:)
REAL, ALLOCATABLE :: tem3d1(:,:,:)
REAL, ALLOCATABLE :: tem3d2(:,:,:)
REAL, ALLOCATABLE :: tem3d3(:,:,:)
REAL, ALLOCATABLE :: tem3d4(:,:,:)
INTEGER, ALLOCATABLE :: tem2dint(:,:)
INTEGER, ALLOCATABLE :: soiltyp(:,:,:)
REAL, ALLOCATABLE :: stypfrct(:,:,:)
REAL, ALLOCATABLE :: tem3dsoil(:,:,:)
REAL, ALLOCATABLE :: tem4dsoilns(:,:,:,:)
REAL, ALLOCATABLE :: zpsoil(:,:,:)
REAL, ALLOCATABLE :: j3soil(:,:,:)
REAL, ALLOCATABLE :: j3soilinv(:,:,:)
REAL, ALLOCATABLE :: u(:,:,:)
REAL, ALLOCATABLE :: v(:,:,:)
REAL, ALLOCATABLE :: qv(:,:,:)
REAL, ALLOCATABLE :: ubar(:,:,:)
REAL, ALLOCATABLE :: vbar(:,:,:)
REAL, ALLOCATABLE :: ptbar(:,:,:)
REAL, ALLOCATABLE :: uprt(:,:,:)
REAL, ALLOCATABLE :: vprt(:,:,:)
REAL, ALLOCATABLE :: ptprt(:,:,:)
REAL, ALLOCATABLE :: pprt(:,:,:)
REAL, ALLOCATABLE :: pbar(:,:,:)
REAL, ALLOCATABLE :: qvbar(:,:,:)
REAL, ALLOCATABLE :: rhostr(:,:,:)
REAL, ALLOCATABLE :: trigs1(:)
REAL, ALLOCATABLE :: trigs2(:)
INTEGER, ALLOCATABLE :: ifax(:)
REAL, ALLOCATABLE :: wsave1(:)
REAL, ALLOCATABLE :: wsave2(:)
REAL, ALLOCATABLE :: vwork1(:,:)
REAL, ALLOCATABLE :: vwork2(:,:)
REAL, ALLOCATABLE :: qcumsrc(:,:,:,:)
REAL, ALLOCATABLE :: prcrate(:,:,:)
REAL, ALLOCATABLE :: exbcbuf(:)
REAL, ALLOCATABLE :: grdvel2(:,:,:)
REAL, ALLOCATABLE :: grdref2(:,:,:)
REAL, ALLOCATABLE :: wgtvel(:,:,:)
REAL, ALLOCATABLE :: wgtref(:,:,:)
REAL, ALLOCATABLE :: wpotvel(:,:,:)
REAL, ALLOCATABLE :: wpotref(:,:,:)
REAL, ALLOCATABLE :: gridazm(:,:)
REAL, ALLOCATABLE :: gridrng(:,:)
REAL(KIND=8), ALLOCATABLE :: tem_double(:)
!------------------------------------------------------------------------
!
! Data dimensions, determined from the data files
!
!------------------------------------------------------------------------
INTEGER :: maxgate,maxazim,maxelev
!------------------------------------------------------------------------
!
! Misc. variables
!
!------------------------------------------------------------------------
CHARACTER(LEN=256) :: infilname
CHARACTER(LEN=80) :: indir
CHARACTER(LEN=80) :: dir_name
CHARACTER(LEN=256) :: fname
CHARACTER(LEN=256) :: full_fname
CHARACTER(LEN=6) :: velid = 'radv3d'
CHARACTER(LEN=6) :: refid = 'refl3d'
CHARACTER(LEN=20) :: refname = 'Reflectivity'
CHARACTER(LEN=20) :: velname = 'RawVelocity'
CHARACTER(LEN=6) :: refunits = 'dBZ'
CHARACTER(LEN=6) :: velunits = 'm/s'
CHARACTER(LEN=NF_MAX_NAME), ALLOCATABLE :: ncvarname(:)
CHARACTER(LEN=266) :: syscall
INTEGER :: fvartype(maxfile)
LOGICAL :: gzipped(maxfile)
REAL :: felv(maxfile)
INTEGER :: fitime(maxfile)
INTEGER :: n, i, j, k, istatus
INTEGER :: nazim,ngate,nvar
INTEGER :: ivcp,itimcdf,initime,mintimcdf
INTEGER :: ivar,varfill,dmpfmt,hdf4cmpr
INTEGER :: len_dir,len_fname
INTEGER :: arbvaropt,velopt
LOGICAL :: velproc
LOGICAL :: vel2d,vel3d
LOGICAL :: ref2d,ref3d
LOGICAL :: qc_on
LOGICAL :: rho_match
INTEGER :: iyear,iyr,imonth,iday,ihr,imin,isec
INTEGER :: iiyr,iimonth,iiday,iihr,iimin,iisec
INTEGER :: velknt,refknt,rhoknt,kntt2a,kntt2b
INTEGER :: isource,icount,iradtype
INTEGER :: nyqset,timeset,vardump
INTEGER :: xscale
REAL :: rdx,dtime,rmin,rmax
REAL :: rmisval,rngfval,frtime
REAL :: refelvmin,refelvmax
REAL :: refrngmin,refrngmax
INTEGER :: remapopt
REAL :: radius,vconst,bmwidth
REAL :: xrad,yrad
INTEGER :: maxinfile,infile,lenfn
INTEGER :: iarg,jarg,iargc,narg,ifile,jfile,kfile
CHARACTER(LEN=4) tsthead
CHARACTER(LEN=256) charg
CHARACTER(LEN=256) listfile
CHARACTER (LEN=6) :: varid
CHARACTER (LEN=20) :: varname
CHARACTER (LEN=20) :: varunits
!------------------------------------------------------------------------
!
! Some tunable parameters - convert to Namelist in future
!
!------------------------------------------------------------------------
INTEGER, PARAMETER :: maxncradfile=90
INTEGER, PARAMETER :: iordref = 2
INTEGER, PARAMETER :: iordvel = 1
REAL, PARAMETER :: refmedl = 15.
REAL, PARAMETER :: velmedl = 15.
REAL, PARAMETER :: refdazl = 360.
REAL, PARAMETER :: veldazl = 20.
REAL, PARAMETER :: dazim = 1.0
REAL, PARAMETER :: rngmin = 500.
REAL, PARAMETER :: rngmax = 230.E03
REAL, PARAMETER :: rhohvlim = 0.50
CHARACTER(LEN=4) :: radname
CHARACTER(LEN=32) :: radarname
INTEGER :: nncradfn = 0
CHARACTER(LEN=256) :: ncradfn(maxncradfile)
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!------------------------------------------------------------------------
!
! Initializations
!
!------------------------------------------------------------------------
isource = 3 ! 1 - WSR-88D raw 2 - WSR-88D NIDS
refelvmin = 91.0
refelvmax = -91.0
refrngmin = rngmin
refrngmax = rngmax
dmpfmt=1
hdf4cmpr=0
kntt2a=0
kntt2b=0
refknt=0
rhoknt=0
velknt=0
velproc = .TRUE.
qc_on = .TRUE.
ref2d = .FALSE.
ref3d = .FALSE.
vel2d = .FALSE.
vel3d = .FALSE.
!------------------------------------------------------------------------
!
! Initialize ncrad_data NAMELIST variables
!
!------------------------------------------------------------------------
radname = 'DMMY'
nncradfn = 0
dir_name = './'
len_dir = LEN(TRIM(dir_name))
arbvaropt = 0
DO ifile=1,maxncradfile
ncradfn(ifile)='dummy'
END DO
narg=iargc()
WRITE(stdout,'(a,i4)') ' Number of command-line arguments: ',narg
IF(narg < 2 ) THEN
WRITE(stdout,'(a,a,a,a)') &
' Usage: ncrad2arps RADAR_ID radarfile_listfile [-novel] [-hdf 2]', &
' [-binary] [-ref2d] [-ref3d] [-reffile] [-vel2d] [-vel3d] [-velfile]',&
' [-noqc]', &
' < arps.input'
STOP
END IF
!------------------------------------------------------------------------
!
! Read grid info from the arps.input
!
!------------------------------------------------------------------------
CALL initpara
(nx, ny, nz, nzsoil, nstyps)! Reads standard ARPS namelists
!------------------------------------------------------------------------
!
! Process command line
! For backward comptability allow input of radar file names via
! files specified on the command line
!
!------------------------------------------------------------------------
IF(narg > 1) THEN
CALL getarg(1,charg)
radname=charg(1:4)
WRITE(stdout,'(a,a)') ' radname = ', radname
kfile=0
nncradfn=0
iarg=2
DO jarg=2,narg
CALL getarg(iarg,charg)
IF(charg(1:6) == '-novel') THEN
velproc=.FALSE.
ELSE IF(charg(1:4) == '-hdf') THEN
dmpfmt=2
hdf4cmpr=0
IF(iarg < narg) THEN
iarg=iarg+1
CALL getarg(iarg,charg)
READ(charg,'(i1)') hdf4cmpr
hdf4cmpr=min(max(hdf4cmpr,0),7)
END IF
WRITE(stdout,'(a,i2)') &
' Output in hdf format with compression level: ',hdf4cmpr
ELSE IF(charg(1:7) == '-binary') THEN
dmpfmt=1
hdf4cmpr=0
WRITE(stdout,'(a)') ' Output in binary format'
ELSE IF(charg(1:6) == '-ref2d') THEN
ref2d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 2d reflectivity file of lowest tilt'
ELSE IF(charg(1:6) == '-ref3d') THEN
ref3d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 3d reflectivity file for plotting'
ELSE IF(charg(1:8) == '-reffile') THEN
ref3d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 3d reflectivity file for plotting'
ELSE IF(charg(1:6) == '-vel2d') THEN
vel2d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 2d velocity file of lowest tilt'
ELSE IF(charg(1:6) == '-vel3d') THEN
vel3d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 3d velocity file for plotting'
ELSE IF(charg(1:8) == '-velfile') THEN
vel3d=.TRUE.
WRITE(stdout,'(a)') &
' Will produce 3d velocity file for plotting'
ELSE IF(charg(1:5) == '-noqc') THEN
qc_on=.FALSE.
WRITE(stdout,'(a)') &
' Will skip qc steps in processing'
ELSE
listfile=charg
WRITE(stdout,'(a,a)') ' Reading file lists ',charg
OPEN(31,file=listfile,status='old',form='formatted')
DO ifile=(nncradfn+1),maxncradfile
READ(31,'(a)',END=101) ncradfn(ifile)
kfile=kfile+1
WRITE(stdout,'(a,i4,a,a)') &
' kfile:',kfile,' ncradfn: ',TRIM(ncradfn(ifile))
END DO
101 CONTINUE
nncradfn=kfile
END IF
iarg=iarg+1
IF(iarg > narg) EXIT
END DO
WRITE(stdout,'(i4,a)') nncradfn,' file names read'
END IF
istatus = 0
icount = 0
rfirstg=1000.
gatesp = 1000.
irfirstg = NINT(rfirstg)
igatesp = NINT(gatesp)
!------------------------------------------------------------------------
!
! Loop through the NetCDF radar files.
!
!------------------------------------------------------------------------
IF (nncradfn > 0 .AND. nncradfn <= maxncradfile) THEN
maxinfile = nncradfn
ELSE IF (nncradfn > maxncradfile) THEN
WRITE(stdout,'(a,i4,a,i4)') &
' WARNING: nncradfn=', nncradfn,' > maxncradfile=',maxncradfile
WRITE(stdout,'(a,i4,a)') &
' Will only read in ', maxncradfile,' NIDS files.'
maxinfile = maxncradfile
ELSE
maxinfile = 0
END IF
IF (maxinfile > 0) THEN
!
!------------------------------------------------------------------------
!
! First check each file for variable stored and dimensions of data
!
!------------------------------------------------------------------------
!
maxgate=1
maxazim=1
maxelev=1
mintimcdf=-1
DO infile = 1,maxinfile
fvartype(infile)=0
gzipped(infile)=.FALSE.
END DO
DO infile = 1,maxinfile
nazim = 0
ngate = 0
nvar = 0
infilname = ncradfn(infile)
lenfn=len(infilname)
CALL strlnth(infilname,lenfn)
IF(infilname((lenfn-2):lenfn) == ".gz") THEN
gzipped(infile)=.TRUE.
WRITE(syscall,'(a,1x,a)') 'gunzip',TRIM(infilname)
CALL system(syscall)
infilname((lenfn-2):lenfn)=' '
ncradfn(infile)=infilname
WRITE(6,'(a,a)') 'Uncompressed ',TRIM(infilname)
END IF
WRITE(stdout,'(a,a)') 'Checking dims in file:',TRIM(infilname)
CALL extdirn(infilname,indir,fname)
print *, ' infile: ',TRIM(infilname)
print *, ' indir: ',TRIM(indir)
print *, ' fname: ',TRIM(fname)
CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, &
istatus)
WRITE(6,'(a,i6,a,i6)') ' istatus: ',istatus,' nvar in file: ',nvar
IF(istatus == NF_NOERR .AND. nvar > 0) THEN
ALLOCATE(ncvarname(nvar))
IF(iradtype == 21) THEN
kntt2a=kntt2a+1
ALLOCATE(tem_double(nazim))
CALL get_ncrad2ainfo(fname,indir,nazim,nvar,tem_double, &
radarname,radlat,radlon,radelv,itimcdf,frtime,elv, &
ncvarname,istatus)
DEALLOCATE(tem_double)
fvartype(infile)=100
refknt=refknt+1
velknt=velknt+1
rhoknt=rhoknt+1
ELSE IF(iradtype == 22) THEN
kntt2b=kntt2b+1
CALL get_ncrad2binfo(fname,indir,nazim,nvar, &
radarname,radlat,radlon,radelv,itimcdf,frtime, &
ivcp,elv,ncvarname,istatus)
DO k=1,nvar
print *, ' variable(',k,') = ',TRIM(ncvarname(k))
IF(ncvarname(k)(1:7) == 'Reflect') THEN
fvartype(infile)=1
refknt=refknt+1
EXIT
ELSE IF(ncvarname(k)(1:11) == 'RawVelocity' .OR. &
ncvarname(k)(1:15) == 'AliasedVelocity' ) THEN
fvartype(infile)=2
velknt=velknt+1
EXIT
ELSE IF(ncvarname(k)(1:5) == 'RhoHV') THEN
fvartype(infile)=3
rhoknt=rhoknt+1
EXIT
END IF
END DO
ELSE
WRITE(stdout,'(a,a)') 'Unknown radar file type:',iradtype
fvartype(infile)=0
STOP
END IF
IF(radarname(1:5) == 'cyril') THEN
IF(radelv < 100.) radelv=432.9
IF(elv < 0.1) elv=1.0
ELSE IF(radarname(1:5) == 'chick') THEN
IF(radelv < 100.) radelv=354.0
IF(elv < 0.1) elv=1.0
ELSE IF(radarname(1:5) == 'lawto') THEN
IF(radelv < 100.) radelv=377.4
IF(elv < 0.1) elv=1.0
ELSE IF(radarname(1:5) == 'rushs') THEN
IF(radelv < 100.) radelv=414.8
IF(elv < 0.1) elv=1.0
END IF
felv(infile)=elv
fitime(infile)=itimcdf
maxazim=max(nazim,maxazim)
maxgate=max(ngate,maxgate)
IF(mintimcdf < 0) THEN
mintimcdf=itimcdf
ELSE
mintimcdf=min(mintimcdf,itimcdf)
END IF
DEALLOCATE(ncvarname)
ELSE ! error return from get dims
fvartype(infile)=0
END IF
END DO
maxelev=max(maxelev,refknt)
maxelev=max(maxelev,velknt)
IF(velknt == 0) velproc=.false.
WRITE(6,'(/a,i6)') ' Number of reflectivity files:',refknt
WRITE(6,'(a,i6)') ' Number of Rho-HV files:',rhoknt
WRITE(6,'(a,i6)') ' Number of velocity files:',velknt
WRITE(6,'(/a,i6)') ' Max azimuth dimension of NetCDF data:',maxazim
WRITE(6,'(a,i6)') ' Max gate dimension of NetCDF data:',maxgate
WRITE(6,'(a,i6/)') ' Max elevations among data:',maxelev
IF(kntt2a > 0 .AND. kntt2b > 0) THEN
WRITE(6,'(a)') ' Mixed Tier 2a and Tier 2b files in list, stopping'
WRITE(6,'(a,i4,a,i4)')' N Tier 2a=',kntt2a,' N Tier 2b=',kntt2b
STOP
END IF
ALLOCATE(x(nx))
ALLOCATE(y(ny))
ALLOCATE(z(nz))
ALLOCATE(zp(nx,ny,nz))
ALLOCATE(hterain(nx,ny))
ALLOCATE(mapfct(nx,ny,8))
ALLOCATE(j1(nx,ny,nz))
ALLOCATE(j2(nx,ny,nz))
ALLOCATE(j3(nx,ny,nz))
ALLOCATE(tem1d1(nz))
ALLOCATE(tem1d2(nz))
ALLOCATE(tem3d1(nx,ny,nz))
IF(qc_on .AND. velknt > 0) THEN
!------------------------------------------------------------------------
!
! Fill refractivity sounding with constant lapse rate
! This is only used when rfropt > 1, but for completeness and
! future upgrade (when an actual sounding made from gridded data
! would be used) this is filled-in.
!
!------------------------------------------------------------------------
DO k=1,nzsnd
zsnd(k)=(k-1)*dzsnd
rfrsnd(k)=325.+rfrconst*zsnd(k)
END DO
!------------------------------------------------------------------------
!
! ARPS grid array allocations and variable initializations.
!
!------------------------------------------------------------------------
ALLOCATE(u(nx,ny,nz))
ALLOCATE(v(nx,ny,nz))
ALLOCATE(qv(nx,ny,nz))
ALLOCATE(uprt(nx,ny,nz))
ALLOCATE(vprt(nx,ny,nz))
ALLOCATE(ptprt(nx,ny,nz))
ALLOCATE(pprt(nx,ny,nz))
ALLOCATE(ubar(nx,ny,nz))
ALLOCATE(vbar(nx,ny,nz))
ALLOCATE(ptbar(nx,ny,nz))
ALLOCATE(pbar(nx,ny,nz))
ALLOCATE(qvbar(nx,ny,nz))
ALLOCATE(rhostr(nx,ny,nz))
ALLOCATE(tem2d(nx,ny))
ALLOCATE(tem2dyz(ny,nz))
ALLOCATE(tem2dxz(nx,nz))
ALLOCATE(tem2dns(nx,ny,0:nstyps))
ALLOCATE(trigs1(3*(nx-1)/2+1))
ALLOCATE(trigs2(3*(ny-1)/2+1))
ALLOCATE(ifax(13))
ALLOCATE(wsave1(3*(ny-1)+15))
ALLOCATE(wsave2(3*(nx-1)+15))
ALLOCATE(vwork1(nx+1,ny+1))
ALLOCATE(vwork2(ny,nx+1))
ALLOCATE(qcumsrc(nx,ny,nz,5))
ALLOCATE(prcrate(nx,ny,4))
ALLOCATE(exbcbuf(exbcbufsz))
ALLOCATE(soiltyp(nx,ny,nstyps))
ALLOCATE(stypfrct(nx,ny,nstyps))
ALLOCATE(tem2dint(nx,ny))
ALLOCATE(tem3d2(nx,ny,nz))
ALLOCATE(tem3d3(nx,ny,nz))
ALLOCATE(tem3d4(nx,ny,nz))
ALLOCATE(tem3dsoil(nx,ny,nzsoil))
ALLOCATE(tem4dsoilns(nx,ny,nzsoil,0:nstyps))
CALL set_lbcopt(1)
CALL initgrdvar(nx,ny,nz,nzsoil,nt,nstyps,exbcbufsz, &
x,y,z,zp,tem3dsoil,hterain,mapfct, &
j1,j2,j3,tem3dsoil,tem3d1,tem3d1,tem3d1,tem3d1,tem3dsoil, &
u,v,tem3d1,tem3d1,ptprt,pprt, &
qv,tem3d1,tem3d1,tem3d1,tem3d1,tem3d1,tem3d1, &
tem2dyz,tem2dyz,tem2dxz,tem2dxz, &
tem2dyz,tem2dyz,tem2dxz,tem2dxz, &
trigs1,trigs2,ifax,ifax, &
wsave1,wsave2,vwork1,vwork2, &
ubar,vbar,ptbar,pbar,tem3d1,tem3d1, &
rhostr,tem3d1,qvbar,tem3d1,tem3d1, &
soiltyp,stypfrct,tem2dint,tem2d,tem2d,tem2d, &
tem4dsoilns,tem4dsoilns,tem2dns,tem2d,tem2dns, &
tem3d1,qcumsrc,tem3d1,tem2dint,tem2d, &
tem2d,tem2d,tem2d, &
tem2d,tem2d,prcrate,exbcbuf, &
tem3d1,tem2d,tem2d,tem2d,tem2d, &
tem2d,tem2d,tem2d,tem2d, &
tem3dsoil,tem3dsoil,tem3dsoil,tem3dsoil,tem3dsoil, &
tem3d2,tem3d3,tem3d4,tem3d4,tem3d4, &
tem3d4,tem3d4,tem3d4,tem3d4)
print *, ' Back from initgrdvar'
DEALLOCATE(j1)
DEALLOCATE(j2)
DEALLOCATE(j3)
DEALLOCATE(trigs1)
DEALLOCATE(trigs2)
DEALLOCATE(ifax)
DEALLOCATE(wsave1)
DEALLOCATE(wsave2)
DEALLOCATE(vwork1)
DEALLOCATE(vwork2)
DEALLOCATE(qcumsrc)
DEALLOCATE(prcrate)
DEALLOCATE(soiltyp)
DEALLOCATE(stypfrct)
DEALLOCATE(exbcbuf)
DEALLOCATE(tem2dint)
DEALLOCATE(tem2dyz)
DEALLOCATE(tem2dxz)
DEALLOCATE(tem2dns)
DEALLOCATE(tem4dsoilns)
DEALLOCATE(tem3dsoil)
ALLOCATE(xs(nx))
ALLOCATE(ys(ny))
ALLOCATE(zps(nx,ny,nz))
print *, ' Setting radar coordinate: '
CALL radcoord(nx,ny,nz,x,y,z,zp,xs,ys,zps, &
radlat,radlon,radarx,radary)
print *, ' Radar x: ',(0.001*radarx),' km'
print *, ' Radar y: ',(0.001*radary),' km'
print *, ' Environmental wind averaging radius: ',envavgr
CALL extenvprf(nx,ny,nz,nzsnd, &
x,y,zp,xs,ys,zps, &
u,v,ptprt,pprt,qv,ptbar,pbar,tem3d1,tem3d2,tem3d3, &
radarx,radary,radelv,envavgr, &
zsnd,ktsnd,usnd,vsnd,rfrsnd)
print *, ' Back from extenvprf'
DEALLOCATE(ubar)
DEALLOCATE(vbar)
DEALLOCATE(pbar)
DEALLOCATE(ptbar)
DEALLOCATE(qvbar)
DEALLOCATE(rhostr)
DEALLOCATE(uprt)
DEALLOCATE(vprt)
DEALLOCATE(ptprt)
DEALLOCATE(pprt)
DEALLOCATE(qv)
DEALLOCATE(tem3d2)
DEALLOCATE(tem3d3)
DEALLOCATE(tem3d4)
ELSE
ALLOCATE(xs(nx))
ALLOCATE(ys(ny))
ALLOCATE(zps(nx,ny,nz))
ALLOCATE(zpsoil(nx,ny,nzsoil))
ALLOCATE(j3soil(nx,ny,nzsoil))
ALLOCATE(j3soilinv(nx,ny,nzsoil))
print *, ' Setting ARPS grid coordinates ...'
CALL inigrd(nx,ny,nz,nzsoil,x,y,z,zp,zpsoil, &
hterain,mapfct,j1,j2,j3,j3soil,j3soilinv, &
tem1d1,tem1d2,tem3d1)
print *, ' Setting radar coordinate...'
CALL radcoord(nx,ny,nz,x,y,z,zp,xs,ys,zps, &
radlat,radlon,radarx,radary)
DEALLOCATE(j1)
DEALLOCATE(j2)
DEALLOCATE(j3)
DEALLOCATE(zpsoil)
DEALLOCATE(j3soil)
DEALLOCATE(j3soilinv)
END IF
DEALLOCATE(hterain)
DEALLOCATE(mapfct)
DEALLOCATE(tem3d1)
ALLOCATE(gridvel(nx,ny,nz))
ALLOCATE(gridref(nx,ny,nz))
ALLOCATE(gridnyq(nx,ny,nz))
ALLOCATE(gridtim(nx,ny,nz))
ALLOCATE(grdvel2(nx,ny,nz))
ALLOCATE(grdref2(nx,ny,nz))
ALLOCATE(wgtvel(nx,ny,nz))
ALLOCATE(wgtref(nx,ny,nz))
ALLOCATE(wpotvel(nx,ny,nz))
ALLOCATE(wpotref(nx,ny,nz))
ALLOCATE(gridazm(nx,ny))
ALLOCATE(gridrng(nx,ny))
ALLOCATE(varsort(maxsort))
!------------------------------------------------------------------------
!
! Radar data array allocations
!
!------------------------------------------------------------------------
IF(refknt > 0) THEN
ALLOCATE(refdata(maxgate,maxazim))
ALLOCATE(arfdata(maxgate,maxazim))
END IF
IF(velknt > 0) THEN
ALLOCATE(veldata(maxgate,maxazim))
IF(qc_on) THEN
ALLOCATE(spwdata(maxgate,maxazim))
ALLOCATE(unfvdata(maxgate,maxazim))
ALLOCATE(bkgvel(maxgate,maxazim))
ALLOCATE(bgate(maxazim))
ALLOCATE(egate(maxazim))
spwdata=0.
END IF
END IF
IF(rhoknt > 0) ALLOCATE(rhodata(maxgate,maxazim))
ALLOCATE(snstvty(maxgate))
ALLOCATE(rtem(maxgate,maxazim))
ALLOCATE(time(maxazim))
ALLOCATE(azim(maxazim))
ALLOCATE(beamw(maxazim))
ALLOCATE(gtspc(maxazim))
ALLOCATE(vnyq(maxazim))
ALLOCATE(elev(maxazim))
ALLOCATE(istrgate(maxazim))
ALLOCATE(kntrgat(maxazim,maxelev))
ALLOCATE(kntrazm(maxelev))
ALLOCATE(kntvgat(maxazim,maxelev))
ALLOCATE(kntvazm(maxelev))
ALLOCATE(nyqvvol(maxelev))
ALLOCATE(timevol(maxazim,maxelev))
ALLOCATE(rngrvol(maxgate,maxelev))
ALLOCATE(azmrvol(maxazim,maxelev))
ALLOCATE(elvrvol(maxazim,maxelev))
ALLOCATE(refvol(maxgate,maxazim,maxelev))
ALLOCATE(rngvvol(maxgate,maxelev))
ALLOCATE(azmvvol(maxazim,maxelev))
ALLOCATE(elvvvol(maxazim,maxelev))
ALLOCATE(velvol(maxgate,maxazim,maxelev))
ALLOCATE(elvmean(maxelev))
ALLOCATE(rxvol(maxgate,maxazim,maxelev))
ALLOCATE(ryvol(maxgate,maxazim,maxelev))
ALLOCATE(rzvol(maxgate,maxazim,maxelev))
!------------------------------------------------------------------------
!
! Initialize 3-d radar data arrays
!
!------------------------------------------------------------------------
itimcdf=mintimcdf+itime1970
time1st=itimcdf
CALL abss2ctim(itimcdf, iyear, imonth, iday, ihr, imin, isec)
iyr = MOD(iyear, 100)
CALL rmpinit(nx,ny,nz,maxgate,maxgate,maxazim,maxelev, &
kntrgat,kntrazm,kntrelv, &
kntvgat,kntvazm,kntvelv, &
nyqvvol,timevol, &
rngrvol,azmrvol,elvrvol,refvol, &
rngvvol,azmvvol,elvvvol,velvol, &
gridvel,gridref,gridnyq,gridtim)
iiyr = iyr
iimonth = imonth
iiday = iday
iihr = ihr
iimin = imin
iisec = isec
IF(kntt2a > 0) THEN
DO infile = 1,maxinfile
infilname = ncradfn(infile)
WRITE(6,'(a,i4,a,i4)') ' infile= ',infile,' type=',fvartype(infile)
CALL extdirn(infilname,indir,fname)
CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, &
istatus)
ALLOCATE(tem_double(nazim))
ALLOCATE(ncvarname(nvar))
CALL get_ncrad2ainfo(fname,indir,nazim,nvar,tem_double, &
radarname,radlat,radlon,radelv,itimcdf,frtime,elv, &
ncvarname,istatus)
DEALLOCATE(ncvarname)
ALLOCATE(attref(ngate,nazim))
ALLOCATE(refl(ngate,nazim))
ALLOCATE(radv(ngate,nazim))
ALLOCATE(rhohv(ngate,nazim))
CALL rd2atiltcdf(nazim,ngate,fname,indir, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyquist,rfirstg,istrgate, &
azim,beamw,gtspc,attref,refl,radv,rhohv, &
tem_double)
DEALLOCATE(tem_double)
bmwidth=beamw(1)
arfdata=-999.0
refdata=-999.0
DO j=1,nazim
DO i=istrgate(j),ngate
IF(abs(refl(i,j)) < 99. ) THEN
arfdata(i,j)=attref(i,j)
refdata(i,j)=refl(i,j)
END IF
END DO
END DO
veldata=-999.0
DO j=1,nazim
DO i=istrgate(j),ngate
IF(abs(radv(i,j)) < 199. ) veldata(i,j)=radv(i,j)
END DO
END DO
rhodata=-999.0
DO j=1,nazim
DO i=istrgate(j),ngate
IF(abs(rhohv(i,j)) < 1.1 ) rhodata(i,j)=rhohv(i,j)
END DO
END DO
DEALLOCATE(attref)
DEALLOCATE(refl)
DEALLOCATE(radv)
DEALLOCATE(rhohv)
!------------------------------------------------------------------------
!
! Process reflectivity data.
!
!------------------------------------------------------------------------
WRITE(stdout,'(a,i4)') &
' Processing base reflectivity data... ', icount
WRITE(stdout,*)'Transferring ',nazim,' reflectivity radials.'
dtime=float(itime-time1st)
refelvmin=min(elv,refelvmin)
refelvmax=max(elv,refelvmax)
DO j = 1, maxazim
time(j) = dtime
elev(j) = elv
END DO
WRITE(stdout,'(a,i5,a,4f9.1)') &
' Ref j azim elev refmin refmax'
DO j = 1, nazim
rmin=999.
rmax=-999.
DO i=1,ngate
IF(refdata(i,j) > refchek) THEN
rmin=min(refdata(i,j),rmin)
rmax=max(refdata(i,j),rmax)
END IF
END DO
IF(mod(j,5) == 0) &
WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax
END DO
gatesp = gtspc(1)
igatesp = NINT(gatesp)
irfirstg = NINT(rfirstg)
IF( qc_on) THEN
print *, ' calling SNR sensitivity screen'
CALL snrchk(maxgate,maxazim,ngate,nazim,refchek,rfirstg,gatesp, &
snstvty,refdata,veldata,arfdata)
print *, ' calling RhoHV screen'
CALL rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, &
refdata,veldata,rhodata)
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,refchek, &
refdata,rtem)
END IF
print *, ' calling volbuild reflectivity'
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,ngate,nazim, &
nyqset,timeset, &
kntrgat,kntrazm,kntrelv, &
igatesp,irfirstg,refchek, &
vnyquist,time, &
azim,elev,refdata, &
nyqvvol,timevol, &
rngrvol,azmrvol,elvrvol,refvol)
!------------------------------------------------------------------------
!
! Process velocity data.
!
!------------------------------------------------------------------------
icount = icount + 1
WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp
WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') &
' DATE: ', iyear, '/', imonth, '/', iday
WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') &
' TIME: ', ihr, ':', imin, ':', isec
WRITE(stdout,'(a,i4)')' Processing base velocity data... ', icount
WRITE(stdout,'(a,i6,a)')'Transferring',nazim,' velocity radials.'
WRITE(stdout,'(a,f10.2)') ' Nyquist velocity: ',vnyquist
dtime=float(itime-time1st)
DO j = 1, maxazim
time(j) = dtime
elev(j) = elv
END DO
WRITE(stdout,'(a,i5,a,4f9.1)') &
' Ref j azim elev vr min vr max'
DO j = 1, nazim
rmin=999.
rmax=-999.
DO i=1,ngate
IF(veldata(i,j) > velchek) THEN
rmin=min(veldata(i,j),rmin)
rmax=max(veldata(i,j),rmax)
END IF
END DO
IF(mod(j,5) == 0) &
WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax
END DO
gatesp = gtspc(1)
igatesp = NINT(gatesp)
irfirstg = NINT(rfirstg)
!------------------------------------------------------------------------
!
! Call radar volume builder.
!
!------------------------------------------------------------------------
ivelgatsp=igatesp
nyqset=1
timeset=1
IF( qc_on ) THEN
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,velchek, &
veldata,rtem)
print *,' calling unfnqc'
print *,' vnyquist=',vnyquist
CALL unfnqc(maxgate,maxazim,maxgate,nazim, &
nzsnd,zsnd,usnd,vsnd,rfrsnd, &
bkgopt,shropt,rfropt,igatesp,irfirstg, &
radlat,radlon,radelv, &
veldata,spwdata,elev,azim,vnyquist, &
unfvdata,bkgvel,bgate,egate,rtem)
print *, ' calling volbuild velocity '
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, &
nyqset,timeset, &
kntvgat,kntvazm,kntvelv, &
igatesp,irfirstg,velchek, &
vnyquist,time, &
azim,elev,unfvdata, &
nyqvvol,timevol, &
rngvvol,azmvvol,elvvvol,velvol)
ELSE
print *, ' calling volbuild velocity '
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, &
nyqset,timeset, &
kntvgat,kntvazm,kntvelv, &
igatesp,irfirstg,velchek, &
vnyquist,time, &
azim,elev,veldata, &
nyqvvol,timevol, &
rngvvol,azmvvol,elvvvol,velvol)
END IF
END DO ! infile
ELSE IF(kntt2b > 0) THEN
!
!------------------------------------------------------------------------
!
! Process reflectivity files
!
!------------------------------------------------------------------------
!
DO infile = 1,maxinfile
WRITE(6,'(a,i4,a,i4)') ' infile= ',infile, &
' type=',fvartype(infile)
IF(fvartype(infile) == 1) THEN ! type is reflectivity
infilname = ncradfn(infile)
WRITE(stdout,'(a,/a)') &
'Reading reflectivity file:',TRIM(infilname)
CALL extdirn(infilname,indir,fname)
CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, &
istatus)
ALLOCATE(ncvarname(nvar))
CALL get_ncrad2binfo(fname,indir,nazim,nvar, &
radarname,radlat,radlon,radelv,itimcdf,frtime, &
ivcp,elv,ncvarname,istatus)
DEALLOCATE(ncvarname)
ALLOCATE(refl(ngate,nazim))
CALL rdrftiltcdf(nazim,ngate,fname,indir, &
rmisval,rngfval,itimcdf,frtime,initime,rfirstg, &
azim,beamw,gtspc,refl)
bmwidth=beamw(1)
refdata=rmisval
DO j=1,nazim
DO i=21,ngate
refdata(i,j)=refl(i,j)
END DO
END DO
DEALLOCATE(refl)
IF(radarname(1:5) == 'cyril') THEN
IF(radelv < 0.1) radelv=445.3
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'chick') THEN
IF(radelv < 0.1) radelv=355.4
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'lawto') THEN
IF(radelv < 0.1) radelv=295.0
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'rushs') THEN
IF(radelv < 0.1) radelv=365.8
IF(elv < 0.1) elv=1.5
END IF
!------------------------------------------------------------------------
!
! Now begin processing the data.
!
!------------------------------------------------------------------------
icount = icount + 1
WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp
WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') &
' DATE: ', iyear, '/', imonth, '/', iday
WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') &
' TIME: ', ihr, ':', imin, ':', isec
!------------------------------------------------------------------------
!
! Process reflectivity data.
!
!------------------------------------------------------------------------
WRITE(stdout,'(a,i4)') &
' Processing base reflectivity data... ', icount
WRITE(stdout,*)'Transferring ',nazim,' reflectivity radials.'
dtime=float(itime-time1st)
refelvmin=min(elv,refelvmin)
refelvmax=max(elv,refelvmax)
DO j = 1, maxazim
time(j) = dtime
elev(j) = elv
END DO
WRITE(stdout,'(a,i5,a,4f9.1)') &
' Ref j azim elev refmin refmax'
DO j = 1, nazim, 5
rmin=999.
rmax=-999.
DO i=1,ngate
IF(refdata(i,j) > refchek) THEN
rmin=min(refdata(i,j),rmin)
rmax=max(refdata(i,j),rmax)
END IF
END DO
WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax
END DO
gatesp = gtspc(1)
igatesp = NINT(gatesp)
rfirstg = NINT(rfirstg)
!------------------------------------------------------------------------
!
! Call radar volume builder.
!
!------------------------------------------------------------------------
irefgatsp=igatesp
print *, ' gatesp,rfirstg,refchek=',gatesp,rfirstg,refchek
!
!------------------------------------------------------------------------
!
! Get matching rho-HV data
!
!------------------------------------------------------------------------
!
rho_match=.FALSE.
DO jfile=1,maxinfile
IF(fvartype(jfile) == 3 .AND. fitime(jfile) == itimcdf .AND. &
abs(felv(jfile)-elv) < 0.01 ) THEN
rho_match=.TRUE.
infilname = ncradfn(jfile)
WRITE(stdout,'(a,a)') 'Reading Rho-HV file:',TRIM(infilname)
CALL extdirn(infilname,indir,fname)
ALLOCATE(rhohv(ngate,nazim))
CALL rdrhotiltcdf(nazim,ngate,fname,indir,radarname, &
radlat,radlon,radelv,elv, &
rmisval,rngfval,itimcdf,frtime,initime, &
rfirstg,azim,rhohv)
rhodata=rmisval
print *, ' ngate =',ngate
print *, ' max range =',(ngate*0.024)
DO j=1,nazim
DO i=21,ngate
rhodata(i,j)=rhohv(i,j)
END DO
END DO
DEALLOCATE(rhohv)
EXIT
END IF
END DO ! jfile loop
IF(radarname(1:5) == 'cyril') THEN
IF(radelv < 0.1) radelv=445.3
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'chick') THEN
IF(radelv < 0.1) radelv=355.4
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'lawto') THEN
IF(radelv < 0.1) radelv=295.0
IF(elv < 0.1) elv=1.5
ELSE IF(radarname(1:5) == 'rushs') THEN
IF(radelv < 0.1) radelv=365.8
IF(elv < 0.1) elv=1.5
END IF
IF( qc_on) THEN
IF(rho_match) THEN
print *, ' calling RhoHV screen'
CALL rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, &
refdata,veldata,rhodata)
END IF
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,refchek, &
refdata,rtem)
END IF
nyqset=0
IF ( velproc ) THEN
timeset=0
ELSE
timeset=1
END IF
vnyquist=0.
print *, ' calling volbuild reflectivity'
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,ngate,nazim, &
nyqset,timeset, &
kntrgat,kntrazm,kntrelv, &
igatesp,irfirstg,refchek, &
vnyquist,time, &
azim,elev,refdata, &
nyqvvol,timevol, &
rngrvol,azmrvol,elvrvol,refvol)
END IF ! type is reflectivity
END DO ! file loop
!------------------------------------------------------------------------
!
! Process radial velocity files
!
!------------------------------------------------------------------------S
!
DO infile = 1,maxinfile
IF(fvartype(infile) == 2) THEN ! type is velocity
infilname = ncradfn(infile)
WRITE(stdout,'(a,a)') 'Reading velocity file:',TRIM(infilname)
CALL extdirn(infilname,indir,fname)
CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, &
istatus)
ALLOCATE(ncvarname(nvar))
CALL get_ncrad2binfo(fname,indir,nazim,nvar, &
radarname,radlat,radlon,radelv,itimcdf,frtime, &
ivcp,elv,ncvarname,istatus)
DEALLOCATE(ncvarname)
ALLOCATE(radv(ngate,nazim))
CALL rdvrtiltcdf(nazim,ngate,fname,indir, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyquist,rfirstg, &
azim,beamw,gtspc,vnyq,radv)
bmwidth=beamw(1)
veldata=-999.
DO j=1,nazim
DO i=1,ngate
veldata(i,j)=radv(i,j)
END DO
END DO
DEALLOCATE(radv)
!------------------------------------------------------------------------
!
! Correction for early experimental CASA data
!
!------------------------------------------------------------------------
IF(radarname(1:5) == 'cyril' .OR. &
radarname(1:5) == 'chick' .OR. &
radarname(1:5) == 'lawto' .OR. &
radarname(1:5) == 'rushs' ) THEN
vconst=15./acos(-1.)
DO j=1,nazim
DO i=1,ngate
IF(veldata(i,j) > -4.0) THEN
veldata(i,j)=veldata(i,j)*vconst
END IF
END DO
END DO
vnyquist=15.0
vnyq=15.0
END IF
!------------------------------------------------------------------------
!
! Process velocity data.
!
!------------------------------------------------------------------------
icount = icount + 1
WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp
WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') &
' DATE: ', iyear, '/', imonth, '/', iday
WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') &
' TIME: ', ihr, ':', imin, ':', isec
WRITE(stdout,'(a,i4)')' Processing base velocity data... ', icount
WRITE(stdout,'(a,i6,a)')'Transferring',nazim,' velocity radials.'
WRITE(stdout,'(a,f10.2)') ' Nyquist velocity: ',vnyquist
dtime=float(itime-time1st)
DO j = 1, maxazim
time(j) = dtime
elev(j) = elv
END DO
DO j = 1, nazim, 20
WRITE(stdout,'(a,i5,a,f8.1)') &
' Vel j =',j,' azim =',azim(j)
END DO
gatesp = gtspc(1)
igatesp = NINT(gtspc(1))
irfirstg = NINT(rfirstg)
!------------------------------------------------------------------------
!
! Call radar volume builder.
!
!------------------------------------------------------------------------
ivelgatsp=igatesp
nyqset=1
timeset=1
IF( qc_on ) THEN
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,velchek, &
veldata,rtem)
print *, ' calling unfnqc'
print *, ' Nyquist velocity: ',vnyquist
CALL unfnqc(maxgate,maxazim,maxgate,nazim, &
nzsnd,zsnd,usnd,vsnd,rfrsnd, &
bkgopt,shropt,rfropt,igatesp,irfirstg, &
radlat,radlon,radelv, &
veldata,spwdata,elev,azim,vnyquist, &
unfvdata,bkgvel,bgate,egate,rtem)
print *, ' calling volbuild velocity '
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, &
nyqset,timeset, &
kntvgat,kntvazm,kntvelv, &
igatesp,irfirstg,velchek, &
vnyquist,time, &
azim,elev,unfvdata, &
nyqvvol,timevol, &
rngvvol,azmvvol,elvvvol,velvol)
ELSE
print *, ' calling volbuild velocity '
print *, ' elv=',elev(1),' nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, &
nyqset,timeset, &
kntvgat,kntvazm,kntvelv, &
igatesp,irfirstg,velchek, &
vnyquist,time, &
azim,elev,veldata, &
nyqvvol,timevol, &
rngvvol,azmvvol,elvvvol,velvol)
END IF
END IF ! velocity
END DO ! DO ifile = 1,maxinfile
END IF ! knt2b > 0
!
! Restore any gzipped files
!
DO infile=1,maxinfile
IF(gzipped(infile)) THEN
WRITE(syscall,'(a,1x,a)') 'gzip',TRIM(ncradfn(infile))
CALL system(syscall)
END IF
END DO
END IF ! IF (maxinfile > 0)
!------------------------------------------------------------------------
!
! Create filename for output of remapped reflectivity and velocity.
!
!------------------------------------------------------------------------
IF (icount > 0) THEN
!------------------------------------------------------------------------
!
! Call remapping routines
!
!------------------------------------------------------------------------
IF( qc_on ) THEN
print *, ' Calling apdetect '
CALL apdetect(maxgate,maxgate,maxazim,maxelev, &
kntrgat,kntrazm,kntrelv, &
kntvgat,kntvazm,kntvelv, &
refchek,velchek, &
irefgatsp,ivelgatsp, &
winszrad,winszazim,ivcp, &
rngrvol,azmrvol,elvrvol, &
rngvvol,azmvvol,elvvvol, &
refvol,velvol,rtem, &
istatus)
END IF
nyqset=0
IF( velproc ) THEN
timeset=0
ELSE
timeset=1
END IF
IF(ref3d) THEN
vardump=1
ELSE
vardump=0
END IF
varfill=0
ivar=1
print *, ' Calling remapvol for reflectivity '
CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz,nzsnd,maxsort, &
vardump,varfill,ivar,nyqset,timeset,rfropt, &
refid,refname,refunits, &
refchek,refmiss,bmwidth,refmedl,refdazl,iordref, &
kntrgat,kntrazm,kntrelv, &
radlat,radlon,radarx,radary,radelv,dazim, &
rngmin,rngmax,time1st, &
rngrvol,azmrvol,elvrvol, &
refvol,timevol,nyqvvol,rxvol,ryvol,rzvol, &
xs,ys,zps,zsnd,rfrsnd,varsort,elvmean, &
gridref,gridtim,gridnyq,istatus)
IF( ref2d ) THEN
vardump=1
ivar=1
varid='refl2d'
varname='Low-level reflect'
varunits='dBZ'
print *, ' Calling remap2d for reflectivity'
CALL remap2d(maxgate,maxazim,maxelev,nx,ny,nzsnd,maxsort, &
vardump,ivar,rfropt,varid,varname,varunits,dmpfmt,hdf4cmpr, &
refchek,refmiss,refmedl,refdazl,iordref, &
kntrgat,kntrazm,kntrelv, &
radlat,radlon,radarx,radary,radelv,dazim, &
rngmin,rngmax,rngrvol,azmrvol,elvrvol, &
refvol,rxvol,ryvol,xs,ys,zsnd,rfrsnd,varsort, &
tem2d,istatus)
END IF
IF( velproc ) THEN
nyqset=1
timeset=1
IF(vel3d) THEN
vardump=1
ELSE
vardump=0
END IF
vardump=0
varfill=0
ivar=2
WRITE(stdout,'(a)') ' Calling remapvol for velocity '
CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz,nzsnd,maxsort, &
vardump,varfill,ivar,nyqset,timeset,rfropt, &
velid,velname,velunits, &
velchek,velmiss,bmwidth,velmedl,veldazl,iordvel, &
kntvgat,kntvazm,kntvelv, &
radlat,radlon,radarx,radary,radelv,dazim, &
rngmin,rngmax,time1st, &
rngvvol,azmvvol,elvvvol, &
velvol,timevol,nyqvvol,rxvol,ryvol,rzvol, &
xs,ys,zps,zsnd,rfrsnd,varsort,elvmean, &
gridvel,gridtim,gridnyq,istatus)
END IF
IF( vel2d ) THEN
vardump=1
ivar=1
varid='radv2d'
varname='Low-level Velocity'
varunits='m/s'
print *, ' Calling remap2d for velocity '
CALL remap2d(maxgate,maxazim,maxelev,nx,ny,nzsnd,maxsort, &
vardump,ivar,rfropt,varid,varname,varunits,dmpfmt,hdf4cmpr, &
velchek,velmiss,velmedl,veldazl,iordvel, &
kntvgat,kntvazm,kntvelv, &
radlat,radlon,radarx,radary,radelv,dazim, &
rngmin,rngmax,rngvvol,azmvvol,elvvvol, &
velvol,rxvol,ryvol,xs,ys,zsnd,rfrsnd,varsort, &
tem2d,istatus)
END IF
!------------------------------------------------------------------------
!
! Create filename and write file.
!
!------------------------------------------------------------------------
full_fname=' '
CALL mkradfnm(dmpfmt,dir_name,len_dir,radname,iiyr,iimonth,iiday, &
iihr, iimin, iisec, full_fname, len_fname)
WRITE(stdout,'(a,a)') ' Filename for this volume: ',TRIM(full_fname)
! CALL wtradcol(nx, ny, nz, dmpfmt, 1, hdf4cmpr, full_fname, &
! radname, radlat, radlon, radelv, &
! iyr,imonth,iday,iihr,iimin,iisec,ivcp,isource, &
! refelvmin,refelvmax,refrngmin,refrngmax, &
! xs, ys, zps, gridvel, gridref, gridnyq, gridtim, &
! tem1d1,tem1d2,tem1d3,tem1d4,tem1d5,tem1d6)
CALL wtradcol(nx,ny,nz,dmpfmt,1,hdf4cmpr,full_fname, &
radname,radlat,radlon,radelv, &
iyr,imonth,iday,iihr,iimin,iisec,ivcp,isource, &
refelvmin,refelvmax,refrngmin,refrngmax, &
xs,ys,zps,gridvel,gridref,gridnyq,gridtim,tem2d);
END IF
!------------------------------------------------------------------------
!
! The End.
!
!------------------------------------------------------------------------
WRITE(stdout,'(/a)') ' Normal termination of ncrad2arps.'
STOP
END PROGRAM ncrad2arps
SUBROUTINE rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, & 2
refl,radv,rhohv)
IMPLICIT NONE
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
REAL, INTENT(IN) :: refchek
REAL, INTENT(IN) :: rhohvlim
REAL, INTENT(INOUT) :: refl(maxgate,maxazim)
REAL, INTENT(INOUT) :: radv(maxgate,maxazim)
REAL, INTENT(IN) :: rhohv(maxgate,maxazim)
!
REAL, PARAMETER :: refmiss=-901.0
REAL, PARAMETER :: velmiss=-901.0
INTEGER, PARAMETER :: nhist=20
INTEGER :: histknt(nhist)
INTEGER :: igate,jazim,ihist
INTEGER :: kntvalid,kntflag
REAL :: histcst,histinv,histbin,pctflg,pcthist
histknt=0
histcst=1.0/float(nhist)
histinv=float(nhist)
kntvalid=0
kntflag=0
DO jazim=1,nazim
DO igate=1,ngate
IF(refl(igate,jazim) > refchek ) THEN
kntvalid=kntvalid+1
IF( rhohv(igate,jazim) < rhohvlim) THEN
kntflag=kntflag+1
refl(igate,jazim) = refmiss
radv(igate,jazim) = velmiss
END IF
ihist=1+INT(rhohv(igate,jazim)*histinv)
ihist=min(max(ihist,1),nhist)
histknt(ihist)=histknt(ihist)+1
END IF
END DO
END DO
!
pctflg=0.
IF(kntvalid > 0) THEN
pctflg=100.*float(kntflag)/float(kntvalid)
END IF
WRITE(6,'(a,i9,a,i9,a,f5.1,a)') &
' Rho-HV screened ',kntflag,' of ',kntvalid,' = ',pctflg, &
'% reflectivity gates'
!
WRITE(6,'(/a)') ' Histogram of Rho-HV (count, percent):'
DO ihist=1,nhist
histbin=histcst*(ihist-1)
pcthist=100.*float(histknt(ihist))/float(kntvalid)
WRITE(6,'(f7.3,i9,f10.1)') histbin,histknt(ihist),pcthist
END DO
RETURN
END SUBROUTINE rhohvchk
SUBROUTINE snrchk(maxgate,maxazim,ngate,nazim,refchek,rfirstg,gatesp, & 1,1
snstvty,refl,radv,attref)
IMPLICIT NONE
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: ngate
INTEGER, INTENT(IN) :: nazim
REAL, INTENT(IN) :: refchek
REAL, INTENT(IN) :: rfirstg
REAL, INTENT(IN) :: gatesp
REAL, INTENT(OUT) :: snstvty(maxgate)
REAL, INTENT(INOUT) :: refl(maxgate,maxazim)
REAL, INTENT(INOUT) :: radv(maxgate,maxazim)
REAL, INTENT(IN) :: attref(maxgate,maxazim)
!
REAL, PARAMETER :: refmiss=-911.0
REAL, PARAMETER :: velmiss=-911.0
INTEGER :: igate,jazim
INTEGER :: kntvalid,kntflag
REAL :: snr,pctflg
CALL xbsnstvty(maxgate,rfirstg,gatesp,snstvty)
DO igate=1,maxgate,10
print *, ' igate: ',igate,' range=',(0.001*(rfirstg+(igate-1)*gatesp)), &
' sens= ',snstvty(igate)
END DO
kntvalid=0
kntflag=0
DO jazim=1,nazim
DO igate=1,ngate
IF(refl(igate,jazim) > refchek ) THEN
kntvalid=kntvalid+1
IF( attref(igate,jazim) < snstvty(igate)) THEN
kntflag=kntflag+1
refl(igate,jazim) = refmiss
radv(igate,jazim) = velmiss
END IF
END IF
END DO
END DO
!
pctflg=0.
IF(kntvalid > 0) THEN
pctflg=100.*float(kntflag)/float(kntvalid)
END IF
WRITE(6,'(a,i9,a,i9,a,f5.1,a)') &
' SNR screened ',kntflag,' of ',kntvalid,' = ',pctflg, &
'% reflectivity gates'
RETURN
END SUBROUTINE snrchk
SUBROUTINE xbsnstvty(ngate,rfirstg,gatesp,sens) 1
!
!-----------------------------------------------------------------------
! Calculate X-band radar sensitivity (dBZ) as a function of range.
!
! Based on information provided by Francesc Joyvount
!
! Keith Brewster and Erin Fay, CAPS
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: ngate
REAL, INTENT(IN) :: rfirstg
REAL, INTENT(IN) :: gatesp
REAL, INTENT (OUT):: sens(ngate) ! Sensitivity values for all dist from radar.
!
! Parameters for CASA radar
!
! REAL, PARAMETER :: lmbda = 3.0E-02 ! Wavelength (m)
REAL, PARAMETER :: freqkhz=11505000.0 ! kHz
REAL, PARAMETER :: ptxdbm=66.68 ! Peak power (dBm)
! REAL, PARAMETER :: pt = 12500 ! Peak power (Watts)
REAL, PARAMETER :: G = 38.0 ! Antenna gain (dB)
REAL, PARAMETER :: F = 5.5 ! Noise figure (dB)
! REAL, PARAMETER :: tau = 0.67E-06 ! Radar pulse length (s)
REAL, PARAMETER :: tau = 4.06E-07 ! Radar pulse length (s)
REAL, PARAMETER :: theta = 1.8 ! Antenna half-power beamwidth (deg)
REAL, PARAMETER :: B = 2.5 ! Bandwidth (MHz)
REAL, PARAMETER :: lm = 1.0 ! Receiver mis-match loss(dB)
REAL, PARAMETER :: Kw2 = 0.91 ! Refractive Index of water squarred
REAL, PARAMETER :: T0 = 300. ! Rx temperature (K)
REAL, PARAMETER :: Ta = 200. ! Antenna temperature (K)
!
! Physical constants
!
REAL, PARAMETER :: K = 1.38E-23 ! Boltsmann's Constant (J/K)
REAL, PARAMETER :: c = 2.99792E08 ! Speed of light (m/s)
REAL, PARAMETER :: rconst =1.0E18 ! m6/m3 to mm6/m3
!
! Misc internal variables
!
INTEGER :: igate
REAL :: ln2,pi,pi3,four3,bwrad,rnoise,BHz,Ni,sconst,rlm,rG
REAL :: pt,lmbda
REAL :: range
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
pt=0.001*(10.0**(0.1*ptxdbm))
print *, ' pt= ',pt
lmbda=c/(freqkhz*1000.)
print *, ' lmbda (m)= ',lmbda
ln2=log(2.0)
pi=acos(-1.)
pi3=pi*pi*pi
four3=4.*4.*4.
bwrad=theta*pi/180.
rnoise=10.**(0.1*F)
rlm=10.**(0.1*lm)
rG=10.**(0.1*G)
BHz=1.0E06*b
Ni=K*(Ta+(rnoise-1.0)*T0)*BHz
sconst=rconst *(Ni/(pi3*Kw2)) * ((8.*ln2)/(bwrad*bwrad)) * &
(2./(c*tau)) * ((lmbda*lmbda*four3*rlm)/(Pt*rG*rG))
DO igate = 1, ngate
range=rfirstg+((igate-1)*gatesp)
sens(igate) = 10.*LOG10(range*range*sconst)
END DO
END SUBROUTINE xbsnstvty
SUBROUTINE extdirn(filename,dirname,fname) 5
!
! Extract directory name from complete filename and output
! directory name and filename in two character string variables
!
IMPLICIT NONE
CHARACTER (LEN=*) :: filename
CHARACTER (LEN=*) :: dirname
CHARACTER (LEN=*) :: fname
INTEGER :: lenf
INTEGER :: i
lenf=LEN(filename)
dirname='./'
fname=filename
DO i=(lenf-1),1,-1
IF(filename(i:i) == '/') THEN
dirname=filename(1:(i-1))
fname=filename((i+1):lenf)
EXIT
END IF
END DO
RETURN
END SUBROUTINE extdirn