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