!########################################################################
!########################################################################
!######### #########
!######### PROGRAM nids2arps #########
!######### #########
!######### Developed by #########
!######### Center for Analysis and Prediction of Storms #########
!######### University of Oklahoma #########
!######### #########
!########################################################################
!########################################################################
PROGRAM nids2arps,31
!------------------------------------------------------------------------
!
! PURPOSE:
!
! Reads WSR 88D radar data from NIDS raw data files (e.g., provided
! by WSI) and passes the data to ARPS remapping routines. To run,
! use:
! nids2arps < nids2arps.input
!
!------------------------------------------------------------------------
!
! AUTHOR:
!
! MODIFICATIONS:
!
! Yunheng Wang (July 2001)
! Converted from nids2arps.c and using NWS public code instead of
! wsibase.c. Based on the similar work on Echo Top and VIL by
! Eric Kemp 2001
!
! Eric Kemp, 9 August 2001
! Added VIL, Echo Top, and DPA product remapping -- these data are
! currently written in ARPSPLT "arbitrary variable" format to allow for
! plotting. In the future, a better file format will be developed and
! added for used with ADAS. Also, changed code to use INTEGER*4 and
! similar variable types instead of using the KIND statement. Finally,
! consolidated parameters into module N2ACONST, and cleaned up code.
!
! Eric Kemp, 17 August 2001
! Modified code to avoid use of INTEGER*4, INTEGER*1, etc. Based on
! work by Yunheng Wang.
!
! Eric Kemp, 10 September 2001
! Modified driver code to avoid use of unit 11. This fixes conflict
! with subroutine 'setlookup', which uses unit 11 for reading/writing the
! radar map parameter file.
!
! Eric Kemp, 17 September 2001
! Consolidated code, allowing for all NIDS data files to be listed
! in the nids2arps.input file in the new 'nidsfn' variable array. Also,
! added an option to adjust the remapped 3D reflectivity field using
! the NIDS Echo Top and VIL products.
!
!------------------------------------------------------------------------
!
! REMARKS:
!
! DIFFERENCE from the previous C code
!
! 1. Radar lat, lon, alt decoded from NIDS data instead of inputing
! from external file (i.e. radarinfo.dat)
! 2. itime is the actual time recorded in the NIDS data, while the
! C code use only time in the first file for every input file.
!
! The above two difference will make the output file a little
! different from that got by running the legacy (nids2arps.c).
!
!------------------------------------------------------------------------
!
! Use modules.
!
!------------------------------------------------------------------------
USE n2aconst
!------------------------------------------------------------------------
!
! Force explicit declarations.
!
!------------------------------------------------------------------------
IMPLICIT NONE
!------------------------------------------------------------------------
!
! Include files
!
!------------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
!------------------------------------------------------------------------
!
! Variable Declarations.
!
!------------------------------------------------------------------------
INTEGER :: nx, ny, nz, nzsoil, itime, nstyps
REAL :: radar_alt, radar_lat, radar_lon
REAL :: radarx, radary
!------------------------------------------------------------------------
!
! NWS decoder variables
!
!------------------------------------------------------------------------
INTEGER :: prodID ! product ID code
INTEGER :: prodtype ! product type code
INTEGER :: ival0, ival1
INTEGER :: ival(maxlen)
INTEGER :: iheader(iheaderdim)
INTEGER :: ifield(maxgate,maxazim)
INTEGER :: ifield_raster(ndim1,ndim2)
INTEGER :: ifield_dpa(ndim1,ndim2)
CHARACTER(LEN=1) :: itrailer(itrailerdim)
INTEGER :: msglen,icode,isite,isite_lat,isite_lon,isite_elev, &
iyr,imonth,iday,ihr,iminit,isec,iprod,ngate,nazim, &
maxval2,ihr1,iminit1,ihr2,iminit2,istm_spd,istm_dir, &
len_header,len_trailer
INTEGER :: iazmuth(maxazim)
INTEGER :: maxval1, ivcp, ielev
INTEGER :: iyr1,imonth1,iday1,iyr2,imonth2,iday2
INTEGER :: icatt(0:15)
INTEGER :: icats(0:15)
REAL :: rcats(0:15)
INTEGER :: ndims_p,nrowsp,ncolsp
INTEGER :: iyr_save
!------------------------------------------------------------------------
!
! Variables for subroutine remaptilt.
!
!------------------------------------------------------------------------
REAL, ALLOCATABLE :: rdata(:,:) ! refl (dBZ) or velocity data (m/s)
REAL, ALLOCATABLE :: rtem(:,:) ! temporary array for radar data processing
REAL, ALLOCATABLE :: time(:) ! time offset from time1st
REAL, ALLOCATABLE :: azim(:) ! azimuth angle for each radial(degree)
REAL, ALLOCATABLE :: elev(:) ! elevation angle for each radial (degree)
INTEGER :: time1st ! time of first radial in volume
INTEGER :: rfirstg ! range to first gate (meters)
INTEGER :: gatesp ! gate spacing (meters)
REAL :: eleva ! elevation angle for this tilt
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 :: 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 :: zpsoil(:,:,:)
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 :: j3soil(:,:,:)
REAL, ALLOCATABLE :: j3soilinv(:,:,:)
REAL, ALLOCATABLE :: tem1d1(:)
REAL, ALLOCATABLE :: tem1d2(:)
REAL, ALLOCATABLE :: tem1d3(:)
REAL, ALLOCATABLE :: tem1d4(:)
REAL, ALLOCATABLE :: tem1d5(:)
REAL, ALLOCATABLE :: tem1d6(:)
REAL, ALLOCATABLE :: tem3d(:,:,:)
REAL, ALLOCATABLE :: grdvel2(:,:,:)
REAL, ALLOCATABLE :: grdref2(:,:,:)
REAL, ALLOCATABLE :: wgtvel(:,:,:)
REAL, ALLOCATABLE :: wgtref(:,:,:)
REAL, ALLOCATABLE :: wpotvel(:,:,:)
REAL, ALLOCATABLE :: wpotref(:,:,:)
REAL, ALLOCATABLE :: gridazm(:,:)
REAL, ALLOCATABLE :: gridrng(:,:)
!------------------------------------------------------------------------
!
! Arrays for raster (Echo Top and VIL) and DPA data.
!
!------------------------------------------------------------------------
REAL, ALLOCATABLE :: rasterdata(:,:)
REAL, ALLOCATABLE :: rasterx(:,:)
REAL, ALLOCATABLE :: rastery(:,:)
REAL, ALLOCATABLE :: rasterlat(:,:)
REAL, ALLOCATABLE :: rasterlon(:,:)
REAL, ALLOCATABLE :: remapped_raster(:,:)
REAL, ALLOCATABLE :: dpadata(:,:)
REAL, ALLOCATABLE :: dpax(:,:)
REAL, ALLOCATABLE :: dpay(:,:)
REAL, ALLOCATABLE :: dpalat(:,:)
REAL, ALLOCATABLE :: dpalon(:,:)
REAL, ALLOCATABLE :: remapped_dpa(:,:)
REAL, ALLOCATABLE :: arpslat(:,:)
REAL, ALLOCATABLE :: arpslon(:,:)
REAL, ALLOCATABLE :: arpsrasterx(:,:)
REAL, ALLOCATABLE :: arpsrastery(:,:)
REAL, ALLOCATABLE :: arpsdpax(:,:)
REAL, ALLOCATABLE :: arpsdpay(:,:)
!------------------------------------------------------------------------
!
! New variables for adjusting reflectivity (6 September 2001)
!
!------------------------------------------------------------------------
REAL, ALLOCATABLE :: nidset(:,:)
REAL, ALLOCATABLE :: nidsvil(:,:)
!------------------------------------------------------------------------
!
! Misc. variables
!
!------------------------------------------------------------------------
CHARACTER(LEN=1) :: cval(maxlen)
CHARACTER(LEN=256) :: infile
INTEGER :: n, i, j, k, ios, istatus, istatus1, istatus2
INTEGER :: irefgatsp,ivelgatsp
INTEGER :: unitin, unitout, unitl, llbound, lubound
INTEGER :: dmpfmt,hdf4cmpr
INTEGER :: lenmpflstr
CHARACTER(LEN=100) full_fname
INTEGER :: len_dir, len_fname
LOGICAL :: velproc
LOGICAL :: i_first_scan
INTEGER :: iiyr,iimonth,iiday,iihr,iiminit,iisec
INTEGER :: isource, count
INTEGER :: nyqset,timeset,vardump
INTEGER :: xscale
REAL :: rdx,dtime
INTEGER :: remapopt
REAL :: radius
REAL :: xrad,yrad
REAL :: dBA
REAL :: dpatrulat(2)
REAL :: dpa_x_center, dpa_y_center
INTEGER :: raddattype
INTEGER :: maxinfile,infilenum
LOGICAL :: assigned_vil, assigned_et
INTEGER :: iarg,jarg,iargc,narg,ifile,kfile
CHARACTER(LEN=132) charg
CHARACTER(LEN=132) listfile
!------------------------------------------------------------------------
!
! Some tunable parameters - convert to Namelist in future
!
!------------------------------------------------------------------------
INTEGER, PARAMETER :: iordref = 2
INTEGER, PARAMETER :: iordvel = 2
REAL, PARAMETER :: refmedl = 20.
REAL, PARAMETER :: velmedl = 15.
REAL, PARAMETER :: refdazl = 360.
REAL, PARAMETER :: veldazl = 30.
REAL, PARAMETER :: dazim = 1.0
REAL, PARAMETER :: rngmin = 10.E03
REAL, PARAMETER :: rngmax = 230.E03
!------------------------------------------------------------------------
!
! Namelist declaration
!
!------------------------------------------------------------------------
CHARACTER(LEN=4) :: radar_name
INTEGER :: nnidsfn = 0
CHARACTER(LEN=256) :: radar_map_file
CHARACTER(LEN=256) :: nidsfn(maxnidsfile)
INTEGER :: map_flg
CHARACTER(LEN=120) dir_name
CHARACTER(LEN=256) :: etfn,vilfn,dpafn
INTEGER :: et_remapopt,vil_remapopt,dpa_remapopt
REAL :: et_radius,vil_radius,dpa_radius
INTEGER :: arbvaropt
INTEGER :: adjreflopt
NAMELIST /nids_data/ radar_name,nnidsfn,nidsfn, &
radar_map_file, map_flg, dir_name, &
et_remapopt, et_radius, &
vil_remapopt, vil_radius, &
dpa_remapopt,dpa_radius,arbvaropt,adjreflopt
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!------------------------------------------------------------------------
!
! Initializations
!
!------------------------------------------------------------------------
isource = 2 ! 1 - WSR-88D raw 2 - WSR-88D NIDS
dmpfmt=1
hdf4cmpr=0
unitin = 1
i_first_scan = .TRUE.
velproc = .TRUE.
!------------------------------------------------------------------------
!
! Initialize nids_data NAMELIST variables
!
!------------------------------------------------------------------------
radar_name = 'DMMY'
nnidsfn = 0
radar_map_file = 'dummy'
map_flg = 0
dir_name = './'
etfn = 'dummy'
vilfn = 'dummy'
dpafn = 'dummy'
et_radius = 10.0E03
vil_radius = 10.0E03
dpa_radius = 10.0E03
arbvaropt = 0
adjreflopt = 0
DO ifile=1,maxnidsfile
nidsfn(ifile)='dummy'
END DO
!------------------------------------------------------------------------
!
! Read nids2arps.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
!
!------------------------------------------------------------------------
narg=iargc()
WRITE(stdout,'(a,i4)') ' Number of command-line arguments: ',narg
IF(narg > 1 ) THEN
CALL getarg(1,charg)
radar_name=charg(1:4)
WRITE(stdout,'(a,a)') ' radar_name = ', radar_name
kfile=0
nnidsfn=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
listfile=charg
WRITE(stdout,'(a)') ' Reading file lists from command line'
OPEN(31,file=listfile,status='old',form='formatted')
DO ifile=(nnidsfn+1),maxnidsfile
READ(31,'(a)',END=101) nidsfn(ifile)
kfile=kfile+1
WRITE(stdout,'(a,i4,a,a)') &
' kfile:',kfile,' nidsfn: ',TRIM(nidsfn(ifile))
END DO
101 CONTINUE
nnidsfn=kfile
END IF
iarg=iarg+1
IF(iarg > narg) EXIT
END DO
WRITE(stdout,'(i4,a)') nnidsfn,' file names read'
END IF
!
! If file lists have not been obtained from command line
! then process the nids block in the input file
!
IF( nnidsfn == 0 ) THEN
WRITE(stdout,'(a)') ' Getting file lists from input file'
READ(5,nids_data,END=999)
WRITE(stdout,'(/a)') ' Namelist block nids_data successfully read.'
lenmpflstr = LEN(TRIM(radar_map_file))
len_dir = LEN(TRIM(dir_name))
WRITE(stdout,'(a,a)') ' radar_name= ', radar_name
WRITE(stdout,'(a,a)') ' remap directory= ',TRIM(dir_name)
IF(map_flg == 1) THEN
WRITE(stdout,'(a,a)')' Radar map file (input)= ',radar_map_file
ELSE IF(map_flg == 2) THEN
WRITE(stdout,'(a,a)')' Radar map file (output)= ',radar_map_file
ELSE
WRITE(stdout,'(a)') ' No radar map file specified.'
END IF
END IF
!------------------------------------------------------------------------
!
! Array allocations and variable initializations.
!
!------------------------------------------------------------------------
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(rxvol(maxgate,maxazim,maxelev))
ALLOCATE(ryvol(maxgate,maxazim,maxelev))
ALLOCATE(rzvol(maxgate,maxazim,maxelev))
ALLOCATE(gridvel(nx,ny,nz))
ALLOCATE(gridref(nx,ny,nz))
ALLOCATE(gridnyq(nx,ny,nz))
ALLOCATE(gridtim(nx,ny,nz))
ALLOCATE( ALLOCATE( ALLOCATE( ALLOCATE(zp(nx,ny,nz))
ALLOCATE(zpsoil(nx,ny,nzsoil))
ALLOCATE(xs(nx))
ALLOCATE(ys(ny))
ALLOCATE(zps(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(j3soil(nx,ny,nzsoil))
ALLOCATE(j3soilinv(nx,ny,nzsoil))
ALLOCATE(tem1d1(nz))
ALLOCATE(tem1d2(nz))
ALLOCATE(tem1d3(nz))
ALLOCATE(tem1d4(nz))
ALLOCATE(tem1d5(nz))
ALLOCATE(tem1d6(nz))
ALLOCATE(tem3d(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(rdata(maxgate,maxazim))
ALLOCATE(rtem(maxgate,maxazim))
ALLOCATE(time(maxazim))
ALLOCATE(azim(maxazim))
ALLOCATE(elev(maxazim))
IF (adjreflopt == 1) THEN
ALLOCATE(nidset(nx,ny))
ALLOCATE(nidsvil(nx,ny))
nidset(:,:) = -9999.
nidsvil(:,:) = -9999.
END IF
assigned_vil = .FALSE.
assigned_et = .FALSE.
istatus = 0
count = 0
rfirstg = 1000
gatesp = 1000
CALL inigrd(nx,ny,nz,nzsoil, &
x,y,z,zp,zpsoil,hterain,mapfct,j1,j2,j3,j3soil,j3soilinv,&
tem3d)
DEALLOCATE(mapfct)
DEALLOCATE(j1)
DEALLOCATE(j2)
DEALLOCATE(j3)
DEALLOCATE(j3soil)
DEALLOCATE(j3soilinv)
DEALLOCATE(zpsoil)
ALLOCATE(arpslat(nx,ny))
ALLOCATE(arpslon(nx,ny))
ALLOCATE(arpsrasterx(nx,ny))
ALLOCATE(arpsrastery(nx,ny))
ALLOCATE(arpsdpax(nx,ny))
ALLOCATE(arpsdpay(nx,ny))
!------------------------------------------------------------------------
!
! Loop through the NIDS files.
!
!------------------------------------------------------------------------
IF (nnidsfn > 0 .AND. nnidsfn <= maxnidsfile) THEN
maxinfile = nnidsfn
ELSE IF (nnidsfn > maxnidsfile) THEN
WRITE(stdout,'(a,i4,a,i4)') &
' WARNING: nnidsfn=', nnidsfn,' > maxnidsfile=',maxnidsfile
WRITE(stdout,'(a,i4,a)') &
' Will only read in ', maxnidsfile,' NIDS files.'
maxinfile = maxnidsfile
ELSE
maxinfile = 0
END IF
IF (maxinfile > 0) THEN
DO infilenum = 1,maxinfile
infile = nidsfn(infilenum)
nazim = 0
ielev = 0
!------------------------------------------------------------------------
!
! Read product into array ival() below.
!
!------------------------------------------------------------------------
ival(:) = 0 ! Initialize array
unitin = 12
OPEN (UNIT=unitin, FILE=TRIM(infile), STATUS='OLD', &
ACCESS='DIRECT',RECL=1, FORM='UNFORMATTED', ERR=998, &
IOSTAT=ios)
DO n = 1, maxlen
READ (UNIT=unitin, REC=n, ERR=10, IOSTAT=ios) cval(n)
IF (ios /= 0) EXIT
END DO
10 CONTINUE
CLOSE(UNIT=unitin)
msglen = n - 1
DO i = 1, msglen
ival(i) = ICHAR(cval(i))
END DO
WRITE(stdout,'(a,a)') ' Finished reading file ',TRIM(infile)
WRITE (stdout,'(a,i5)', ADVANCE='NO') ' system code: ',ios
WRITE (stdout,'(a,i8)') ' Message length: ',msglen
!------------------------------------------------------------------------
!
! Get product ID code.
!
!------------------------------------------------------------------------
ival0 = ival(1); ival1 = ival(2)
IF (ival1 < 0) ival1 = 256 + ival1
IF (ival0 < 0) ival0 = 256 + ival0
prodID = ival1 + (256*ival0)
WRITE(stdout,'(a,i5)',ADVANCE='NO') ' Product ID: ',prodID
!------------------------------------------------------------------------
!
! Get product type code.
!
!------------------------------------------------------------------------
ival0 = ival(137); ival1 = ival(138)
IF (ival1 < 0) ival1 = 256 + ival1
IF (ival0 < 0) ival0 = 256 + ival0
prodtype = ival1 + (256*ival0)
WRITE (stdout,FMT='(a,i10)') ' Type code: ', prodtype
!------------------------------------------------------------------------
!
! Call appropriate decoding subroutine based on product type
! (radial, raster, or DPA). Note that the hybrid scan reflectivity
! routine could be added here in the future.
!
!------------------------------------------------------------------------
IF (prodtype == ipradial) THEN
WRITE(stdout,'(a)')'About to call get_radial_data...'
CALL get_radial_data(ival,msglen,isite, &
isite_lat,isite_lon,isite_elev, &
iyr,imonth,iday,ihr,iminit,isec, &
iprod,maxval1,ivcp,ifield,maxgate,maxazim, &
iazmuth,ielev,ngate,nazim,icatt,icats, &
maxval2,iyr1,imonth1,iday1,ihr1,iminit1, &
iyr2,imonth2,iday2,ihr2,iminit2,istm_spd, &
istm_dir,iheader,len_header,itrailer, &
len_trailer,icode)
IF (icode == 2) THEN
WRITE(stdout,'(a)')' ERROR -- Not a radial graphic product.'
CYCLE
ELSE IF (icode == 3) THEN
WRITE(stdout,'(a)')' ERROR -- ncols/nrows too small for product.'
CYCLE
ELSE IF (icode == 4) THEN
WRITE(stdout,'(a)')' ERROR -- msglen too small to store product.'
CYCLE
END IF
ELSE IF (prodtype == ipraster1 .OR. prodtype == ipraster2) THEN
WRITE(stdout,'(a)')'About to call get_raster_data...'
CALL get_raster_data(ival,msglen,isite, &
isite_lat,isite_lon,isite_elev, &
iyr,imonth,iday,ihr,iminit,isec, &
iprod,maxval1,ivcp,ifield_raster,ndim1,ndim2,&
icatt,icats,ndims_p,iheader,len_header, &
itrailer,len_trailer,icode)
IF (icode == 2) THEN
WRITE(stdout,'(a)')'ERROR -- Not a raster graphic product.'
CYCLE
ELSE IF (icode == 3) THEN
WRITE(stdout,'(a)')'ERROR -- ncols/nrows too small for product.'
CYCLE
ELSE IF (icode == 4) THEN
WRITE(stdout,'(a)')'ERROR -- msglen too small to store product.'
CYCLE
END IF
ELSE IF (prodID == DPAprodID) THEN
WRITE(stdout,'(a)')'About to call get_dpa_data...'
CALL get_dpa_data(ival,msglen,isite, &
isite_lat,isite_lon,isite_elev, &
iyr,imonth,iday,ihr,iminit,isec, &
iprod,maxval1,ivcp,ifield_dpa,ndim1,ndim2, &
ncolsp,nrowsp,iheader,len_header,itrailer, &
len_trailer,icode)
IF (icode == 2) THEN
WRITE(stdout,'(a)')'ERROR -- Not a DPA product.'
CYCLE
ELSE IF (icode == 3) THEN
WRITE(stdout,'(a)') &
'ERROR -- ncols/nrows too small for product.'
CYCLE
ELSE IF (icode == 4) THEN
WRITE(stdout,'(a)') &
'ERROR -- msglen too small to store product.'
CYCLE
END IF
ELSE
WRITE(stdout,'(a,/a)') &
'ERROR: Current file does not contain', &
'NIDS data. Moving to next file...'
CYCLE
END IF
!------------------------------------------------------------------------
!
! Save radar lat/lon, elevation, and number of radials.
!
!------------------------------------------------------------------------
radar_alt = isite_elev*mpfoot
radar_lat = isite_lat*0.001
radar_lon = isite_lon*0.001
eleva = ielev*0.1
WRITE(stdout,'(a)') 'Radar information:'
WRITE(stdout,'(a,a,a,f8.1)') ' NAME ', radar_name, &
' ALTITUDE ', radar_alt
WRITE(stdout,'(a,f10.4,a,f10.4)') ' LATITUDE ', radar_lat, &
' LONGITUDE ', radar_lon
WRITE(stdout,'(a,i3,a,i4)') &
(' Threshold ', i, ' icats = ', icats(i), i=0,15)
DO i=0,15
rcats(i)=float(icats(i))
END DO
!------------------------------------------------------------------------
!
! Initialize 3-d radar data arrays
!
!------------------------------------------------------------------------
CALL ctim2abss(iyr, imonth, iday, ihr, iminit, isec, itime)
iyr = MOD(iyr, 100)
IF (i_first_scan) THEN
CALL radcoord(nx,ny,nz, &
x,y,z,zp,xs,ys,zps, &
radar_lat,radar_lon,radarx,radary)
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
iiminit = iminit
iisec = isec
time1st=itime
i_first_scan = .FALSE.
END IF
!------------------------------------------------------------------------
!
! Now begin processing the data.
!
!------------------------------------------------------------------------
IF (prodtype == ipradial) THEN
count = count + 1
WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp
WRITE(stdout,'(1x,a,i2.2,a1,i2.2,a1,i2.2)', ADVANCE='NO') &
' DATE: ', iyr, '/', imonth, '/', iday
WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') &
' TIME: ', ihr, ':', iminit, ':', isec
!------------------------------------------------------------------------
!
! Process base reflectivity data.
!
!------------------------------------------------------------------------
IF ( (prodID >= BREFprodID1) .AND. (prodID <= BREFprodID2) ) THEN
WRITE(stdout,'(a,i4)') &
' Processing base reflectivity data... ', count
WRITE(stdout,*)'Transferring ', nazim, &
' reflectivity radials.'
DO j = 1, maxazim
DO i = 1, maxgate
IF (ifield(i,j) > 0 .AND. ifield(i,j) < 16) THEN
rdata(i,j) = rcats(ifield(i,j))
ELSE
rdata(i,j) = r_missing
END IF
END DO
END DO
dtime=float(itime-time1st)
DO j = 1, maxazim
time(j) = dtime
azim(j) = iazmuth(j)*0.10
elev(j) = eleva
END DO
DO j = 1, nazim, 20
WRITE(stdout,'(a,i5,a,f8.1)') &
' Ref j =',j,' azim =',azim(j)
END DO
IF (iprod < 0 .OR. iprod > 109) THEN
gatesp = 0 ! iprod: product code number
rfirstg = 0 ! should be product->code
ELSE ! in struct Proddesc?
gatesp = INT(1000*res(iprod))
rfirstg = INT(1000*res(iprod))
END IF
!------------------------------------------------------------------------
!
! Call radar volume builder.
!
!------------------------------------------------------------------------
irefgatsp=gatesp
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,refchek, &
rdata,rtem)
nyqset=0
IF ( velproc ) THEN
timeset=0
ELSE
timeset=1
END IF
vnyquist=0.
print *, ' calling volbuild reflectivity - nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,ngate,nazim, &
nyqset,timeset, &
kntrgat,kntrazm,kntrelv, &
gatesp,rfirstg,refchek, &
vnyquist,time, &
azim,elev,rdata, &
nyqvvol,timevol, &
rngrvol,azmrvol,elvrvol,refvol)
!------------------------------------------------------------------------
!
! Process velocity data.
!
!------------------------------------------------------------------------
ELSE IF ( (prodID >= BVELprodID1) .AND. &
(prodID <= BVELprodID2) ) THEN
WRITE(stdout,'(a,i4)')' Processing base velocity data... ', count
WRITE(stdout,'(a,i6,a)')'Transferring',nazim,' velocity radials.'
DO i = 1,14
rcats(i)=kts2ms*rcats(i)
END DO
vnyquist=rcats(14)
WRITE(stdout,'(a,f10.2)') ' Nyquist velocity: ',vnyquist
DO j = 1, maxazim
DO i = 1, maxgate
IF (ifield(i,j) > 0 .AND. ifield(i,j) < 15) THEN
rdata(i,j) = rcats(ifield(i,j))
ELSE
rdata(i,j) = r_missing
END IF
END DO
END DO
dtime=float(itime-time1st)
DO j = 1, maxazim
time(j) = dtime
azim(j) = 0.1* float(iazmuth(j))
elev(j) = eleva
END DO
DO j = 1, nazim, 20
WRITE(stdout,'(a,i5,a,f8.1)') &
' Vel j =',j,' azim =',azim(j)
END DO
IF (iprod < 0 .OR. iprod > 109) THEN
gatesp = 0 ! iprod: product code number
rfirstg = 0 ! should be product->code
ELSE ! in struct Proddesc?
gatesp = INT(1000*res(iprod))
rfirstg = INT(1000*res(iprod))
END IF
!------------------------------------------------------------------------
!
! Call radar volume builder.
!
!------------------------------------------------------------------------
ivelgatsp=gatesp
print *, ' calling despekl '
CALL despekl(maxgate,maxazim,maxgate,maxazim,velchek, &
rdata,rtem)
nyqset=1
timeset=1
print *, ' calling volbuild velocity - nazim =',nazim
CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, &
nyqset,timeset, &
kntvgat,kntvazm,kntvelv, &
gatesp,rfirstg,velchek, &
vnyquist,time, &
azim,elev,rdata, &
nyqvvol,timevol, &
rngvvol,azmvvol,elvvvol,velvol)
!------------------------------------------------------------------------
!
! Other product ID code for radial data.
!
!------------------------------------------------------------------------
ELSE
WRITE(stdout,'(a,i6)') &
' ERROR -- Unknown NIDS product ID code: ', prodID
CYCLE
END IF
iyr_save = iyr
ELSE IF (prodtype == ipraster1 .OR. prodtype == ipraster2) THEN
!------------------------------------------------------------------------
!
! Save resolution of raster grid.
!
!------------------------------------------------------------------------
ival1 = ival(148)
IF (ival1 < 0) ival1 = 256 + ival1
ival0 = ival(147)
IF (ival0 < 0) ival0 = 256 + ival0
xscale = ival1 + (256*ival0)
rdx = REAL(xscale)*1000.
!------------------------------------------------------------------------
!
! Allocate raster arrays and save decoded data.
!
!------------------------------------------------------------------------
ALLOCATE(rasterdata(ndims_p,ndims_p),STAT=istatus)
ALLOCATE(rasterx(ndims_p,ndims_p),STAT=istatus)
ALLOCATE(rastery(ndims_p,ndims_p),STAT=istatus)
ALLOCATE(rasterlat(ndims_p,ndims_p),STAT=istatus)
ALLOCATE(rasterlon(ndims_p,ndims_p),STAT=istatus)
ALLOCATE(remapped_raster(nx,ny),STAT=istatus)
rasterdata(:,:) = rasterchek
remapped_raster(:,:) = rastermiss
IF (prodID == ETprodID) THEN ! Echo tops
WRITE(stdout,'(a)') ' Processing NIDS Echo Top data...'
remapopt = et_remapopt
radius = et_radius
DO j = 1,ndims_p ! Column
DO i = 1,ndims_p ! Row
! rasterdata(i,j)=icats(ifield_raster(i,ndims_p-j+1))*mpfoot
! rasterdata(i,j)=icats(ifield_raster(i,ndims_p-j+1)) ! Test
rasterdata(i,j)=icats(ifield_raster(i,ndims_p-j+1))* &
mpfoot*1000.
END DO ! Row
END DO ! Column
ELSE IF (prodID == VILprodID) THEN ! VIL
WRITE(stdout,'(a)') ' Processing NIDS VIL data...'
remapopt = vil_remapopt
radius = vil_radius
DO j = 1,ndims_p ! Column
DO i = 1,ndims_p ! Row
rasterdata(i,j) = icats(ifield_raster(i,ndims_p-j+1))
END DO ! Row
END DO ! Column
ELSE
WRITE(stdout,'(a)') ' ERROR -- Unknown raster data type.'
WRITE(stdout,'(a)') ' Moving to next file.'
CYCLE
END IF
!------------------------------------------------------------------------
!
! Calculate Cartesian coordinates of raster data on raster grid
! relative to the radar, and set rasterdata to missing if it is
! beyond radar range.
!
!------------------------------------------------------------------------
CALL setrastergrid(ndims_p,radar_lat,radar_lon,rdx, &
rasterlat,rasterlon,rasterx,rastery, &
rasterdata,rasterchek)
xrad = 0. ! Radar at center of raster grid
yrad = 0. ! Radar at center of raster grid
!------------------------------------------------------------------------
!
! Determine Cartesian coordinates of ARPS scalar points on radar
! grid relative to radar.
!
!------------------------------------------------------------------------
CALL xytoll(nx,ny,xs,ys,arpslat,arpslon)
CALL getarpsrasterxy(nx,ny,arpslat,arpslon,radar_lat,radar_lon, &
arpsrasterx,arpsrastery)
!------------------------------------------------------------------------
!
! Remap the data and write to file.
!
!------------------------------------------------------------------------
CALL remapraster(nx,ny,arpsrasterx,arpsrastery,dx,rasterchek, &
rastermiss,xrad,yrad, &
remapopt,radius,ndims_p,ndims_p,rasterdata, &
rasterx,rastery,remapped_raster,istatus)
IF (prodID == ETprodID) THEN ! Echo tops
IF (arbvaropt == 1) THEN
CALL wrtvar1(nx,ny,1,remapped_raster,radar_name//'et', &
'NIDS Echo Tops','m',000000,runname(1:lfnkey), &
TRIM(dir_name),istatus)
END IF
IF (adjreflopt == 1 .AND. .NOT. assigned_et) THEN
nidset(:,:) = remapped_raster(:,:)
assigned_et = .TRUE.
ELSE IF (adjreflopt == 1 .AND. assigned_et) THEN
WRITE(stdout,'(a)') &
' WARNING: NIDS Echo Top product already remapped!'
END IF
ELSE IF (prodID == VILprodID) THEN ! VIL
IF (arbvaropt == 1) THEN
CALL wrtvar1(nx,ny,1,remapped_raster,radar_name//'vl', &
'NIDS VIL','kg m^-2',000000,runname(1:lfnkey), &
TRIM(dir_name),istatus)
END IF
IF (adjreflopt == 1 .AND. .NOT. assigned_vil) THEN
nidsvil(:,:) = remapped_raster(:,:)
assigned_vil = .TRUE.
ELSE IF (adjreflopt == 1 .AND. assigned_vil) THEN
WRITE(stdout,'(a)') &
' WARNING: NIDS VIL product already remapped!'
END IF
END IF
DEALLOCATE(rasterdata,rasterx,rastery,rasterlat,rasterlon, &
remapped_raster,STAT=istatus)
ELSE IF (prodID == DPAprodID) THEN
!------------------------------------------------------------------------
!
! Allocate arrays and save decoded data.
!
!------------------------------------------------------------------------
WRITE(stdout,'(a)') ' Processing NIDS DPA data...'
ALLOCATE(dpadata(nrowsp,ncolsp),STAT=istatus)
ALLOCATE(dpax(nrowsp,ncolsp),STAT=istatus)
ALLOCATE(dpay(nrowsp,ncolsp),STAT=istatus)
ALLOCATE(dpalat(nrowsp,ncolsp),STAT=istatus)
ALLOCATE(dpalon(nrowsp,ncolsp),STAT=istatus)
ALLOCATE(remapped_dpa(nx,ny),STAT=istatus)
dpadata(:,:) = dpachek
remapped_dpa(:,:) = dpamiss
remapopt = dpa_remapopt
radius = dpa_radius
DO j = 1,ncolsp
DO i = 1,nrowsp
IF (ifield_dpa(i,ncolsp-j+1) == 0) THEN
dpadata(i,j) = 0.
ELSE IF (ifield_dpa(i,ncolsp-j+1) /= 255) THEN
dBA = -6.125 + (REAL(ifield_dpa(i,ncolsp-j+1))*0.125)
! dpadata(i,j) = 10.**(dBA*0.1) ! Rainfall in mm.
dpadata(i,j) = mm2in*10.**(dBA*0.1) ! Rainfall in inches
dpadata(i,j) = MAX(0.,dpadata(i,j))
END IF
END DO ! i loop
END DO ! j loop
!------------------------------------------------------------------------
!
! Determine Cartesian coordinates of DPA points on DPA grid
! relative to the radar.
!
!------------------------------------------------------------------------
DO j = 1,nrowsp
DO i = 1,ncolsp
! dpax(i,j) = DPAdx*(REAL(i) - 66.)
! dpay(i,j) = DPAdx*(REAL(j) - 66.)
dpax(i,j) = DPAdx*(REAL(i) - 67.) ! Better match with NWS
dpay(i,j) = DPAdx*(REAL(j) - 67.) ! graphics.
END DO ! i loop
END DO ! j loop
xrad = 0. ! Radar at center of grid
yrad = 0. ! Radar at center of grid
!------------------------------------------------------------------------
!
! Determine Cartesian coordinates of ARPS scalar points in DPA map
! projection relative to radar.
!
!------------------------------------------------------------------------
CALL xytoll(nx,ny,xs,ys,arpslat,arpslon)
dpatrulat(1) = DPAtrulat1
dpatrulat(2) = DPAtrulat1 ! Not used for polar stereographic
CALL setmapr(DPAiproj,DPAscale,dpatrulat,DPAtrulon)
CALL setorig(2,radar_lat,radar_lon)
CALL lltoxy(nx,ny,arpslat,arpslon,arpsdpax,arpsdpay)
!------------------------------------------------------------------------
!
! Remap the DPA data and write to file.
!
!------------------------------------------------------------------------
CALL remapraster(nx,ny,arpsdpax,arpsdpay,dx,dpachek, &
dpamiss,xrad,yrad, &
remapopt,radius,nrowsp,ncolsp,dpadata, &
dpax,dpay,remapped_dpa,istatus)
IF (arbvaropt == 1) THEN
CALL wrtvar1(nx,ny,1,remapped_dpa,radar_name//'dp','NIDS DPA', &
'mm',000000,runname(1:lfnkey),TRIM(dir_name), &
istatus)
END IF
END IF
END DO ! DO ifilenum = 1,maxinfile
END IF ! IF (maxinfile > 0)
!------------------------------------------------------------------------
!
! Create filename for output of remapped reflectivity and velocity.
!
!------------------------------------------------------------------------
IF (count > 0) THEN
!------------------------------------------------------------------------
!
! Call remapping routines
!
!------------------------------------------------------------------------
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)
nyqset=0
IF( velproc ) THEN
timeset=0
ELSE
timeset=1
END IF
vardump=0
print *, ' Calling remapvol for reflectivity '
CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz, &
vardump,nyqset,timeset, &
refid,refname,refunits, &
refchek,refmiss,refmedl,refdazl,iordref, &
kntrgat,kntrazm,kntrelv, &
radarx,radary,radar_alt,dazim, &
rngmin,rngmax,time1st, &
rngrvol,azmrvol,elvrvol, &
refvol,timevol,nyqvvol,rxvol,ryvol,rzvol, &
xs,ys,zps,gridref,gridtim,gridnyq,istatus)
IF( velproc ) THEN
nyqset=1
timeset=1
vardump=0
WRITE(stdout,'(a)') ' Calling remapvol for velocity '
CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz, &
vardump,nyqset,timeset, &
velid,velname,velunits, &
velchek,velmiss,velmedl,veldazl,iordvel, &
kntvgat,kntvazm,kntvelv, &
radarx,radary,radar_alt,dazim, &
rngmin,rngmax,time1st, &
rngvvol,azmvvol,elvvvol, &
velvol,timevol,nyqvvol,rxvol,ryvol,rzvol, &
xs,ys,zps,gridvel,gridtim,gridnyq,istatus)
END IF
!------------------------------------------------------------------------
!
! Now use NIDS echo top and VIL fields to adjust remapped reflectivity.
!
!------------------------------------------------------------------------
! Developed for WDT.
IF (adjreflopt == 1) THEN
WRITE(stdout,'(a)') ' Adjusting reflectivity field...'
CALL adjust_refl(nx,ny,nz,gridref,zps,nidset,nidsvil)
END IF
!------------------------------------------------------------------------
!
! Create filename and write file.
!
!------------------------------------------------------------------------
CALL mkradfnm(dmpfmt,dir_name,len_dir,radar_name,iiyr,iimonth,iiday, &
iihr, iiminit, 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, &
radar_name, radar_lat, radar_lon, radar_alt, &
iyr, imonth, iday, iihr, iiminit, iisec, &
ivcp, isource, xs, ys, zps, &
gridvel, gridref, gridnyq, gridtim, &
tem1d1,tem1d2,tem1d3,tem1d4,tem1d5,tem1d6)
END IF
!------------------------------------------------------------------------
!
! The End.
!
!------------------------------------------------------------------------
WRITE(stdout,'(/a)') ' Normal termination of nids2arps.'
STOP
!------------------------------------------------------------------------
!
! Error reading namelist
!
!------------------------------------------------------------------------
999 CONTINUE
WRITE(stdout,'(/a)') ' ERROR: Cannot read NAMELIST block. Aborting...'
STOP
!------------------------------------------------------------------------
!
! Error opening file.
!
!------------------------------------------------------------------------
998 CONTINUE
WRITE(stdout,'(/a,a,/a)')' ERROR: Cannot open file ',TRIM(infile), &
' Aborting...'
STOP
END PROGRAM nids2arps