!
!##################################################################
!##################################################################
!###### ######
!###### 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, & 1,2
fname,lfnm)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=120) :: 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=100) :: 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,41
rfname,radid,latrad,lonrad,elvrad, &
iyr,imon,iday,ihr,imin,isec,vcpnum,isource, &
refelvmin,refelvmax,refrngmin,refrngmax, &
xsc,ysc,zpsc,gridvel,gridref,gridnyq,gridtim, &
outk,outhgt,outref,outvel,outnyq,outtim)
!
!-----------------------------------------------------------------------
!
! 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=100) :: 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 :: 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 :: mxradvr,nradvr
PARAMETER(mxradvr=10,nradvr=6)
INTEGER :: iradvr(mxradvr)
DATA iradvr /1,2,3,4,5,6,0,0,0,0/
!
!-----------------------------------------------------------------------
!
! Radar output thresholds
!
!-----------------------------------------------------------------------
!
REAL :: refmin,refmax,velmin,velmax
PARAMETER(refmin=-5.0, refmax=100., &
velmin=-200.,velmax=200.)
REAL :: misval
PARAMETER(misval=-999.0)
!
!-----------------------------------------------------------------------
!
! Radar output variables
!
!-----------------------------------------------------------------------
!
REAL :: outarray(nz,6)
REAL :: grdlatc(nx,ny)
REAL :: grdlonc(nx,ny)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iunit,myr,itime
INTEGER :: i,j,k,klev,kk,kntcol,nn
INTEGER :: idummy
INTEGER :: istat,sd_id
INTEGER :: irngmin,irngmax
REAL :: gridlat,gridlon,elev,rdummy
INTEGER(2), allocatable :: itmp(:,:,:) ! Temporary array
REAL, allocatable :: hmax(:), hmin(:) ! Temporary array
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! 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)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmp, returning"
RETURN
END IF
ALLOCATE (hmax(nz),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmax, returning"
RETURN
END IF
ALLOCATE (hmin(nz),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmin, returning"
RETURN
END IF
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)
IF(dmpfmt == 1)THEN
CALL getunit(iunit)
!
!-----------------------------------------------------------------------
!
! Open file for output
!
!-----------------------------------------------------------------------
!
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,idummy,idummy,idummy
!
!-----------------------------------------------------------------------
!
! 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
WRITE(iunit) dx,dy,dz,dzmin,ctrlat, &
ctrlon,trulat1,trulat2,trulon,sclfct, &
latrad,lonrad,elvrad,refelvmin,refelvmax
WRITE(iunit) nradvr,iradvr
ELSE !HDF4 format
!
!-----------------------------------------------------------------------
!
! Open file for output
!
!-----------------------------------------------------------------------
!
CALL hdfopen(trim(rfname), 2, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "WTRADCOL: ERROR opening ", &
trim(rfname)," for writing."
istat = 1
STOP
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 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 hdfwrt1d(iradvr,mxradvr,sd_id,'iradvr', 'iradvr','')
ENDIF
!
!-----------------------------------------------------------------------
!
! For each horizontal grid point form a column of remapped
! data containing the non-missing grid points
!
!-----------------------------------------------------------------------
!
IF(dmpfmt==1)THEN
kntcol=0
DO j=1,ny
DO i=1,nx
DO k=1,nz
outk(k)=misval
outhgt(k)=misval
outref(k)=misval
outvel(k)=misval
outnyq(k)=misval
outtim(k)=misval
END DO
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
outk(klev)=FLOAT(k)
outhgt(klev)=zpsc(i,j,k)
outref(klev)=gridref(i,j,k)
outvel(klev)=gridvel(i,j,k)
outnyq(klev)=gridnyq(i,j,k)
outtim(klev)=gridtim(i,j,k)
END IF
END DO
!
!-----------------------------------------------------------------------
!
! If there are data in this column, write them to the file.
!
!-----------------------------------------------------------------------
!
IF(klev > 0) THEN
kntcol=kntcol+1
CALL xytoll(1,1,xsc(i),ysc(j),gridlat,gridlon)
elev=0.5*(zpsc(i,j,1)+zpsc(i,j,2))
WRITE(iunit) i,j,xsc(i),ysc(j), &
gridlat,gridlon,elev,klev
WRITE(iunit) (outk(k),k=1,klev)
WRITE(iunit) (outhgt(k),k=1,klev)
WRITE(iunit) (outref(k),k=1,klev)
WRITE(iunit) (outvel(k),k=1,klev)
WRITE(iunit) (outnyq(k),k=1,klev)
WRITE(iunit) (outtim(k),k=1,klev)
END IF
END DO
END DO
ELSE !HDF4 format
!
! Improve compression by setting Time and Nyquist values to misval
! where corresponding velocity and reflectivity values are missing.
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
IF(gridvel(i,j,k) < velmin .OR. gridvel(i,j,k) > velmax) THEN
gridnyq(i,j,k)=misval
IF(gridref(i,j,k) < refmin .OR. gridref(i,j,k) > refmax) &
gridtim(i,j,k)=misval
END IF
END DO
END DO
END DO
!
CALL xytoll(nx,ny,xsc,ysc,grdlatc,grdlonc)
CALL hdfwrt2d(grdlatc,nx,ny,sd_id,0,hdf4cmpr, &
'grdlatc','grid latitude','degree',itmp)
CALL hdfwrt2d(grdlonc,nx,ny,sd_id,0,hdf4cmpr, &
'grdlonc','grid lontitude','degree',itmp)
! write(19,'(10f10.1)') gridref
CALL hdfwrt3d(zpsc,nx,ny,nz,sd_id,0,hdf4cmpr, &
'hgtrad','height','m', &
itmp,hmax,hmin)
CALL hdfwrt3d(gridref,nx,ny,nz,sd_id,0,hdf4cmpr, &
'gridref','reflectivity','dbz', &
itmp,hmax,hmin)
CALL hdfwrt3d(gridvel,nx,ny,nz,sd_id,0,hdf4cmpr, &
'gridvel','radial velocity','m/s', &
itmp,hmax,hmin)
CALL hdfwrt3d(gridnyq,nx,ny,nz,sd_id,0,hdf4cmpr, &
'gridnyq','nyquist velocity','m/s', &
itmp,hmax,hmin)
CALL hdfwrt3d(gridtim,nx,ny,nz,sd_id,0,hdf4cmpr, &
'gridtim','time','', &
itmp,hmax,hmin)
!
! Although all columns are actually written, find the non-missing
! data columns for statistics reporting purposes.
!
kntcol=0
DO j=1,ny
DO i=1,nx
DO k=1,nz
IF(gridref(i,j,k)>refmin .AND. gridref(i,j,k)<refmax) THEN
kntcol=kntcol+1
EXIT
END IF
END DO
END DO
END DO
ENDIF
!
IF(dmpfmt==1)THEN
CLOSE(iunit)
CALL retunit(iunit)
ELSE
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) THEN
DEALLOCATE (itmp,stat=istat)
DEALLOCATE (hmax,stat=istat)
DEALLOCATE (hmin,stat=istat)
END IF
ENDIF
!
!-----------------------------------------------------------------------
!
! 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 ',kntcol,' columns written ', &
' of a total ',(nx*ny),' possible.'
!
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=120) :: dir
INTEGER :: ldir
CHARACTER (LEN=4) :: radar
INTEGER :: iyr
INTEGER :: imo
INTEGER :: ida
INTEGER :: ihr
INTEGER :: imin
INTEGER :: isec
!
!-----------------------------------------------------------------------
!
! Output arguments
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=100) :: 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=100) :: 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