!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKRADFNM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkradfnm(dmpfmt,dir,ldir,radar,iyr,imo,ida,ihr,imin,isec, & 2,2
fname,lfnm)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=256) :: dir
INTEGER :: ldir
CHARACTER (LEN=4) :: radar
INTEGER :: iyr
INTEGER :: imo
INTEGER :: ida
INTEGER :: ihr
INTEGER :: imin
INTEGER :: isec
INTEGER :: dmpfmt
!
!-----------------------------------------------------------------------
!
! Output arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=256) :: fname
INTEGER :: lfnm
!
INTEGER :: rdtime,jyr,jmo,jda,jhr,jmin,jsec
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Round to nearest minute
!
!-----------------------------------------------------------------------
!
jyr=iyr+1900
IF(iyr < 50) jyr=jyr+100
CALL ctim2abss
( jyr,imo,ida,ihr,imin,isec, rdtime )
IF(isec >= 30) rdtime=rdtime+(60-isec)
CALL abss2ctim( rdtime, jyr, jmo, jda, jhr, jmin, jsec )
jyr=MOD(jyr,100)
IF(dmpfmt==1)THEN
WRITE(fname,'(a,a4,a1,3(i2.2),a1,2(i2.2))') &
dir(1:ldir),radar,'.', &
jyr,jmo,jda,'.',jhr,jmin
lfnm=ldir+16
ELSE
WRITE(fname,'(a,a4,a1,3(i2.2),a1,2(i2.2),a5)') &
dir(1:ldir),radar,'.', &
jyr,jmo,jda,'.',jhr,jmin,'.hdf4'
lfnm=ldir+21
ENDIF
RETURN
END SUBROUTINE mkradfnm
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WTRADCOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wtradcol(nx,ny,nz,dmpfmt,iradfmt,hdf4cmpr, & 3,55
rfname,radid,latrad,lonrad,elvrad, &
iyr,imon,iday,ihr,imin,isec,vcpnum,isource, &
refelvmin,refelvmax,refrngmin,refrngmax, &
xsc,ysc,zpsc,gridvel,gridref,gridnyq,gridtim,kcol)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Writes gridded radar data to a file as columns with
! individual lat,lons.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 06/22/95
!
! MODIFICATION HISTORY:
!
! 2000/09/11 (Gene Bassett)
! Use only reflectivity to accept or reject a column (thus allowing one
! to output nids columns without processing velocity data).
!
! 04/29/02 Leilei Wang and Keith Brewster
! Added hdf option, including two new variables in the argument list.
!
! 02/02/04 Keith Brewster
! Modified to replace Time and Nyquist values with missing at points
! where vel and ref are missing for hdf write. Should save compressed
! file space.
!
!-----------------------------------------------------------------------
!
! INPUT:
! dmpfmt file format (1:binary, 2:hdf)
! iradfmt binary format
! hdf4cmpr hdf4 compression level
! rfname radar file name (character*80)
! radid radar id (character*4)
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of radar (m MSL)
! iyr year
! imon month
! iday day
! ihr hour
! imin min
! isec sec
! vcpnum VCP (scan type) number
! isource) source number
! 1= WSR-88D raw
! 2= WSR-88D NIDS
!
! OUTPUT:
! data are written to file
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz
!
INTEGER :: dmpfmt
INTEGER :: iradfmt
INTEGER :: hdf4cmpr
CHARACTER (LEN=256) :: rfname
CHARACTER (LEN=4) :: radid
REAL :: latrad
REAL :: lonrad
REAL :: elvrad
REAL :: refelvmin
REAL :: refelvmax
REAL :: refrngmin
REAL :: refrngmax
INTEGER :: iyr,imon,iday,ihr,imin,isec
INTEGER :: vcpnum
INTEGER :: isource
!
REAL :: xsc(nx)
REAL :: ysc(ny)
REAL :: zpsc(nx,ny,nz)
REAL :: gridref(nx,ny,nz)
REAL :: gridvel(nx,ny,nz)
REAL :: gridnyq(nx,ny,nz)
REAL :: gridtim(nx,ny,nz)
REAL :: kcol(nx,ny)
!
! REAL :: outk(nz)
! REAL :: outhgt(nz)
! REAL :: outref(nz)
! REAL :: outvel(nz)
! REAL :: outnyq(nz)
! REAL :: outtim(nz)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'remap.inc'
!
!-----------------------------------------------------------------------
!
! Radar output descriptors
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: mxradvr=10, nradvr=6
INTEGER :: iradvr(mxradvr)
DATA iradvr /1,2,3,4,5,6,0,0,0,0/
!
!-----------------------------------------------------------------------
!
! Radar output thresholds
!
!-----------------------------------------------------------------------
!
! REAL, PARAMETER :: refmin = -100.0, refmax=100., velmin=-200., velmax=200.
REAL, PARAMETER :: refmin = -5.0, refmax=100., velmin=-200., velmax=200.
REAL, PARAMETER :: misval = -999.0
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iunit,myr,itime
INTEGER :: i,j,k,klev,kntcol
INTEGER :: idummy
INTEGER :: istat,sd_id
INTEGER :: irngmin,irngmax
REAL :: gridlat,gridlon,rdummy
INTEGER(2), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Radar output variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: typelev = 1
INTEGER :: numradcol, numelev
REAL :: xmin, xmax, ymin, ymax
REAL, ALLOCATABLE :: colx(:), coly(:)
REAL, ALLOCATABLE :: collat(:), collon(:)
INTEGER, ALLOCATABLE :: coli(:), colj(:), colk(:,:)
REAL, ALLOCATABLE :: elevsfc(:)
INTEGER, ALLOCATABLE :: colnumelev(:)
REAL, ALLOCATABLE :: radcolhgt(:,:), radcolref(:,:)
REAL, ALLOCATABLE :: radcolvel(:,:), radcolnyq(:,:), radcoltim(:,:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
print *, ' inside wtradcol'
idummy=-999
rdummy=-999.
!
print *, ' dmpfmt = ',dmpfmt
print *, ' hdf4cmpr= ',hdf4cmpr
IF (dmpfmt > 1 .AND. hdf4cmpr > 3) THEN
ALLOCATE (itmp(nx,ny,nz),stat=istat)
CALL check_alloc_status(istat, "wtradcol:itmp")
END IF
myr=1900+iyr
IF(myr < 1960) myr=myr+100
CALL ctim2abss(myr,imon,iday,ihr,imin,isec,itime)
irngmin=nint(refrngmin)
irngmax=nint(refrngmax)
!-----------------------------------------------------------------------
!
! First, go through the 3D reflectivity arrays to find the right
! number of radar columns and the number vertical levels for
! allocating working arrays.
!
! xmin, xmax, ymin, ymax - Valid range for this radar
!
! numradcol - Number of radar columns
! numelev - Maximun number of vertical observation among all radar columns
! typelev - Vertical elevation type
! 0: unknow
! 1: ARPS vertical levels (regular grid)
! 2: Actual radar observation elevation (irregular elevations)
!
!-----------------------------------------------------------------------
xmin = 9.E37
ymin = 9.E37
xmax = -9.E37
ymax = -9.E37
kcol(:,:) = 0
numelev = 0
numradcol= 0
DO j=1,ny
DO i=1,nx
klev=1 ! Note 1 not 0
DO k=2,nz-1 ! Note started from 2 instead of 1
IF((gridref(i,j,k)>refmin .AND. gridref(i,j,k)<refmax) .OR. &
(gridvel(i,j,k)>velmin .AND. gridvel(i,j,k)<velmax)) THEN
klev=klev+1
END IF
END DO
IF(klev > 1) THEN ! Note 1 not 0
xmin = MIN(xmin,xsc(i))
xmax = MAX(xmax,xsc(i))
ymin = MIN(ymin,ysc(j))
ymax = MAX(ymax,ysc(j))
kcol(i,j) = klev
numelev = max(numelev,klev)
numradcol=numradcol+1
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Assign output arrays
!
! radcollat(:) - radar column latitude
! radcollon(:) - radar column longitude
! radnumelev(:) - Number of observation in each column
!
! radcolelev(:,:) - physical height of each observation
! radcolref(:,:) - observations in each radar column at each level
! radcolvel(:,:)
! radcolnyq(:,:)
! radcoltim(:,:)
!
!-----------------------------------------------------------------------
!
ALLOCATE(colx(numradcol), STAT = istat)
ALLOCATE(coly(numradcol), STAT = istat)
ALLOCATE(coli(numradcol), STAT = istat)
ALLOCATE(colj(numradcol), STAT = istat)
ALLOCATE(colk(numelev,numradcol), STAT = istat)
ALLOCATE(collat(numradcol), STAT = istat)
ALLOCATE(collon(numradcol), STAT = istat)
ALLOCATE(elevsfc(numradcol), STAT = istat)
ALLOCATE(colnumelev(numradcol), STAT = istat)
colx(:) = misval
coly(:) = misval
coli(:) = 0
colj(:) = 0
colk(:,:) = 0
elevsfc(:) = misval
collat(:) = misval
collon(:) = misval
colnumelev(:) = 0
ALLOCATE(radcolhgt(numelev,numradcol), STAT = istat)
ALLOCATE(radcolref(numelev,numradcol), STAT = istat)
ALLOCATE(radcolvel(numelev,numradcol), STAT = istat)
ALLOCATE(radcolnyq(numelev,numradcol), STAT = istat)
ALLOCATE(radcoltim(numelev,numradcol), STAT = istat)
radcolhgt(:,:) = misval
radcolref(:,:) = misval
radcolvel(:,:) = misval
radcolnyq(:,:) = misval
radcoltim(:,:) = misval
kntcol = 1
DO j=1,ny
DO i=1,nx
IF (kcol(i,j) > 0) THEN
klev=0
DO k=1,nz-1
IF((gridref(i,j,k)>refmin .AND. gridref(i,j,k)<refmax) .OR. &
(gridvel(i,j,k)>velmin .AND. gridvel(i,j,k)<velmax))THEN
klev=klev+1
colk(klev,kntcol) = k
radcolhgt(klev,kntcol) = zpsc(i,j,k)
radcolref(klev,kntcol) = gridref(i,j,k)
radcolvel(klev,kntcol) = gridvel(i,j,k)
radcolnyq(klev,kntcol) = gridnyq(i,j,k)
radcoltim(klev,kntcol) = gridtim(i,j,k)
END IF
END DO
!
!-----------------------------------------------------------------------
!
! If there are data in this column, write them to the file.
!
!-----------------------------------------------------------------------
!
colnumelev(kntcol) = klev
coli(kntcol) = i ! use in binary format only
colj(kntcol) = j
colx(kntcol) = xsc(i)
coly(kntcol) = ysc(j)
elevsfc(kntcol)=0.5*(zpsc(i,j,1)+zpsc(i,j,2)) ! use in binary format only
CALL xytoll(1,1,xsc(i),ysc(j),gridlat,gridlon)
collat(kntcol) = gridlat
collon(kntcol) = gridlon
kntcol=kntcol+1 ! next radcol NO.
END IF
END DO
END DO
!-----------------------------------------------------------------------
!
! Binary format outputs (be compatible with previous version)
!
!-----------------------------------------------------------------------
IF(dmpfmt == 1)THEN
CALL getunit(iunit)
OPEN(iunit,FILE=TRIM(rfname),STATUS='UNKNOWN',FORM='UNFORMATTED')
!
!-----------------------------------------------------------------------
!
! Write radar description variables
!
!-----------------------------------------------------------------------
!
WRITE(iunit) radid
WRITE(iunit) ireftim,itime,vcpnum,isource,idummy, &
idummy,idummy,nx,ny,nz
!
!-----------------------------------------------------------------------
!
! Write grid description variables
! This should provide enough info to verify that the
! proper grid has been chosen. To recreate the grid,
! icluding elevation information,
! the reading program should get a grid-base-file
! named runname.grdbasfil
!
!-----------------------------------------------------------------------
!
idummy=0
rdummy=0.
WRITE(iunit) runname
WRITE(iunit) iradfmt,strhopt,mapproj,irngmin,irngmax, &
! idummy,idummy,idummy,idummy,idummy
typelev,numradcol,numelev,idummy,idummy
WRITE(iunit) dx,dy,dz,dzmin,ctrlat, &
ctrlon,trulat1,trulat2,trulon,sclfct, &
latrad,lonrad,elvrad,refelvmin,refelvmax
WRITE(iunit) nradvr,iradvr
WRITE(iunit) xmin, xmax, ymin, ymax
!
!-----------------------------------------------------------------------
!
! If there are data in this column, write them to the file.
!
!-----------------------------------------------------------------------
!
DO kntcol = 1, numradcol
klev = colnumelev(kntcol)
WRITE(iunit) coli(kntcol),colj(kntcol),colx(kntcol),coly(kntcol), &
collat(kntcol),collon(kntcol),elevsfc(kntcol),klev
WRITE(iunit) (colk(k,kntcol),k=1,klev)
WRITE(iunit) (radcolhgt(k,kntcol),k=1,klev)
WRITE(iunit) (radcolref(k,kntcol),k=1,klev)
WRITE(iunit) (radcolvel(k,kntcol),k=1,klev)
WRITE(iunit) (radcolnyq(k,kntcol),k=1,klev)
WRITE(iunit) (radcoltim(k,kntcol),k=1,klev)
END DO
CLOSE(iunit)
CALL retunit(iunit)
!
!-----------------------------------------------------------------------
!
! HDF4 format
!
!-----------------------------------------------------------------------
!
ELSE
CALL hdfopen(trim(rfname), 2, sd_id)
IF (sd_id < 0) THEN
WRITE (6,'(1x,a)') 'WTRADCOL: ERROR opening ",trim(rfname)," for writing.'
istat = 1
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Write radar description variables
!
!-----------------------------------------------------------------------
!
CALL hdfwrtc(sd_id, 4, 'radid', radid, istat)
CALL hdfwrti(sd_id, 'ireftim', ireftim, istat)
CALL hdfwrti(sd_id, 'itime', itime, istat)
CALL hdfwrti(sd_id, 'vcpnum', vcpnum, istat)
CALL hdfwrti(sd_id, 'isource', isource, istat)
!
!-----------------------------------------------------------------------
!
! Write grid description variables
! This should provide enough info to verify that the
! proper grid has been chosen. To recreate the grid,
! icluding elevation information,
! the reading program should get a grid-base-file
! named runname.grdbasfil
!
!-----------------------------------------------------------------------
!
CALL hdfwrtc(sd_id, 4, 'runname', runname, istat)
CALL hdfwrti(sd_id, 'iradfmt', iradfmt, istat)
CALL hdfwrti(sd_id, 'strhopt', strhopt, istat)
CALL hdfwrti(sd_id, 'mapproj', mapproj, istat)
CALL hdfwrti(sd_id, 'nx', nx, istat)
CALL hdfwrti(sd_id, 'ny', ny, istat)
CALL hdfwrti(sd_id, 'nz', nz, istat)
CALL hdfwrtr(sd_id, 'dx', dx, istat)
CALL hdfwrtr(sd_id, 'dy', dy, istat)
CALL hdfwrtr(sd_id, 'dz', dz, istat)
CALL hdfwrtr(sd_id, 'dzmin', dzmin, istat)
CALL hdfwrtr(sd_id, 'ctrlat', ctrlat, istat)
CALL hdfwrtr(sd_id, 'ctrlon', ctrlon, istat)
CALL hdfwrtr(sd_id, 'trulat1', trulat1, istat)
CALL hdfwrtr(sd_id, 'trulat2', trulat2, istat)
CALL hdfwrtr(sd_id, 'trulon', trulon, istat)
CALL hdfwrtr(sd_id, 'sclfct', sclfct, istat)
CALL hdfwrtr(sd_id, 'latrad', latrad, istat)
CALL hdfwrtr(sd_id, 'lonrad', lonrad, istat)
CALL hdfwrtr(sd_id, 'elvrad', elvrad, istat)
CALL hdfwrti(sd_id, 'irngmin', irngmin, istat)
CALL hdfwrti(sd_id, 'irngmax', irngmax, istat)
CALL hdfwrtr(sd_id, 'refelvmin', refelvmin, istat)
CALL hdfwrtr(sd_id, 'refelvmax', refelvmax, istat)
CALL hdfwrti(sd_id, 'nradvr', nradvr, istat)
CALL hdfwrt1di(iradvr,mxradvr,sd_id,'iradvr', 'iradvr','index')
CALL hdfwrti(sd_id, 'typelev',typelev, istat)
CALL hdfwrtr(sd_id, 'xmin', xmin, istat)
CALL hdfwrtr(sd_id, 'xmax', xmax, istat)
CALL hdfwrtr(sd_id, 'ymin', xmin, istat)
CALL hdfwrtr(sd_id, 'ymax', ymax, istat)
CALL hdfwrti(sd_id, 'numradcol',numradcol,istat)
CALL hdfwrti(sd_id, 'nummaxelv',numelev, istat)
!
!-----------------------------------------------------------------------
!
! For each horizontal grid point form a column of remapped
! data containing the non-missing grid points
!
!-----------------------------------------------------------------------
!
!
CALL hdfwrt1d(collat,numradcol,sd_id,'radcollat', &
'Latitude array for the radar columns','degree')
CALL hdfwrt1d(collon,numradcol,sd_id,'radcollon', &
'Longitude array for the radar columns','degree')
CALL hdfwrt1di(colnumelev,numradcol,sd_id,'numelev', &
'Number of vertical levels for each column','index')
CALL hdfwrt1di(coli,numradcol,sd_id,'radcoli', &
'i-index in the ARPS grid','index')
CALL hdfwrt1di(colj,numradcol,sd_id,'radcolj', &
'j-index in the ARPS grid','index')
CALL hdfwrt2di(colk,numelev,numradcol,sd_id,0,0, &
'radcolk','k-index in the ARPS grid','index')
CALL hdfwrt2d(radcolhgt,numelev,numradcol,sd_id,0,hdf4cmpr, &
'radcolhgt','height','m',itmp)
CALL hdfwrt2d(radcolref,numelev,numradcol,sd_id,0,hdf4cmpr, &
'radcolref','reflectivity','dbz',itmp)
CALL hdfwrt2d(radcolvel,numelev,numradcol,sd_id,0,hdf4cmpr, &
'radcolvel','radial velocity','m/s',itmp)
CALL hdfwrt2d(radcolnyq,numelev,numradcol,sd_id,0,hdf4cmpr, &
'radcolnyq','nyquist velocity','m/s',itmp)
CALL hdfwrt2d(radcoltim,numelev,numradcol,sd_id,0,hdf4cmpr, &
'radcoltim','relative time','second',itmp)
CALL hdfclose(sd_id,istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR on closing file ",trim(rfname), &
" (status",istat,")"
END IF
IF (hdf4cmpr > 3) DEALLOCATE (itmp,stat=istat)
END IF
!
!-----------------------------------------------------------------------
!
! Report on what data were written
!
!-----------------------------------------------------------------------
!
WRITE(6,'(//a,i4.4,i2.2,i2.2,a1,i2.2,a1,i2.2)') &
' Output statistics for time ', &
iyr,imon,iday,' ',ihr,':',imin
WRITE(6,'(a,i6,a,/a,i6,a//)') &
' There were ',numradcol,' columns written ', &
' of a total ',(nx*ny),' possible.'
DEALLOCATE(colx,coly)
DEALLOCATE(collat, collon)
DEALLOCATE(coli, colj, colk)
DEALLOCATE(elevsfc)
DEALLOCATE(colnumelev)
DEALLOCATE(radcolhgt, radcolref)
DEALLOCATE(radcolvel, radcolnyq, radcoltim)
RETURN
END SUBROUTINE wtradcol
!
SUBROUTINE refract(nx,ny,nz,ptprt,pprt,qv,ptbar,pbar,rfrct) 1
IMPLICIT NONE
!
INTEGER :: nx,ny,nz
REAL :: ptprt(nx,ny,nz)
REAL :: pprt(nx,ny,nz)
REAL :: qv(nx,ny,nz)
REAL :: ptbar(nx,ny,nz)
REAL :: pbar(nx,ny,nz)
REAL :: rfrct(nx,ny,nz)
!
! Constants from Doviak and Zrnic', 1st Ed, Eq 2.18.
! After Bean and Dutton, 1966.
!
REAL, PARAMETER :: cd1=0.776 ! K/Pa
REAL, PARAMETER :: cw1=0.716 ! K/Pa
REAL, PARAMETER :: cw2=3700. ! K*K/Pa
!
! Misc local variables
!
INTEGER :: i,j,k
REAL :: pr,tk,pw,pd,tkinv,refr
!
INCLUDE 'phycst.inc'
!
rfrct =0.
!
! Calculate refractivity index from ARPS pressure,
! potential temp and specific humidity.
!
! pr=total pressure, Pa
! tk=temperature, K
! pw=partial pressure of water vapor, Pa
! pd=partial pressure of dry air, Pa
!
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
pr=pbar(i,j,k)+pprt(i,j,k)
tk=(ptbar(i,j,k)+ptprt(i,j,k))*(pr/p0)**rddcp
pw=pr*(qv(i,j,k)/((0.378*qv(i,j,k))+0.622))
pd=pr-pw
tkinv=1./tk
rfrct(i,j,k)=(cd1*pd*tkinv)+(cw1*pw*tkinv)+ &
(cw2*pw*tkinv*tkinv)
END DO
END DO
END DO
RETURN
END SUBROUTINE refract
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE VADVOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE vadvol(maxgate,maxazim,maxelev, &,1
velchek,rngvad,rngwid,minknt, &
kntvgat,kntvazm,kntvelv, &
rngvvol,azmvvol,elvvvol,velvol, &
vadhgt,vaddir,vadspd)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute vad from volume of radar data from 2-D tilts.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
! November, 2003
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
!
! ngate Number of gates in each radial
! nazim Number of radials in each tilt
!
! velchek Threshold value to determine good vs. flagged data
!
! kntvgat Number of gates in each radial (3-D)
! kntvazm Number of radials in each tilt (3-D)
! kntelev Number of elevation angles (tilts)
!
! kntvazm Number of radials in each tilt (3-D)
! kntvelv Number of elevation angles (tilts)
!
! rngvvol Range to gate
! azmvvol Azimuth angle
! elvvvol Elevation angle
! velvol Radial velocity volume.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: maxgate
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: kntvgat(maxazim,maxelev)
INTEGER, INTENT(IN) :: kntvazm(maxelev)
INTEGER, INTENT(IN) :: kntvelv
REAL, INTENT(IN) :: velchek
REAL, INTENT(IN) :: rngvad
REAL, INTENT(IN) :: rngwid
INTEGER, INTENT(IN) :: minknt
REAL, INTENT(IN) :: rngvvol(maxgate,maxelev)
REAL, INTENT(IN) :: azmvvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvvvol(maxazim,maxelev)
REAL, INTENT(IN) :: velvol(maxgate,maxazim,maxelev)
REAL, INTENT(OUT) :: vadhgt(maxelev)
REAL, INTENT(OUT) :: vaddir(maxelev)
REAL, INTENT(OUT) :: vadspd(maxelev)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: misval = -9999.
REAL, PARAMETER :: mischek = -9990.
INTEGER :: igate,jazim,kelev
INTEGER :: irgate,nwidth,igbgn,igend,igendj,kntpt
REAL :: deg2rad,rad2deg,azmrad,sinaz,cosaz,sin2az,cos2az,delrng,vr
REAL :: sum_q0r,sum_q5r,sum_q5i,sum_q4r,sum_q4i,sum_q3r,sum_q3i
REAL :: flknt,twon,sum_elv,elvavg,height,sfcrng
COMPLEX :: q0,q1,q2,q3,q4,q5,ccj_q4,int_coeff,qq_int,qq
REAL :: cf1,cf2,cf3,dd,ff
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
deg2rad=acos(-1.)/180.
rad2deg=1./deg2rad
vadhgt=misval
vaddir=misval
vadspd=misval
DO kelev = 1, kntvelv
delrng=rngvvol(2,kelev)-rngvvol(1,kelev)
nwidth=MIN(NINT((0.5*rngwid)/delrng),1)
irgate=NINT((rngvad-rngvvol(1,kelev))/delrng)
igbgn=MAX(1,(irgate-nwidth))
igend=(irgate+nwidth)
kntpt=0
sum_elv=0.
sum_q0r=0.
sum_q5r=0.
sum_q5i=0.
sum_q4r=0.
sum_q4i=0.
sum_q3r=0.
sum_q3i=0.
!
! Form sums going around in azim
!
DO jazim = 1, kntvazm(kelev)
azmrad=deg2rad*azmvvol(jazim,kelev)
sinaz=sin(azmrad)
cosaz=cos(azmrad)
sin2az=sin(2.0*azmrad)
cos2az=cos(2.0*azmrad)
igendj=MIN(igend,kntvgat(jazim,kelev))
DO igate=igbgn,igendj
IF(velvol(igate,jazim,kelev) > velchek) THEN
vr=velvol(igate,jazim,kelev)
kntpt=kntpt+1
sum_elv=sum_elv+elvvvol(jazim,kelev)
sum_q0r=sum_q0r+vr
sum_q5r=sum_q5r+cos2az
sum_q5i=sum_q5i+sin2az
sum_q4r=sum_q4r+cosaz
sum_q4i=sum_q4i+sinaz
sum_q3r=sum_q3r+vr*cosaz
sum_q3i=sum_q3i+vr*sinaz
END IF
END DO
END DO
!
cf1=misval
cf2=misval
cf3=misval
!
IF(kntpt > minknt) THEN
flknt=float(kntpt)
twon=float(2*kntpt)
elvavg=sum_elv/flknt
q0=CMPLX(sum_q0r/flknt)
q5=CMPLX((sum_q5r/twon),-(sum_q5i/twon))
q4=CMPLX((sum_q4r/twon),(sum_q4i/twon))
q3=CMPLX((sum_q3r/flknt),-(sum_q3i/flknt))
!
! Complex conjugate of q4
!
ccj_q4=CONJG(q4)
qq=q4-(1./(4.*ccj_q4))
!
IF( qq /= 0.) THEN
q2=(ccj_q4-(q5/(2.*ccj_q4)))/qq
q1=(q0-(q3/(2.*ccj_q4)))/qq
qq_int=1.-(CABS(q2)*CABS(q2))
IF( qq_int /= 0.) THEN
int_coeff=(q1-q2*CONJG(q1))/qq_int
cf3=IMAG(int_coeff)
cf2=REAL(int_coeff)
cf1=REAL(q0)-(2.*REAL(int_coeff*q4))
END IF
END IF
END IF
IF( cf1 > mischek ) THEN
ff=sqrt(cf2*cf2 + cf3*cf3)/cos(deg2rad*elvavg)
dd=180.-rad2deg*(atan2(cf3,cf2))
IF(dd > 360.) dd = dd-360.
CALL beamhgt(elvavg,rngvad,height,sfcrng)
print *, ' elv:',elvavg,' hgt:',height,' dd:',dd,' ff:',ff
vadhgt(kelev)=height
vaddir(kelev)=dd
vadspd(kelev)=ff
END IF
END DO
END SUBROUTINE vadvol
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKVADFNM ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkvadfnm(dir,ldir,radar,iyr,imo,ida,ihr,imin,isec, &,2
fname,lfnm)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=256) :: dir
INTEGER :: ldir
CHARACTER (LEN=4) :: radar
INTEGER :: iyr
INTEGER :: imo
INTEGER :: ida
INTEGER :: ihr
INTEGER :: imin
INTEGER :: isec
!
!-----------------------------------------------------------------------
!
! Output arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=256) :: fname
INTEGER :: lfnm
!
INTEGER :: rdtime,jyr,jmo,jda,jhr,jmin,jsec
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Round to nearest minute
!
!-----------------------------------------------------------------------
!
jyr=iyr+1900
IF(iyr < 50) jyr=jyr+100
CALL ctim2abss( jyr,imo,ida,ihr,imin,isec, rdtime )
IF(isec >= 30) rdtime=rdtime+(60-isec)
CALL abss2ctim( rdtime, jyr, jmo, jda, jhr, jmin, jsec )
jyr=MOD(jyr,100)
WRITE(fname,'(a,a4,a1,3(i2.2),2(i2.2),a4)') &
dir(1:ldir),radar,'.', &
jyr,jmo,jda,jhr,jmin,'.vad'
lfnm=ldir+19
RETURN
END SUBROUTINE mkvadfnm
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WTVADWIND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wtvadprf(maxelev,vfname,radar,radlat,radlon,radelv, &,1
vadhgt,vaddir,vadspd)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write VAD wind profile to file to be read by ADAS.
! Format of wind profiler data is used.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
! November, 2003
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
IMPLICIT NONE
INTEGER :: maxelev
CHARACTER (LEN=256) :: vfname
CHARACTER (LEN=5) :: radar
REAL :: radlat
REAL :: radlon
REAL :: radelv
REAL :: vadhgt(maxelev)
REAL :: vaddir(maxelev)
REAL :: vadspd(maxelev)
!
! Misc local variables
!
INTEGER :: kelev
INTEGER :: kntvalid,numsta,iunit
!
! Determine if there are any valid data in this VAD profile
!
numsta=88888
kntvalid=0
DO kelev=1,maxelev
IF(vaddir(kelev) <= 360. .AND. vaddir(kelev) >= 0.) &
kntvalid=kntvalid+1
END DO
!
WRITE(6,'(i6,a)') kntvalid,' valid heights in VAD profile'
IF ( kntvalid > 0 ) THEN
WRITE(6,'(a,a,a)') ' Opening ',TRIM(vfname),' for writing'
CALL GETUNIT(iunit)
OPEN(iunit,file=TRIM(vfname),status='unknown')
WRITE(iunit,'(i12,i12,f11.4,f15.4,f15.0,6x,a4)') &
numsta,kntvalid,radlat,radlon,radelv,radar(1:4)
DO kelev=1,maxelev
IF(vaddir(kelev) <= 360. .AND. vaddir(kelev) >= 0.) THEN
WRITE(iunit,'(f10.0,f10.0,f10.1)') &
vadhgt(kelev),vaddir(kelev),vadspd(kelev)
END IF
END DO
CLOSE(iunit)
CALL RETUNIT(iunit)
ELSE
WRITE(6,'(//a//)') ' No valid VAD data. Skipping VAD write.'
END IF
RETURN
END SUBROUTINE wtvadprf