!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RDRADCOL                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE rdradcol(nx,ny,nz,nsrc_rad,nvar_radin,                       & 2,37
           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.
!
!-----------------------------------------------------------------------
!
!  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=132) :: 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=2
  CHARACTER (LEN=8) :: srcradin(nsrcradin)
  DATA srcradin /'88D-AII','88D-NIDS'/
!
  CHARACTER (LEN=4) :: stn
  CHARACTER (LEN=80) :: runname
  CHARACTER (LEN=132) :: fname
  INTEGER :: ireftim,itime,vcpnum,isource,idummy
  INTEGER :: iradfmt,strhopt,mapprin,dmpfmt
  INTEGER :: nchanl,ierr,nradvr,iradvr, 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,nradvr

  REAL :: dxin,dyin,dzin,dzminin,ctrlatin
  REAL :: ctrlonin,tlat1in,tlat2in,tlonin,scalin,rdummy
  REAL :: latcin,loncin
  REAL :: xrd,yrd,elev,dummy

  LOGICAL :: hdf_alloc
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  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=2
    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
!
!  Allocate hdf temporary arrays if not already done
!
      IF( .NOT. hdf_alloc ) THEN
        ALLOCATE(latradt(nx,ny),STAT=istatus)
        latradt=-999999.
        ALLOCATE(lonradt(nx,ny),STAT=istatus)
        lonradt=-999999.
        ALLOCATE (itmp(nx,ny,nz),stat=istatus)
        IF (istatus /= 0) THEN
          WRITE (6,*) "HDFREAD: ERROR allocating itmp, returning"
          STOP
        END IF
        ALLOCATE (hmax(nz),stat=istatus)
        IF (istatus /= 0) THEN
          WRITE (6,*) "HDFREAD: ERROR allocating hmax, returning"
          STOP
        END IF
        ALLOCATE (hmin(nz),stat=istatus)
        IF (istatus /= 0) THEN
          WRITE (6,*) "HDFREAD: ERROR allocating hmin, returning"
          STOP
        END IF
        hdf_alloc=.TRUE.
      END IF
!
      CALL hdfopen(trim(fname), 1, sd_id)
      IF (sd_id < 0) THEN
        WRITE (6,*) "RDRADCOL: ERROR opening ",                       &
                  trim(fname)," for reading."
        istatus = 1
        STOP
      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'
      STOP
    END IF
    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)
        GO TO 25
      END IF
    END DO
    WRITE(6,'(/a,i4,a,a)') ' rdradcol read isource: ',isource,          &
                   ' could not find srcrad= ',srcradin(isource)
    WRITE(6,'(a//)') ' Check radar source names'
    STOP
    25   CONTINUE
!
!-----------------------------------------------------------------------
!
!    Read more header data
!
!-----------------------------------------------------------------------
!
    IF (dmpfmt == 1) THEN
!
      READ(nchanl) runname
      PRINT *, ' runname: ',runname
      READ(nchanl) iradfmt,strhopt,mapprin,idummy,idummy,               &
               idummy,idummy,idummy,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, 4, '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 hdfrd1d(sd_id,'iradvr', mxradvr,iradvr,'')

    END IF

    PRINT *, ' iradfmt: ',iradfmt
    PRINT *, ' strhopt: ',strhopt
    PRINT *, ' mapprin: ',mapprin
    PRINT *, ' idummy: ',idummy
!
    PRINT *, ' dxin,dyin,dzin: ',dxin,dyin,dzin
    PRINT *, ' ctrlatin,ctrlonin: ',ctrlatin,ctrlonin
    PRINT *, ' tlat1in,tlat2in,tlonin: ',tlat1in,tlat2in,tlonin
    PRINT *, ' scalin ',scalin
    PRINT *, ' latrad,lonrad,elvrad: ',                                 &
               latrad(krad),lonrad(krad),elvrad(krad)
    PRINT *, ' rdummy: ',rdummy
    PRINT *, ' 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
!
!-----------------------------------------------------------------------
!
!  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

      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
  ELSE  ! iuserad check failed
    WRITE(6,'(a,a,/a/)')                                                &
          '   Data for source ',srcradin(isource),                      &
          ' are not used in successive correction step, skipping'
  END IF
  END DO  ! file loop
  
  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

  IF( hdf_alloc ) THEN
    DEALLOCATE (itmp)
    DEALLOCATE(latradt)
    DEALLOCATE(lonradt)
  END IF 

  RETURN
115 CONTINUE

    WRITE(6,'(/a/)') ' Error reading data in HDFREAD'
    STOP

400 CONTINUE

    WRITE(6,'(a,a,/a)')                                                 &
          '   Error opening radar file ',trim(fname),                   &
          ' Stopping in rdradcol.'
    STOP

END SUBROUTINE rdradcol