!
!##################################################################
!##################################################################
!######                                                      ######
!######                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