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