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