!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDRADCOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdradcol(nx,ny,nz,nsrc_rad,nvar_radin, & 4,77
mx_rad,nz_rdr,mx_colrad,mx_pass, &
radistride,radkstride, &
iuserad,iusechk,npass,nradfil,radfname, &
srcrad,isrcrad,stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,irad,nlevrad,hgtradc,obsrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads radar data stored as columns on a remapped grid.
! This allows the remapping to occur on a different grid
! than the analysis and allows the radar data to be treated
! as "soundings".
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! August, 1995
!
! MODIFICATION HISTORY:
!
! 04/25/02 (Keith Brewster)
! Added stride processing for data thinning, when desired.
!
! 05/04/02 (Leilei Wang and Keith Brewster)
! Added reading of hdf formatted files.
!
! 10/29/2002 (Keith Brewster)
! Improvement to alloc and error handling of hdf i/o.
!
!-----------------------------------------------------------------------
!
! 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
! nsrc_rad number of sources of radar data
! nvar_radin number of variables in the obsrad array
! mx_rad maximum number of radars
! nz_rdr maximum number of levels in a radar column
! mx_colrad maximum number of radar columns
!
! nradfil number of radar files
! radfname file name for radar datasets
! srcrad name of radar sources
!
! OUTPUT:
!
! isrcrad index of radar source
! stnrad radar site name character*4
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of feed horn of radar (m MSL)
! latradc latitude of radar column (degrees N)
! lonradc longitude of radar column (degrees E)
! irad radar number of each column
! nlevrad number of levels of radar data in each column
! hgtradc height (m MSL) of radar observations
! obsrad radar observations
! ncolrad number of radar columns read-in
! istatus status indicator
!
! tem1 Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: nsrc_rad,nvar_radin,mx_rad,nz_rdr,mx_colrad
INTEGER :: mx_pass
!
INTEGER :: radistride
INTEGER :: radkstride
INTEGER :: iuserad(0:nsrc_rad,mx_pass)
INTEGER :: iusechk
INTEGER :: npass
!
INTEGER :: nradfil
CHARACTER (LEN=256) :: radfname(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar site variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: srcrad(nsrc_rad)
INTEGER :: isrcrad(mx_rad)
REAL :: latrad(mx_rad),lonrad(mx_rad)
REAL :: elvrad(mx_rad)
CHARACTER (LEN=5) :: stnrad(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: irad(mx_colrad)
INTEGER :: nlevrad(mx_colrad)
REAL :: latradc(mx_colrad),lonradc(mx_colrad)
REAL :: hgtradc(nz_rdr,mx_colrad)
REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
INTEGER :: ncolrad
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny,nz)
REAL :: tem2(nx,ny,nz)
REAL :: tem3(nx,ny,nz)
REAL :: tem4(nx,ny,nz)
REAL :: tem5(nx,ny,nz)
REAL :: tem6(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! hdf variables and temporary arrays
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: mxradvr=10
INTEGER :: iradvr(mxradvr)
!
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:)
REAL, ALLOCATABLE :: hmax(:)
REAL, ALLOCATABLE :: hmin(:)
REAL, ALLOCATABLE :: latradt(:,:)
REAL, ALLOCATABLE :: lonradt(:,:)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: velchek=-200.
REAL, PARAMETER :: refchek=-20.
!
INTEGER, PARAMETER :: nsrcradin=5
CHARACTER (LEN=8) :: srcradin(nsrcradin)
DATA srcradin /'88D-AII','88D-NIDS','CASA-IP1', &
'TDWR','88D-POL'/
!
CHARACTER (LEN=4) :: stn
CHARACTER (LEN=80) :: runname
CHARACTER (LEN=256) :: fname
INTEGER :: ireftim,itime,vcpnum,isource,idummy
INTEGER :: iradfmt,strhopt,mapprin,dmpfmt
! INTEGER :: nchanl,ierr,nradvr,iradvr, ipass, ncolbgn
INTEGER :: nchanl,ierr,nradvr,ipass, ncolbgn
INTEGER :: iyr, imon, idy, ihr, imin, isec
INTEGER :: i,j,k,ivar,isrc,icstrt,icol,jcol
INTEGER :: klev,kk,krad,nfile,maxk,knt,lens,sd_id
REAL :: dxin,dyin,dzin,dzminin,ctrlatin
REAL :: ctrlonin,tlat1in,tlat2in,tlonin,scalin,rdummy
REAL :: latcin,loncin
REAL :: xrd,yrd,elev,dummy
LOGICAL :: hdf_alloc, fndsrcrad
INTEGER :: num_colrad, num_verlev
INTEGER :: typelev
REAL :: xmin, xmax, ymin, ymax
INTEGER :: numradcol, nummaxlev
INTEGER, ALLOCATABLE :: coli(:), colj(:), colk(:,:)
INTEGER, ALLOCATABLE :: numlev(:)
REAL, ALLOCATABLE :: collat(:), collon(:)
REAL, ALLOCATABLE :: radcolhgt(:,:), radcolref(:,:), radcolvel(:,:), &
radcolnyq(:,:), radcoltim(:,:)
INTEGER :: kcol
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
maxk=0
icol=0
istatus=0
icstrt=0
nfile=nradfil
hdf_alloc=.FALSE.
IF(nradfil == 0 ) THEN
ncolrad = 0
ELSE IF(nradfil > mx_rad) THEN
WRITE(6,'(a,i3,a,i3/a,i3,a)') &
' WARNING nradfil ',nradfil,' exceeds mx_rad dimension ', &
mx_rad,' only ',mx_rad,' files will be read.'
nfile=mx_rad
END IF
!
!-----------------------------------------------------------------------
!
! Loop through all radars
!
!-----------------------------------------------------------------------
!
DO krad=1,nfile
fname=radfname(krad)
CALL asnctl
('NEWLOCAL', 1, ierr)
CALL asnfile(fname, '-F f77 -N ieee', ierr)
lens=LEN(trim(fname))
IF (fname(lens-3:lens) == 'hdf4') THEN
dmpfmt=3
ELSE
dmpfmt=1
ENDIF
write(6,'(a,i4,a,a)') &
' rdradcol: dmpfmt=',dmpfmt,' fname=',TRIM(fname)
IF (dmpfmt == 1) THEN
WRITE(6,'(/a,a)') &
' Opening radar file ',trim(fname)
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(fname),ERR=400, &
FORM='unformatted',STATUS='old')
istatus=1
READ(nchanl) stn
stnrad(krad)=stn
READ(nchanl) ireftim,itime,vcpnum,isource,idummy, &
idummy,idummy,idummy,idummy,idummy
ELSE
IF (.NOT. ALLOCATED(itmp)) THEN
ALLOCATE (itmp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:itmp")
END IF
CALL hdfopen(trim(fname), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,'(3a)') 'RDRADCOL: ERROR opening hdf file:', &
trim(fname),' for reading.'
istatus = 0
goto 700
END IF
CALL hdfrdc(sd_id, 4, 'radid', stn, istatus)
stnrad(krad)=stn
CALL hdfrdi(sd_id, 'ireftim', ireftim, istatus)
CALL hdfrdi(sd_id, 'itime', itime, istatus)
CALL hdfrdi(sd_id, 'vcpnum', vcpnum, istatus)
CALL hdfrdi(sd_id, 'isource', isource, istatus)
END IF
!
!-----------------------------------------------------------------------
!
! Connect source number to source name.
!
!-----------------------------------------------------------------------
!
IF(isource < 1) THEN
WRITE(6,'(a,i5,/a/)') &
' rdradcol: read isource as ',isource, &
' Setting isource to default, 1'
isource=1
END IF
IF(isource > nsrcradin) THEN
WRITE(6,'(a,i3,a,i3,a,/a/)') &
' rdradcol: read isource as ',isource, &
' only ',nsrcradin,' sources known.', &
' If a new radar source has been added, update rdradcol.f90'
CALL arpsstop('nsrcradin is too small.',1)
END IF
fndsrcrad = .FALSE.
DO isrc=1,nsrc_rad
IF(srcrad(isrc) == srcradin(isource)) THEN
isrcrad(krad)=isrc
WRITE(6,'(a,i3,a,i3,a,a)') &
' rdradcol isource: ',isource,' = srcrad(',isrc,') ',srcrad(isrc)
fndsrcrad = .TRUE.
EXIT
END IF
END DO
IF (.NOT. fndsrcrad) THEN
WRITE(6,'(/a,i4,a,a)') ' rdradcol read isource: ',isource, &
' could not find srcrad= ',srcradin(isource)
WRITE(6,'(a//)') ' Check radar source names'
CALL arpsstop('Wrong radar source names.',1)
END IF
!
!-----------------------------------------------------------------------
!
! Read more header data
!
!-----------------------------------------------------------------------
!
IF (dmpfmt == 1) THEN
READ(nchanl) runname
READ(nchanl) iradfmt,strhopt,mapprin,idummy,idummy, &
typelev,numradcol,nummaxlev,idummy,idummy
READ(nchanl) dxin,dyin,dzin,dzminin,ctrlatin, &
ctrlonin,tlat1in,tlat2in,tlonin,scalin, &
latrad(krad),lonrad(krad),elvrad(krad), &
rdummy,rdummy
READ(nchanl) nradvr,iradvr
ELSE
CALL hdfrdc(sd_id, 40, 'runname', runname, istatus)
CALL hdfrdi(sd_id, 'iradfmt', iradfmt, istatus)
CALL hdfrdi(sd_id, 'strhopt', strhopt, istatus)
CALL hdfrdi(sd_id, 'mapproj', mapprin, istatus)
CALL hdfrdr(sd_id, 'dx', dxin, istatus)
CALL hdfrdr(sd_id, 'dy', dyin, istatus)
CALL hdfrdr(sd_id, 'dz', dzin, istatus)
CALL hdfrdr(sd_id, 'dzmin', dzminin, istatus)
CALL hdfrdr(sd_id, 'ctrlat', ctrlatin, istatus)
CALL hdfrdr(sd_id, 'ctrlon', ctrlonin, istatus)
CALL hdfrdr(sd_id, 'trulat1', tlat1in, istatus)
CALL hdfrdr(sd_id, 'trulat2', tlat2in, istatus)
CALL hdfrdr(sd_id, 'trulon', tlonin, istatus)
CALL hdfrdr(sd_id, 'sclfct', scalin, istatus)
CALL hdfrdr(sd_id, 'latrad', latrad(krad), istatus)
CALL hdfrdr(sd_id, 'lonrad', lonrad(krad), istatus)
CALL hdfrdr(sd_id, 'elvrad', elvrad(krad), istatus)
CALL hdfrdi(sd_id, 'nradvr', nradvr, istatus)
CALL hdfrd1di(sd_id,'iradvr', mxradvr,iradvr,istatus)
CALL hdfrdi(sd_id, 'typelev', typelev, istatus)
END IF
WRITE(6,'(1x,2a)') 'runname: ',runname
WRITE(6,'(1x,3(a,I2))') 'iradfmt: ',iradfmt,', strhopt: ',strhopt,', mapprin: ',mapprin
WRITE(6,'(1x,a,3F8.0)') 'dxin,dyin,dzin: ',dxin,dyin,dzin
WRITE(6,'(1x,a,2F8.2)') 'ctrlatin,ctrlonin: ',ctrlatin,ctrlonin
WRITE(6,'(1x,a,3F8.2)') 'tlat1in,tlat2in,tlonin: ',tlat1in,tlat2in,tlonin
WRITE(6,'(1x,a,F8.0)') 'scalin: ',scalin
WRITE(6,'(1x,a,3F8.2)') 'latrad,lonrad,elvrad: ', latrad(krad),lonrad(krad),elvrad(krad)
WRITE(6,'(1x,a,20I4)') 'Got nradvr,iradvr: ',nradvr,iradvr
CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
iyr=MOD(iyr,100)
WRITE(6,'(/a,i2.2,a,i2.2,a,i2.2,1X,i2.2,a,i2.2,a)') &
' Reading remapped raw radar data for: ', &
imon,'-',idy,'-',iyr,ihr,':',imin,' UTC'
!
!-----------------------------------------------------------------------
!
! If this radar source is not going to be used in the successive
! corrections analysis step, then don't bother reading columns.
!
!-----------------------------------------------------------------------
!
IF(iusechk > 0 ) THEN
knt=0
DO ipass=1,npass
knt=knt+iuserad(isource,ipass)
END DO
ELSE
knt=1
END IF
IF( knt > 0 ) THEN
IF (typelev == 1) THEN ! new format with ARPS vertical heights
IF (dmpfmt == 1) THEN
READ(nchanl) xmin, xmax, ymin, ymax
ELSE
CALL hdfrdi(sd_id, 'numradcol', numradcol, istatus)
CALL hdfrdi(sd_id, 'nummaxelv', nummaxlev, istatus)
CALL hdfrdr(sd_id, 'xmin', xmin, istatus)
CALL hdfrdr(sd_id, 'xmax', xmax, istatus)
CALL hdfrdr(sd_id, 'ymin', xmin, istatus)
CALL hdfrdr(sd_id, 'ymax', ymax, istatus)
END IF
! CALL radinf(rlim)
! if this radar is in the influence of radius
IF (icol+numradcol > mx_colrad ) THEN
WRITE(6,'(1x,a,I4,a,/,a,/)') &
'Radar column number mx_colrad (',mx_colrad,') is too small.', &
'Please change the value in adas.inc and recompile the program.'
CALL arpsstop('ERROR: Dimension mx_colrad in adas.inc is too small.',1)
END IF
IF (nummaxlev > nz_rdr) THEN
WRITE(6,'(1x,a,I4,a,/,a,/)') &
'Radar column maximum vertical level nz_rdr (',nz_rdr,') is too small.', &
'Please change the value in adas.inc and recompile the program.'
CALL arpsstop('ERROR: Dimension nz_rdr in adas.inc is too small.',1)
END IF
ALLOCATE(coli(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:coli")
ALLOCATE(colj(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:colj")
ALLOCATE(colk(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:colk")
ALLOCATE(numlev(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:numlev")
ALLOCATE(collat(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:collat")
ALLOCATE(collon(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:collon")
ALLOCATE(radcolhgt(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:radcolhgt")
ALLOCATE(radcolref(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:radcolref")
ALLOCATE(radcolvel(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:radcolvel")
ALLOCATE(radcolnyq(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:radcolnyq")
ALLOCATE(radcoltim(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "rdradcol:radcoltim")
radcolhgt(:,:) = -999.
radcolref(:,:) = -999.
radcolvel(:,:) = -999.
radcolnyq(:,:) = -999.
radcoltim(:,:) = -999.
IF (dmpfmt == 1) THEN ! READ binary files
kcol = 0
DO jcol=1,numradcol
READ(nchanl,END=201) i,j,xrd,yrd, &
latcin,loncin,elev,klev
kcol = kcol + 1
coli(kcol) = i
colj(kcol) = j
collat(kcol) = latcin
collon(kcol) = loncin
numlev(kcol) = klev
READ(nchanl,END=202) (colk(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolhgt(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolref(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolvel(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolnyq(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcoltim(kk,kcol),kk=1,klev)
END DO
201 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached after reading', &
(kcol-icstrt),' columns'
GO TO 205
202 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached while reading', &
kcol,' column'
205 CONTINUE
CLOSE(nchanl)
CALL retunit( nchanl )
nummaxlev = MAX(nummaxlev,klev)
ELSE ! READ hdf file
CALL hdfrd1d (sd_id,'radcollat',numradcol,collat,istatus)
CALL hdfrd1d (sd_id,'radcollon',numradcol,collon,istatus)
CALL hdfrd1di(sd_id,'numelev', numradcol,numlev,istatus)
CALL hdfrd1di(sd_id,'radcoli', numradcol,coli,istatus)
CALL hdfrd1di(sd_id,'radcolj', numradcol,colj,istatus)
CALL hdfrd2di(sd_id,'radcolk', nummaxlev,numradcol,colk,istatus)
CALL hdfrd2d (sd_id,'radcolhgt',nummaxlev,numradcol,radcolhgt,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolref',nummaxlev,numradcol,radcolref,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolvel',nummaxlev,numradcol,radcolvel,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolnyq',nummaxlev,numradcol,radcolnyq,istatus,itmp)
CALL hdfrd2d (sd_id,'radcoltim',nummaxlev,numradcol,radcoltim,istatus,itmp)
CALL hdfclose(sd_id,istatus)
END IF
DO kcol = 1, numradcol
IF (MOD(coli(kcol)+colj(kcol),radistride) == 0) THEN
icol = icol + 1
irad(icol) = krad
latradc(icol) = collat(kcol)
lonradc(icol) = collon(kcol)
nlevrad(icol) = numlev(kcol)
DO kk = 1,numlev(kcol)
IF ( MOD( (colk(kk,kcol)-1), radkstride ) == 0 .OR. &
kk == numlev(kcol) ) THEN
hgtradc (kk,icol) = radcolhgt(kk,kcol)
obsrad(1,kk,icol) = radcolref(kk,kcol)
obsrad(2,kk,icol) = radcolvel(kk,kcol)
obsrad(3,kk,icol) = radcolnyq(kk,kcol)
obsrad(4,kk,icol) = radcoltim(kk,kcol)
END IF
END DO
END IF ! horizontal stride
END DO
WRITE(6,'(a,i6,a)') ' Read ',(icol-icstrt), &
' non-missing columns from hdf file.'
icstrt = icol
maxk = nummaxlev
DEALLOCATE(coli, colj, colk,numlev)
DEALLOCATE(collat, collon)
DEALLOCATE(radcolhgt, radcolref, radcolvel, radcolnyq, radcoltim)
ELSE ! Maybe old format
!
!-----------------------------------------------------------------------
!
! Note here the radar data indices:
!
! 1 Reflectivity
! 2 Radial Velocity
! 3 Nyquist Velocity
! 4 Time in seconds from the reference time
!
!-----------------------------------------------------------------------
!
IF (dmpfmt == 1) THEN
DO jcol=1,999999
READ(nchanl,END=101) i,j,xrd,yrd, &
latcin,loncin,elev,klev
READ(nchanl,END=102) (tem1(1,1,kk),kk=1,klev)
READ(nchanl,END=102) (tem2(1,1,kk),kk=1,klev)
READ(nchanl,END=102) (tem3(1,1,kk),kk=1,klev)
READ(nchanl,END=102) (tem4(1,1,kk),kk=1,klev)
READ(nchanl,END=102) (tem5(1,1,kk),kk=1,klev)
READ(nchanl,END=102) (tem6(1,1,kk),kk=1,klev)
IF(mod((i+j),radistride) == 0 ) THEN
icol=icol+1
IF(icol <= mx_colrad) THEN
irad(icol)=krad
latradc(icol)=latcin
lonradc(icol)=loncin
DO kk=1,nz_rdr
hgtradc(kk,icol)=-999.
DO ivar=1,nvar_radin
obsrad(ivar,kk,icol)=-999.
END DO
END DO
IF(klev > 1) THEN
DO kk=1,(klev-1),radkstride
k=nint(tem1(1,1,kk))
maxk=MAX(maxk,k)
k=MIN(k,nz_rdr)
hgtradc(k,icol)=tem2(1,1,kk)
obsrad(1,k,icol)=tem3(1,1,kk)
obsrad(2,k,icol)=tem4(1,1,kk)
obsrad(3,k,icol)=tem5(1,1,kk)
obsrad(4,k,icol)=tem6(1,1,kk)
nlevrad(icol)=k
END DO
END IF
k=nint(tem1(1,1,klev))
maxk=MAX(maxk,k)
k=MIN(k,nz_rdr)
hgtradc(k,icol)=tem2(1,1,klev)
obsrad(1,k,icol)=tem3(1,1,klev)
obsrad(2,k,icol)=tem4(1,1,klev)
obsrad(3,k,icol)=tem5(1,1,klev)
obsrad(4,k,icol)=tem6(1,1,klev)
nlevrad(icol)=k
ELSE
maxk=MAX(maxk,nint(tem1(1,1,klev)))
END IF
END IF
END DO
101 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached after reading', &
(icol-icstrt),' columns'
GO TO 105
102 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached while reading', &
icol,' column'
icol=icol-1
105 CONTINUE
CLOSE(nchanl)
CALL retunit( nchanl )
icstrt=icol+1
ELSE ! hdf file
!
! Allocate hdf temporary arrays if not already done
!
IF( .NOT. hdf_alloc ) THEN
ALLOCATE(latradt(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:latradt")
latradt=-999999.
ALLOCATE(lonradt(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:lonradt")
lonradt=-999999.
ALLOCATE (hmax(nz),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:hmax")
ALLOCATE (hmin(nz),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:hmin")
hdf_alloc=.TRUE.
END IF
CALL hdfrd2d(sd_id,'grdlatc',nx,ny,latradt,istatus,itmp)
IF (istatus /= 0) GO TO 115
CALL hdfrd2d(sd_id,'grdlonc',nx,ny,lonradt,istatus,itmp)
IF (istatus /= 0) GO TO 115
CALL hdfrd3d(sd_id,'hgtrad',nx,ny,nz,tem2,istatus,itmp,hmax,hmin)
IF (istatus /= 0) GO TO 115
CALL hdfrd3d(sd_id,'gridref',nx,ny,nz,tem3,istatus,itmp,hmax,hmin)
IF (istatus /= 0) GO TO 115
CALL hdfrd3d(sd_id,'gridvel',nx,ny,nz,tem4,istatus,itmp,hmax,hmin)
IF (istatus /= 0) GO TO 115
CALL hdfrd3d(sd_id,'gridnyq',nx,ny,nz,tem5,istatus,itmp,hmax,hmin)
IF (istatus /= 0) GO TO 115
CALL hdfrd3d(sd_id,'gridtim',nx,ny,nz,tem6,istatus,itmp,hmax,hmin)
IF (istatus /= 0) GO TO 115
DO j=1,ny
DO i=1,nx
IF(mod((i+j),radistride) == 0 ) THEN
!
! Check for non-missing data in this column
!
knt=0
DO k=1,(nz-1),radkstride
IF(tem3(i,j,k) > refchek .OR. &
tem4(i,j,k) > velchek ) THEN
maxk=MAX(maxk,k)
knt=knt+1
END IF
END DO
!
! If non-missing data exist increment column counter and
! transfer to column data
!
IF( knt > 0 ) THEN
icol=icol+1
IF(icol <= mx_colrad) THEN
irad(icol)=krad
latradc(icol)=latradt(i,j)
lonradc(icol)=lonradt(i,j)
DO kk=1,nz_rdr
hgtradc(kk,icol)=-999.
DO ivar=1,nvar_radin
obsrad(ivar,kk,icol)=-999.
END DO
END DO
DO k=1,(nz-1),radkstride
IF(tem3(i,j,k) > refchek .OR. &
tem4(i,j,k) > velchek ) THEN
hgtradc(k,icol)=tem2(i,j,k)
obsrad(1,k,icol)=tem3(i,j,k)
obsrad(2,k,icol)=tem4(i,j,k)
obsrad(3,k,icol)=tem5(i,j,k)
obsrad(4,k,icol)=tem6(i,j,k)
nlevrad(icol)=k
END IF
END DO
END IF ! icol less than mx_colrad
END IF ! non-missing data in column
END IF ! stride check
END DO
END DO
CALL hdfclose(sd_id,istatus)
WRITE(6,'(a,i6,a)') ' Read ',(icol-icstrt), &
' non-missing columns from hdf file.'
icstrt=icol+1
END IF ! HDF file
END IF ! New format
ELSE ! iuserad check failed
IF( dmpfmt == 1 ) THEN
WRITE(6,'(a,a,/a/)') &
' Data for source ',srcradin(isource), &
' are not used in successive correction step, skipping'
CLOSE(nchanl)
CALL retunit( nchanl )
ELSE ! hdf format file
CALL hdfclose(sd_id,istatus)
END IF
END IF ! use source ?
700 CONTINUE
END DO ! file loop
IF (ALLOCATED(itmp)) DEALLOCATE(itmp)
IF (ALLOCATED(latradt)) THEN
DEALLOCATE(latradt, lonradt)
DEALLOCATE(hmax, hmin)
END IF
IF(icol > mx_colrad) THEN
WRITE(6,'(//a)') ' #### WARNING max_colrad EXCEEDED ####'
WRITE(6,'(a,i6,a)') ' increase mx_colrad from ',mx_colrad,' columns'
WRITE(6,'(2x,i6,a//)') icol,' columns are needed to read all files'
END IF
ncolrad=MIN(icol,mx_colrad)
WRITE(6,'(a,i5/)') ' Maximum number of vert levels read: ',maxk
IF(maxk > nz_rdr) THEN
WRITE(6,'(a,i5)') ' EXCEEDS nz_rdr, increase nz_rdr: ',nz_rdr
STOP
END IF
RETURN
115 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in RDRADCOL(HDF format).'
STOP
400 CONTINUE
WRITE(6,'(a,a,/a)') ' Error opening radar file ',trim(fname), &
' Stopping in rdradcol.'
STOP
END SUBROUTINE rdradcol
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FILL1RAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE fill1rad(nx,ny,nz,radfname,lvldbg,xs,ys,zs, & 1,67
strhopt_rad,mapproj_rad,dx_rad,dy_rad,dz_rad, &
dzmin_rad,ctrlat_rad,ctrlon_rad, &
tlat1_rad,tlat2_rad,tlon_rad,sclfct_rad, &
isrcrad,stnrad,latrad,lonrad,elvrad, &
gridvel,gridref,gridnyq,gridtim, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads radar data remapped on the ARPS grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yunheng Wang 07/02/2007
!
! 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
! radfname file name for radar datasets
! lvldbg Debug level
!
! OUTPUT:
!
! isrcrad index of radar source
! stnrad radar site name character*4
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of feed horn of radar (m MSL)
!
! gridvel radial velocity on ARPS grid
! gridref reflectivity on ARPS grid
! gridnyq nyquist velocity on ARPS grid
! gridtim observation time at ARPS grid
! temrad temporary radar array
!
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INCLUDE 'mp.inc'
!
INTEGER :: nx,ny,nz ! the ARPS local grid size
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zs(nx,ny,nz)
CHARACTER (LEN=132) :: radfname
INTEGER :: lvldbg
INTEGER :: strhopt_rad,mapproj_rad
REAL :: dx_rad,dy_rad,dz_rad,dzmin_rad
REAL :: ctrlat_rad,ctrlon_rad,tlat1_rad,tlat2_rad,tlon_rad
REAL :: sclfct_rad
!
! OUTPUT: Radar site variables
!
INTEGER :: isrcrad
CHARACTER (LEN=5) :: stnrad
REAL :: latrad
REAL :: lonrad
REAL :: elvrad
!
! OUTPUT: ARPS radar arrays
!
REAL :: gridvel(nx,ny,nz)
REAL :: gridref(nx,ny,nz)
REAL :: gridnyq(nx,ny,nz)
REAL :: gridtim(nx,ny,nz)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: readk(nz)
REAL :: readhgt(nz)
REAL :: readref(nz)
REAL :: readvel(nz)
REAL :: readnyq(nz)
REAL :: readtim(nz)
!
INTEGER, ALLOCATABLE :: kntref(:)
INTEGER, ALLOCATABLE :: kntvel(:)
!
INTEGER :: mxradvr,nradvr
PARAMETER(mxradvr=10)
INTEGER :: iradvr(mxradvr)
CHARACTER (LEN=4) :: stn
CHARACTER (LEN=80) :: runname
INTEGER :: ireftim,itime,vcpnum,idummy
INTEGER :: iradfmt,strhoptin,mapprin
INTEGER :: nchanl,ierr,iostatus
INTEGER :: iyr, imon, idy, ihr, imin, isec
INTEGER :: i,j,k,irad,jrad,krad,kk,ipt,klev
INTEGER :: irngmin,irngmax
INTEGER :: ilen,jlen,istart,iend,jstart,jend
REAL :: height,rmin,rmax,rmax2,xrad,yrad,dist,dist2,elev,range
REAL :: refelvmin,refelvmax,refrngmin,refrngmax
REAL :: xrd,yrd,gridlat,gridlon,rdummy
LOGICAL :: verbose = .FALSE.
INTEGER :: ips, ipe, jps, jpe
!
REAL, PARAMETER :: defelvmin = 0.5
REAL, PARAMETER :: defelvmax = 19.5
REAL, PARAMETER :: defrmin = 10.E03
REAL, PARAMETER :: defrmax = 230.E03
!
!-----------------------------------------------------------------------
!
! Temporary arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: typelev
REAL :: xmin, xmax, ymin, ymax
INTEGER :: numradcol, nummaxlev
INTEGER, ALLOCATABLE :: coli(:), colj(:), colk(:,:)
INTEGER, ALLOCATABLE :: numlev(:)
REAL, ALLOCATABLE :: collat(:), collon(:)
REAL, ALLOCATABLE :: radcolhgt(:,:), radcolref(:,:), radcolvel(:,:), &
radcolnyq(:,:), radcoltim(:,:)
INTEGER :: kcol
INTEGER :: iproc, jproc
!-----------------------------------------------------------------------
!
! hdf arrays
!
!-----------------------------------------------------------------------
!
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:)
INTEGER :: lens,dmpfmt,sd_id,isource
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
istatus=0
refelvmin=-999.
refelvmax=-999.
refrngmin=-999.
refrngmax=-999.
numradcol=0
nummaxlev=0
!
!-----------------------------------------------------------------------
!
! Initializations
!
!-----------------------------------------------------------------------
!
IF (lvldbg > 90) THEN
verbose = .TRUE.
IF (mp_opt > 0) THEN
ips = 2 ! Patch start index
ipe = nx-2 ! Patch end index
jps = 2
jpe = ny-2
IF (loc_x == 1) ips = 1 ! To cover the same domain as non-MPI runs
IF (loc_x == nproc_x) ipe = nx
IF (loc_y == 1) jps = 1
IF (loc_y == nproc_y) jpe = ny
ELSE
ips = 1 ! Domain start index
ipe = nx ! Domain end index
jps = 1
jpe = ny
END IF
ALLOCATE(kntref(nz))
ALLOCATE(kntvel(nz))
DO k=1,nz
kntref(k)=0
kntvel(k)=0
END DO
END IF
!
! OpenMP changed loop order to j,k,i:
!$OMP PARALLEL DO PRIVATE(j,k,i)
DO j=1,ny
DO k=1,nz
DO i=1,nx
gridref(i,j,k)=-999.
gridvel(i,j,k)=-999.
gridnyq(i,j,k)=-999.
gridtim(i,j,k)=-999.
END DO
END DO
END DO
!
!
!-----------------------------------------------------------------------
!
! Open file
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(radfname, '-F f77 -N ieee', ierr)
lens=LEN(trim(radfname))
IF(radfname(lens-3:lens)=='hdf4')THEN
dmpfmt=3
ELSE
dmpfmt=1
ENDIF
print *, myproc,': radfname = ',TRIM(radfname),', dmpfmt = ',dmpfmt
IF( dmpfmt == 1 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(radfname),iostat=iostatus, &
FORM='unformatted',STATUS='old')
IF(iostatus /= 0) THEN
WRITE(6,'(a,a)') 'Error opening the radar file:', &
TRIM(radfname)
istatus=-1
GO TO 999
END IF
!
!-----------------------------------------------------------------------
!
! Read radar description variables
!
!-----------------------------------------------------------------------
!
istatus=1
READ(nchanl) stn
stnrad=stn
READ(nchanl) ireftim,itime,vcpnum,isrcrad,idummy, &
idummy,idummy,idummy,idummy,idummy
READ(nchanl) runname
READ(nchanl) iradfmt,strhopt_rad,mapproj_rad,irngmin,irngmax, &
! idummy,idummy,idummy,idummy,idummy
typelev,numradcol,nummaxlev,idummy,idummy
READ(nchanl) dx_rad,dy_rad,dz_rad,dzmin_rad,ctrlat_rad, &
ctrlon_rad,tlat1_rad,tlat2_rad,tlon_rad,sclfct_rad, &
latrad,lonrad,elvrad,refelvmin,refelvmax
READ(nchanl) nradvr,iradvr
READ(nchanl) xmin, xmax, ymin, ymax
ELSE
! IF (.NOT. ALLOCATED(itmp)) THEN
! ALLOCATE (itmp(nx,ny,nz),stat=istatus)
! CALL check_alloc_status(istatus, "rdradcol:itmp")
! END IF
CALL hdfopen(trim(radfname), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "FILL1RAD: ERROR opening ", &
trim(radfname)," for reading."
istatus = -1
GO TO 999
ELSE
istatus = 0
END IF
CALL hdfrdc(sd_id, 4, 'radid', stn, istatus)
stnrad=stn
CALL hdfrdi(sd_id, 'ireftim', ireftim, istatus)
CALL hdfrdi(sd_id, 'itime', itime, istatus)
CALL hdfrdi(sd_id, 'vcpnum', vcpnum, istatus)
CALL hdfrdi(sd_id, 'isource', isrcrad, istatus)
CALL hdfrdc(sd_id, 40, 'runname', runname, istatus)
CALL hdfrdi(sd_id, 'iradfmt', iradfmt, istatus)
CALL hdfrdi(sd_id, 'strhopt', strhopt_rad, istatus)
CALL hdfrdi(sd_id, 'mapproj', mapproj_rad, istatus)
CALL hdfrdr(sd_id, 'dx', dx_rad, istatus)
CALL hdfrdr(sd_id, 'dy', dy_rad, istatus)
CALL hdfrdr(sd_id, 'dz', dz_rad, istatus)
CALL hdfrdr(sd_id, 'dzmin', dzmin_rad, istatus)
CALL hdfrdr(sd_id, 'ctrlat', ctrlat_rad, istatus)
CALL hdfrdr(sd_id, 'ctrlon', ctrlon_rad, istatus)
CALL hdfrdr(sd_id, 'trulat1', tlat1_rad, istatus)
CALL hdfrdr(sd_id, 'trulat2', tlat2_rad, istatus)
CALL hdfrdr(sd_id, 'trulon', tlon_rad, istatus)
CALL hdfrdr(sd_id, 'sclfct', sclfct_rad, istatus)
CALL hdfrdr(sd_id, 'latrad', latrad, istatus)
CALL hdfrdr(sd_id, 'lonrad', lonrad, istatus)
CALL hdfrdr(sd_id, 'elvrad', elvrad, istatus)
CALL hdfrdi(sd_id, 'irngmin', irngmin, istatus)
CALL hdfrdi(sd_id, 'irngmax', irngmax, istatus)
CALL hdfrdr(sd_id, 'refelvmin', refelvmin, istatus)
CALL hdfrdr(sd_id, 'refelvmax', refelvmax, istatus)
CALL hdfrdi(sd_id, 'nradvr', nradvr, istatus)
CALL hdfrd1di(sd_id,'iradvr', mxradvr,iradvr,istatus)
IF (verbose) PRINT *, ' Got nradvr,iradvr: ',nradvr,iradvr
CALL hdfrdi(sd_id, 'typelev', typelev, istatus)
CALL hdfrdi(sd_id, 'numradcol', numradcol, istatus)
CALL hdfrdi(sd_id, 'nummaxelv', nummaxlev, istatus)
CALL hdfrdr(sd_id, 'xmin', xmin, istatus)
CALL hdfrdr(sd_id, 'xmax', xmax, istatus)
CALL hdfrdr(sd_id, 'ymin', xmin, istatus)
CALL hdfrdr(sd_id, 'ymax', ymax, istatus)
END IF
CALL abss2ctim(itime, iyr, imon, idy, ihr, imin, isec )
iyr=MOD(iyr,100)
WRITE(6,'(/a,i2.2,a,i2.2,a,i2.2,1X,i2.2,a,i2.2,a)') &
'Reading remapped raw radar data for: ', &
imon,'/',idy,'/',iyr,ihr,':',imin,' UTC'
IF (typelev /= 1) THEN
WRITE(6,'(1x,a,/,1x,2a,/)') &
'Sub fill1rad was rewritten to support new radar data format. ', &
'In order to read data in old format, please change all calls of',&
' fill1rad into fill1rad_old in the source code.'
CALL arpsstop('Maybe in old data format?',1)
! istatus = -1
! GO TO 888
END IF
!
!-----------------------------------------------------------------------
!
! Make sure we don't have a corrupt file.
!
!-----------------------------------------------------------------------
!
IF (numradcol < 1 .OR. nummaxlev < 1) THEN
WRITE(6,'(//,a)') 'WARNING: Data corruption. File not processed.'
istatus = -1
GO TO 888
END IF
ALLOCATE(coli(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:coli")
ALLOCATE(colj(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:colj")
ALLOCATE(colk(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:colk")
ALLOCATE(numlev(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:numlev")
ALLOCATE(collat(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:collat")
ALLOCATE(collon(numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:collon")
ALLOCATE(radcolhgt(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:radcolhgt")
ALLOCATE(radcolref(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:radcolref")
ALLOCATE(radcolvel(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:radcolvel")
ALLOCATE(radcolnyq(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:radcolnyq")
ALLOCATE(radcoltim(nummaxlev,numradcol), STAT = istatus)
CALL check_alloc_status(istatus, "fill1rad:radcoltim")
radcolhgt(:,:) = -999.
radcolref(:,:) = -999.
radcolvel(:,:) = -999.
radcolnyq(:,:) = -999.
radcoltim(:,:) = -999.
IF (dmpfmt == 1) THEN ! READ binary files
kcol = 0
DO kcol=1,numradcol
READ(nchanl,END=201) i,j,xrd,yrd, &
latrad,lonrad,elev,klev
coli(kcol) = i
colj(kcol) = j
collat(kcol) = latrad
collon(kcol) = lonrad
numlev(kcol) = klev
READ(nchanl,END=202) (colk(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolhgt(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolref(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolvel(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcolnyq(kk,kcol),kk=1,klev)
READ(nchanl,END=202) (radcoltim(kk,kcol),kk=1,klev)
END DO
201 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached after reading', &
kcol,' columns'
GO TO 205
202 CONTINUE
WRITE(6,'(a,i6,a)') ' End of file reached while reading', &
kcol,' column'
205 CONTINUE
! CLOSE(nchanl)
! CALL retunit( nchanl )
nummaxlev = MAX(nummaxlev,klev)
ELSE ! READ hdf file
CALL hdfrd1d (sd_id,'radcollat',numradcol,collat,istatus)
CALL hdfrd1d (sd_id,'radcollon',numradcol,collon,istatus)
CALL hdfrd1di(sd_id,'numelev', numradcol,numlev,istatus)
CALL hdfrd1di(sd_id,'radcoli', numradcol,coli,istatus)
CALL hdfrd1di(sd_id,'radcolj', numradcol,colj,istatus)
CALL hdfrd2di(sd_id,'radcolk', nummaxlev,numradcol,colk,istatus)
IF (.NOT. ALLOCATED(itmp)) THEN
ALLOCATE (itmp(nummaxlev,numradcol),stat=istatus)
CALL check_alloc_status(istatus, "rdradcol:itmp")
END IF
CALL hdfrd2d (sd_id,'radcolhgt',nummaxlev,numradcol,radcolhgt,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolref',nummaxlev,numradcol,radcolref,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolvel',nummaxlev,numradcol,radcolvel,istatus,itmp)
CALL hdfrd2d (sd_id,'radcolnyq',nummaxlev,numradcol,radcolnyq,istatus,itmp)
CALL hdfrd2d (sd_id,'radcoltim',nummaxlev,numradcol,radcoltim,istatus,itmp)
! CALL hdfclose(sd_id,istatus)
END IF
DO kcol = 1, numradcol
iproc = (coli(kcol)-2)/(nx-3) + 1
jproc = (colj(kcol)-2)/(ny-3) + 1
IF (loc_x == iproc .AND. loc_y == jproc) THEN
i = MOD((coli(kcol)-2),(nx-3)) + 2
j = MOD((colj(kcol)-2),(ny-3)) + 2
klev = numlev(kcol)
DO kk=1,klev
k=colk(kk,kcol)
IF(k <= nz.AND.k >= 1) THEN
gridref(i,j,k)=max(radcolref(kk,kcol),gridref(i,j,k))
gridvel(i,j,k)=radcolvel(kk,kcol)
gridnyq(i,j,k)=radcolnyq(kk,kcol)
gridtim(i,j,k)=radcoltim(kk,kcol)
END IF ! 1 < k < nz
END DO ! kk = 1, klev
END IF ! 1 < i < nxlg & 1 < j < nylg
END DO
IF (verbose) THEN
DO j=jps,jpe
DO k=1,nz
DO i=ips,ipe
IF (gridref(i,j,k) > -200. .AND. gridref(i,j,k) < 200.) &
kntref(k)=kntref(k)+1
IF (gridvel(i,j,k) > -200. .AND. gridvel(i,j,k) < 200.) &
kntvel(k)=kntvel(k)+1
END DO
END DO
END DO
DO k=1,nz
CALL mptotali(kntref(k))
CALL mptotali(kntvel(k))
END DO
IF (myproc == 0) THEN
WRITE(6,'(a)') ' k n ref n vel'
DO k = 1,nz
WRITE(6,'(i3,2i10)') k,kntref(k),kntvel(k)
END DO
END IF
END IF
!
! Label 999 is for opens that failed.
! Label 888 is for all successes and those failures after open works.
!
888 CONTINUE
IF ( dmpfmt == 1 ) THEN
CLOSE(nchanl)
CALL retunit( nchanl )
ELSE
CALL hdfclose(sd_id,istatus)
END IF
999 CONTINUE
IF (verbose) THEN
DEALLOCATE(kntref)
DEALLOCATE(kntvel)
END IF
IF (ALLOCATED(itmp)) DEALLOCATE(itmp)
IF (ALLOCATED(coli)) THEN
DEALLOCATE(coli, colj, colk, numlev)
DEALLOCATE(collat, collon)
DEALLOCATE(radcolhgt, radcolref, radcolvel, radcolnyq, radcoltim)
END IF
RETURN
END SUBROUTINE fill1rad