! !################################################################## !################################################################## !###### ###### !###### Advanced Regional Prediction System (ARPS) ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! PROGRAM PLT_RADAR_TILT,99 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! Plot radar reflectivity and radial velocity in a tilt and ! vertical cross-sections ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Ming Xue ! 03/19/02. ! ! MODIFICATION HISTORY: ! 06/12/2005 Keith Brewster, CAPS ! Added reading of radar emulator volume files (radfmt=2) ! ! 11/15/2005 Keith Brewster, CAPS ! Updated for changes to content of emulator volume files. ! ! 11/23/2005 Keith Brewster, CAPS ! Restructured code to use considerably less memory. ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, PARAMETER :: mxradfile= 100 INTEGER, PARAMETER :: mxmapfile= 10 INTEGER, PARAMETER :: nref_tilts_max=14 INTEGER, PARAMETER :: nvel_tilts_max=14 INTEGER, PARAMETER :: n_elevation_max=14 INTEGER, PARAMETER :: n_azimuth_max=361 ! !----------------------------------------------------------------------- ! ! Namelist ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: filepath INTEGER :: radfmt INTEGER :: numfiles CHARACTER (LEN=256) :: files(mxradfile) NAMELIST /input_data/ filepath,radfmt,numfiles,files CHARACTER(LEN=4) :: radar_symbol NAMELIST /radar_index/ radar_symbol INTEGER :: nx INTEGER :: ny REAL :: dx REAL :: dy REAL :: ctrlat REAL :: ctrlon NAMELIST /grid/ nx,ny,dx,dy,ctrlat,ctrlon INTEGER :: mapproj REAL :: trulat1 REAL :: trulat2 REAL :: trulon REAL :: sclfct NAMELIST /projection/ mapproj,trulat1,trulat2,trulon,sclfct INTEGER :: overlayrefcontour INTEGER :: n_elevation_ref INTEGER :: n_elevation_vel INTEGER :: n_azimuth_ref INTEGER :: n_azimuth_vel REAL :: elevations_ref(n_elevation_max) REAL :: elevations_vel(n_elevation_max) REAL :: azimuth_ref(n_azimuth_max) REAL :: azimuth_vel(n_azimuth_max) REAL :: h_range_max,h_range_max_e,h_range_max_w,h_range_max_s,h_range_max_n REAL :: z_range_max REAL :: elev_error_range INTEGER :: coltab_ref INTEGER :: ibgncol_ref INTEGER :: iendcol_ref CHARACTER(LEN=80) :: coltabfn_ref INTEGER :: coltab_vel INTEGER :: ibgncol_vel INTEGER :: iendcol_vel CHARACTER(LEN=80) :: coltabfn_vel NAMELIST /plotting/ h_range_max_e,h_range_max_w,h_range_max_s, & h_range_max_n,z_range_max,elev_error_range, & overlayrefcontour, & n_elevation_ref,elevations_ref, & n_elevation_vel,elevations_vel, & n_azimuth_ref,azimuth_ref, & n_azimuth_vel,azimuth_vel, & coltab_ref,ibgncol_ref,iendcol_ref,coltabfn_ref, & coltab_vel,ibgncol_vel,iendcol_vel,coltabfn_vel INTEGER :: ovrmap INTEGER :: mapgrid INTEGER :: nmapfile CHARACTER (LEN=256) :: mapfile(mxmapfile) REAL :: latgrid REAL :: longrid NAMELIST /map_plot/ ovrmap,mapgrid,latgrid,longrid,nmapfile,mapfile INTEGER :: lmapfile REAL :: xorig, yorig ! !----------------------------------------------------------------------- ! ! plot variables ! !----------------------------------------------------------------------- ! REAL :: cl(100) REAL, ALLOCATABLE :: xrplt(:,:) REAL, ALLOCATABLE :: yrplt(:,:) REAL, ALLOCATABLE :: xvplt(:,:) REAL, ALLOCATABLE :: yvplt(:,:) REAL, ALLOCATABLE :: rplt(:,:) REAL, ALLOCATABLE :: rpltm(:,:) REAL, ALLOCATABLE :: zplt(:,:) REAL, ALLOCATABLE :: xrw(:) REAL, ALLOCATABLE :: yrw(:) REAL, ALLOCATABLE :: xvw(:) REAL, ALLOCATABLE :: yvw(:) REAL, ALLOCATABLE :: elevmean(:) INTEGER, ALLOCATABLE :: irwxy(:,:) INTEGER, ALLOCATABLE :: ivwxy(:,:) INTEGER, ALLOCATABLE :: iwrz(:,:) INTEGER :: ncl INTEGER :: mode ! !----------------------------------------------------------------------- ! ! INPUT RADAR PARAMETER AND OBSERVATION ! !----------------------------------------------------------------------- ! REAL :: radar_lat REAL :: radar_lon REAL :: radar_elevation INTEGER :: iyear,imon,iday,ihour,imin,isec CHARACTER(LEN=256) :: datafile CHARACTER(LEN=80) :: runname INTEGER :: lfnkey,abstsec INTEGER :: wrtuaref,wrtuavel,wrtvort,idummy,vcp REAL :: beamwid,rmisval,rngfval,curtim REAL, ALLOCATABLE :: elvref(:) REAL, ALLOCATABLE :: elvvel(:) REAL, ALLOCATABLE :: reftilt(:,:) REAL, ALLOCATABLE :: veltilt(:,:) REAL, ALLOCATABLE :: datarz(:,:) INTEGER :: maxvelgate INTEGER :: maxrefgate INTEGER :: maxazim INTEGER :: maxelev INTEGER, ALLOCATABLE :: itimvol(:) REAL, ALLOCATABLE :: rngvvol(:,:) REAL, ALLOCATABLE :: azmvvol(:,:) REAL, ALLOCATABLE :: elvvvol(:,:) REAL, ALLOCATABLE :: velvol(:,:,:) REAL, ALLOCATABLE :: rngrvol(:,:) REAL, ALLOCATABLE :: azmrvol(:,:) REAL, ALLOCATABLE :: elvrvol(:,:) REAL, ALLOCATABLE :: refvol(:,:,:) REAL, ALLOCATABLE :: vnyquist(:) REAL, ALLOCATABLE :: uarefvol(:,:,:) REAL, ALLOCATABLE :: uavelvol(:,:,:) ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: ifiles CHARACTER(LEN=100) :: string CHARACTER(LEN=100) :: stringtime real :: radmax ,rad INTEGER :: j_razim_start,nref_tilts,nvel_tilts INTEGER :: iazimuth real :: razim_start REAL :: elev_to_plot, xscale, yscale INTEGER :: kref,kvel,ktilt INTEGER :: kref_tilts,kvel_tilts REAL :: refmax, refmin, reflmin, reflmax REAL :: velmax, velmin REAL :: eleva, deg2rad REAL :: xradar, yradar, tem0,tem1,tem2,tem3,tem4 LOGICAL :: fexist INTEGER :: zxplot_initialized = 0 REAL :: missing_value = -999.0 INTEGER :: iunit INTEGER :: i,j,k,jj CHARACTER(LEN=256) :: psfilename INTEGER :: ilen REAL :: factorx,factory REAL :: delazm,dazmin REAL :: elvsum,data_azim,last_azim,azmin,xconst,yconst INTEGER :: elvknt,jidx,jazim,nazim ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! deg2rad=acos(-1.0)/180.0 !----------------------------------------------------------------------- ! Input data file name and grid !----------------------------------------------------------------------- READ(5,input_data,ERR=100) WRITE(6,'(a)') 'Namelist input_data was sucessfully read in.' READ(5,radar_index,ERR=100) WRITE(6,'(a)') 'Namelist radar_index was sucessfully read in.' READ(5,grid,ERR=100) WRITE(6,'(a)') 'Namelist grid was sucessfully read in.' !----------------------------------------------------------------------- ! Input map projection parameters !----------------------------------------------------------------------- sclfct = 1.0 READ (5,projection,ERR=100) WRITE(6,'(a)') 'Namelist projection was sucessfully read in.' READ(5,plotting,ERR=100) WRITE(6,'(a)') 'Namelist plotting sucessfully read in.' READ(5,map_plot,ERR=100) WRITE(6,'(a)') 'Namelist map_plot sucessfully read in.' ! h_range_max = (h_range_max_e+h_range_max_w)/2.0 !----------------------------------------------------------------------- ! Read in radar data volume !----------------------------------------------------------------------- DO ifiles=1,numfiles ! ! read in data on a tilt ! WRITE(datafile,'(a,a,a)') TRIM(filepath),'/',TRIM(files(ifiles)) INQUIRE(FILE=TRIM(datafile),EXIST = fexist) IF( .not. fexist) THEN WRITE(6,'(a,a,a)') 'Radar data file ',TRIM(datafile),' not found.' CYCLE END IF WRITE(6,'(a,i4)') ' Radar format = ',radfmt iunit=27 IF( radfmt == 1 ) THEN WRITE(6,'(a,a)') ' Opening: ',TRIM(datafile) OPEN(iunit,file=TRIM(datafile),form='unformatted',status='old') READ(iunit) maxvelgate,maxrefgate,maxazim,maxelev READ(iunit) radar_elevation,radar_lat, radar_lon READ(iunit) iyear,imon,iday,ihour,imin,isec write(6,'(3i10)') maxvelgate,maxrefgate,maxazim,maxelev write(6,'(3f10.2)') radar_elevation,radar_lat,radar_lon write(6,'(6i6)') iyear,imon,iday,ihour,imin,isec ALLOCATE( rngrvol(maxrefgate,maxelev) ) ALLOCATE( azmrvol(maxazim,maxelev) ) ALLOCATE( elvrvol(maxazim,maxelev) ) ALLOCATE( refvol(maxrefgate,maxazim,maxelev) ) ALLOCATE( rngvvol(maxvelgate,maxelev) ) ALLOCATE( azmvvol(maxazim,maxelev) ) ALLOCATE( elvvvol(maxazim,maxelev) ) ALLOCATE( velvol(maxvelgate,maxazim,maxelev) ) READ(iunit) rngrvol READ(iunit) azmrvol READ(iunit) elvrvol READ(iunit) refvol READ(iunit) rngvvol READ(iunit) azmvvol READ(iunit) elvvvol READ(iunit) velvol close(iunit) psfilename=TRIM(files(ifiles)) ilen=LEN(TRIM(psfilename)) iday=ICHAR(psfilename(ilen-12:ilen-12))-ICHAR("0") iday=iday*10+(ICHAR(psfilename(ilen-11:ilen-11))-ICHAR("0")) ihour=ICHAR(psfilename(ilen-9:ilen-9))-ICHAR("0") ihour=ihour*10+(ICHAR(psfilename(ilen-8:ilen-8))-ICHAR("0")) imin=ICHAR(psfilename(ilen-7:ilen-7))-ICHAR("0") imin=imin*10+(ICHAR(psfilename(ilen-6:ilen-6))-ICHAR("0")) ELSE IF (radfmt == 2) THEN ! radar emulator binary file WRITE(6,'(a,a)') ' Opening: ',TRIM(datafile) OPEN(iunit,FILE=TRIM(datafile),FORM='unformatted',STATUS='old') READ(iunit) runname lfnkey=80 CALL gtlfnkey(runname,lfnkey) WRITE(6,'(a,a)') ' Runname:',runname(1:lfnkey) READ(iunit) radar_symbol WRITE(6,'(a,a)') ' Radar name:',radar_symbol READ(iunit) maxvelgate,maxrefgate,maxazim,maxelev WRITE(6,'(a,i6,a,i6,a,i6)') & ' Ngate : ',maxvelgate,' Nazim: ',maxazim,' maxelev: ',maxelev READ(iunit) iyear,imon,iday,ihour,imin,isec WRITE(6,'(2(a,i2.2),a,i4.4,3(a,i2.2))') & ' Date: ',imon,'/',iday,'/',iyear,' Time: ',ihour,':',imin,':',isec READ(iunit) radar_elevation,radar_lat,radar_lon WRITE(6,'(a,f9.2,a,f9.2,a,f9.2)') & ' Latitude:',radar_lat,' Longitude:',radar_lon, & ' Elevation:',radar_elevation READ(iunit) vcp,beamwid,rmisval,rngfval,curtim WRITE(6,'(a,i5,a,f9.2,a,f10.1,a,f10.1)') & ' VCP:',vcp,' Beamwidth:',beamwid, & ' MissVal: ',rmisval,' RngfVal :',rngfval ! CALL ctim2abss(iyear,imon,iday,ihour,imin,isec,abstsec) READ(iunit) wrtuaref,wrtuavel,wrtvort, & idummy,idummy,idummy,idummy,idummy WRITE(6,'(a)') ' Allocating volume memory ' ALLOCATE(itimvol(maxelev)) ALLOCATE(rngrvol(maxrefgate,maxelev)) ALLOCATE(azmrvol(maxazim,maxelev)) ALLOCATE(elvrvol(maxazim,maxelev)) ALLOCATE(refvol(maxrefgate,maxazim,maxelev)) IF(wrtuaref /= 0) ALLOCATE(uarefvol(maxrefgate,maxazim,maxelev)) ALLOCATE(vnyquist(maxelev)) ALLOCATE(rngvvol(maxvelgate,maxelev)) ALLOCATE(azmvvol(maxazim,maxelev)) ALLOCATE(elvvvol(maxazim,maxelev)) ALLOCATE(velvol(maxvelgate,maxazim,maxelev)) IF(wrtuavel /= 0) ALLOCATE(uavelvol(maxvelgate,maxazim,maxelev)) WRITE(6,'(a)') ' Reading Time Array' READ(iunit) itimvol WRITE(6,'(a)') ' Reading Reflectivity Arrays' READ(iunit) rngrvol READ(iunit) azmrvol READ(iunit) elvrvol READ(iunit) refvol IF( wrtuaref /= 0 ) READ(iunit) uarefvol WRITE(6,'(a)') ' Reading Velocity Arrays' READ(iunit) vnyquist READ(iunit) rngvvol READ(iunit) azmvvol READ(iunit) elvvvol READ(iunit) velvol IF( wrtuavel /= 0 ) READ(iunit) uavelvol ! WRITE(6,'(a)') ' Reading successfully completed' CLOSE(iunit) ELSE IF (radfmt == 4) THEN ! CASA-formatted NetCDF file WRITE(6,'(a,i6,a)') ' Radar format number ',radfmt,' not yet implemented' STOP ELSE WRITE(6,'(a,i6,a)') ' Radar format number ',radfmt,' not supported.' STOP END IF IF ( mapproj == 0 ) THEN trulat1 = radar_lat trulat2 = radar_lat trulon = radar_lon END IF !----------------------------------------------------------------------- ! plot set-up !----------------------------------------------------------------------- print *, ' Begin plotting' ALLOCATE(elvref(maxelev)) ALLOCATE(elvvel(maxelev)) factorx=0.3 factory=0.3 IF(h_range_max_w + h_range_max_e > h_range_max_s + h_range_max_n ) THEN h_range_max=(h_range_max_e+h_range_max_w)/2 factory=factorx*(h_range_max_s + h_range_max_n)/(h_range_max_w + h_range_max_e) ELSE factorx=factory*(h_range_max_w + h_range_max_e)/(h_range_max_s + h_range_max_n) h_range_max=(h_range_max_n+h_range_max_s)/2 END IF IF( n_elevation_ref > 0 .or. n_azimuth_ref > 0 .or. & n_elevation_vel > 0 .or. n_azimuth_vel > 0 ) THEN IF( zxplot_initialized == 0 ) THEN CALL xpsfn('zxout'//'.ps', 5) ! PS output will be called zxout.ps CALL xpaprlnth( 1.0 ) CALL xdevic ! begin graphic zxplot_initialized = 1 END IF END IF elvref=-99. DO ktilt=1,maxelev elvsum=0. elvknt=0 DO j=1,maxazim IF(elvrvol(j,ktilt) > -9.) THEN elvknt=elvknt+1 elvsum=elvsum+elvrvol(j,ktilt) END IF END DO IF(elvknt > 0) THEN elvref(ktilt)=elvsum/float(elvknt) nref_tilts=ktilt ELSE elvref(ktilt)=-99. END IF WRITE(6,'(a,i4,a,f10.2)') ' Elevation(',ktilt,') = ',elvref(ktilt) END DO WRITE(6,'(a,i4)') ' nref_tilts=',nref_tilts elvvel=-99. DO ktilt=1,maxelev elvsum=0. elvknt=0 DO j=1,maxazim IF(elvvvol(j,ktilt) > -9.) THEN elvknt=elvknt+1 elvsum=elvsum+elvrvol(j,ktilt) END IF END DO IF(elvknt > 0) THEN elvvel(ktilt)=elvsum/float(elvknt) nref_tilts=ktilt ELSE elvvel(ktilt)=-99. END IF WRITE(6,'(a,i4,a,f10.2)') ' Elevation(',ktilt,') = ',elvvel(ktilt) END DO WRITE(6,'(a,i4)') ' nvel_tilts=',nvel_tilts !----------------------------------------------------------------------- ! plot reflectivity variables !----------------------------------------------------------------------- IF(n_elevation_ref > 0) THEN WRITE(6,'(a)') ' Allocating xrplt,yrplt,reftilt ...' ALLOCATE( xrplt(maxrefgate,(maxazim+1)) ) ALLOCATE( yrplt(maxrefgate,(maxazim+1)) ) ALLOCATE( reftilt(maxrefgate,(maxazim+1)) ) ALLOCATE( irwxy(maxrefgate,(maxazim+1)) ) ALLOCATE( xrw(8*maxrefgate) ) ALLOCATE( yrw(8*maxrefgate) ) WRITE(6,'(a)') ' Reflectivity tilt plotting...' DO ktilt = 1,n_elevation_ref elev_to_plot = elevations_ref(ktilt) kref = 0 DO k = 1,maxelev IF( abs( elvref(k)-elev_to_plot ) & < elev_error_range ) THEN kref = k EXIT END IF END DO IF(kref == 0 ) CYCLE ! no elevation angle ! within error range was found refmin = 1000.0 refmax =-1000.0 velmin = 1000.0 refmax =-1000.0 eleva=elvref(kref) last_azim=-999. jazim=0 DO j=1,maxazim jidx=j IF(jidx > maxazim) jidx=jidx-maxazim data_azim=azmrvol(jidx,kref) IF(data_azim >= 0. .AND. (abs(data_azim-last_azim) > 1.0E-04)) THEN jazim=jazim+1 last_azim=data_azim xconst=cos( (90.-data_azim)*deg2rad ) yconst=sin( (90.-data_azim)*deg2rad ) reflmax=-999. reflmin=999. DO i=1,maxrefgate rad = 0.001 * rngrvol(i,kref) xrplt(i,jazim)= rad * xconst yrplt(i,jazim)= rad * yconst reftilt(i,jazim)=max(refvol(i,jidx,kref),-30.) IF(refvol(i,jidx,kref) > missing_value ) THEN refmax=max(refvol(i,jidx,kref),refmax) refmin=min(refvol(i,jidx,kref),refmin) reflmax=max(refvol(i,jidx,kref),reflmax) reflmin=min(refvol(i,jidx,kref),reflmin) END IF END DO print *, ' Azimuth,reflmin,reflmax: ',data_azim,reflmin,reflmax END IF END DO nazim=jazim WRITE(6,'(a,2f10.2)') ' Refmin, Refmax on tilt= ', refmin, refmax WRITE(6,'(a,i5)') 'N unique azimuths: ',nazim ! !----------------------------------------------------------------------- ! Define plotting space. ! CALL xpspac(0.5-factorx,0.5+factorx,0.5-factory,0.5+factory) CALL xmap(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) CALL xlabmask(1) IF( coltab_ref == -1 ) CALL xstctfn(coltabfn_ref) CALL xsetclrs(coltab_ref) ! Use internally defined color map number 6 CALL xcolor(1) ! Set color to black CALL xcfont(3) ! call xaxfmt('(I4)') ! ! Set tentative color interval and limit on the min/max number of contours ! cl(1)=0.0 cl(2)=5.0 ! Set tentative contour interval CALL xnctrs(2,100) ! Set lower and upper limits ! on the number of contours allowed ! It's used by both xcolfil and xclimit CALL xnctrs(5,20) ! Set lower and upper limits mode=2 call xctrbadv(1) ! Turn on bad value checking ! in contouring routine call xbadval( missing_value ) ! Specify bad value flag ! ! plot a color filled contour field for positive values ! CALL xctrclr(ibgncol_ref,iendcol_ref) ! Use colors between 10 and 39 inclusive CALL xctrlim(5.0, 75.0) ! Only contours above zero are plotted CALL xwindw(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) WRITE(6,'(a,f10.2)') & ' Plotting colfil reflectivity field at elevation angle = ', eleva CALL xcolfil(reftilt,xrplt,yrplt,irwxy,xrw,yrw, & maxrefgate,maxrefgate,nazim, cl, ncl,mode) IF ( overlayrefcontour == 1 ) THEN call XCLTYP(0) CALL xctrclr(14,15) ! Use colors between 10 and 39 inclusive mode=3 ncl=1 cl(1)=15.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=30.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=45.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=55.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) END IF CALL xwdwof CALL xaxsca(-h_range_max_w, h_range_max_e, 0.1*h_range_max, & -h_range_max_s, h_range_max_n, 0.1*h_range_max) WRITE(string,'(a,f5.2)') 'Reflectivity at Elevation Angle=',eleva CALL xchsiz(0.025*2*h_range_max ) ! Set character size CALL xcharc((-h_range_max_w+h_range_max_e)*0.5, & h_range_max_n+2*h_range_max*0.04, trim(string) ) WRITE(stringtime,'(a,i2.2,a,i2.2,a,i4,a,i2.2,a,i2.2)') 'Scan Time:', & imon,'/',iday,'/',iyear, ' ',ihour,':',imin CALL xchsiz(0.025*2*h_range_max ) ! Set character size CALL xcharc((-h_range_max_w+h_range_max_e)*0.5, & h_range_max_n+2*h_range_max*0.01, trim(stringtime) ) CALL xchsiz(0.02* 2*h_range_max ) ! Set character size CALL xcpalet(1) ! Plot horizontal color palette xscale = h_range_max yscale = h_range_max xorig = -(h_range_max_e+h_range_max_w)/2.0 yorig = -(h_range_max_s+h_range_max_n)/2.0 xradar = 0.0 yradar = 0.0 CALL xthick(3) CALL xpenup( xradar-xscale*0.02, yradar ) CALL xpendn( xradar+xscale*0.02, yradar ) CALL xpenup( xradar, yradar-yscale*0.02 ) CALL xpendn( xradar, yradar+yscale*0.02 ) call xthick(1) CALL xchsiz(0.02* 2*h_range_max ) ! Set character size CALL xcharc(xradar, yradar-0.05*yscale, radar_symbol) xscale = h_range_max_e+h_range_max_w yscale = h_range_max_s+h_range_max_n CALL XSTPJGRD(mapproj,trulat1,trulat2,trulon,ctrlat,ctrlon, & xscale,yscale,xorig,yorig) CALL xwindw(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) IF( ovrmap /= 0 ) THEN DO i=1,MIN(mxmapfile,nmapfile) lmapfile=len_trim(mapfile(i)) WRITE(6,'(1x,a,a)') 'Input was ',trim(mapfile(i)) INQUIRE(FILE=trim(mapfile(i)), EXIST = fexist ) IF( .NOT.fexist) THEN WRITE(6,'(a)') 'Map file '//mapfile(i)(1:lmapfile)// & ' not found. Corresponding map no plotted.' ELSE CALL xdrawmap(10,mapfile(i)(1:lmapfile),latgrid,longrid) END IF END DO END IF CALL xwdwof CALL xframe END DO ! ktilt DEALLOCATE( xrplt ) DEALLOCATE( yrplt ) DEALLOCATE( reftilt ) DEALLOCATE( irwxy ) DEALLOCATE( xrw ) DEALLOCATE( yrw ) END IF ! n_elevation_ref>0 ! !----------------------------------------------------------------------- ! To plot vertical cross-sections of reflectivity !----------------------------------------------------------------------- ! IF(n_azimuth_ref > 0) THEN ALLOCATE( rplt(maxrefgate,nref_tilts)) ALLOCATE( rpltm(maxrefgate,nref_tilts)) ALLOCATE( zplt(maxrefgate,nref_tilts)) ALLOCATE( datarz(maxrefgate,nref_tilts) ) ALLOCATE( iwrz(maxrefgate,nref_tilts) ) ALLOCATE( xrw(8*maxrefgate) ) ALLOCATE( yrw(8*maxrefgate) ) DO k=1,nref_tilts DO i=1,maxrefgate rad = 0.001 * rngrvol(i,k) rplt(i,k)= rad*cos( elvref(k)*deg2rad ) rpltm(i,k)=-rad*cos( elvref(k)*deg2rad ) zplt(i,k)= rad*sin( elvref(k)*deg2rad ) END DO END DO DO iazimuth= 1,n_azimuth_ref WRITE(6,'(a,f10.2)') & ' Plotting vertical cross-section at azimuth angle = ', & azimuth_ref(iazimuth) CALL xpspac(0.1,0.9,0.2,0.6) z_range_max = 12.0 CALL xmap(-h_range_max, h_range_max, 0.0, z_range_max) CALL xlabmask(1) IF( coltab_ref == -1 ) CALL xstctfn(coltabfn_ref) CALL xsetclrs(coltab_ref) ! Use internally defined color map number 6 CALL xcolor(1) ! Set color to black CALL xcfont(3) ! ! Set tentative color interval and limit on the min/max number of contours ! cl(1)=0.0 cl(2)=5.0 ! Set tentative contour interval CALL xnctrs(2,100) ! Set lower and upper limits ! on the number of contours allowed ! It's used by both xcolfil and xclimit CALL xnctrs(5,20) ! Set lower and upper limits mode=2 CALL xctrbadv(1) ! Turn on bad value checking ! in contouring routine CALL xbadval( missing_value ) ! Specify bad value flag ! ! plot a color filled contour field for positive values ! CALL xctrclr(ibgncol_ref,iendcol_ref) ! Use colors between 10 and 39 inclusive CALL xctrlim(5.0, 75.0) ! Only contours above zero are plotted CALL xwindw(-h_range_max, h_range_max, 0.0, z_range_max) DO k=1,nref_tilts jj=1 dazmin=720. DO j=1,maxazim delazm=azmrvol(j,k)-azimuth_ref(iazimuth) IF(delazm > 180.) delazm=delazm-360. IF(delazm < -180.) delazm=delazm+360. delazm=abs(delazm) IF(delazm < dazmin) THEN jj=j dazmin=delazm END IF END DO DO i=1,maxrefgate datarz(i,k) = max(refvol(i,jj,k),-30.) END DO END DO CALL xcolfil(datarz,rplt,zplt,iwrz,xrw,yrw, & maxrefgate,maxrefgate,nref_tilts,cl,ncl,mode) tem0 = azimuth_ref(iazimuth)+180.0 IF( tem0 > 360.) tem0=tem0-360. DO k=1,nref_tilts jj=1 dazmin=720. DO j=1,maxazim delazm=azmrvol(j,k)-tem0 IF(delazm > 180.) delazm=delazm-360. IF(delazm < -180.) delazm=delazm+360. delazm=abs(delazm) IF(delazm < dazmin) THEN jj=j dazmin=delazm END IF END DO DO i=1,maxrefgate datarz(i,k) = max(refvol(i,jj,k),-30.) END DO END DO CALL xcolfil(datarz,rpltm,zplt,iwrz,xrw,yrw, & maxrefgate,maxrefgate,nref_tilts,cl,ncl,mode) CALL xwdwof CALL xaxsca(-h_range_max, h_range_max, 0.1*h_range_max, & 0.0, z_range_max, 0.5 ) WRITE(string,'(a,f6.2)') & ' Reflectivity at Azimuth Angle=',azimuth_ref(iazimuth) CALL xchsiz(0.06*z_range_max ) ! Set character size CALL xcharc(0.0, z_range_max*1.10, trim(string) ) WRITE(stringtime,'(a,i2.2,a,i2.2,a,i4,a,i2.2,a,i2.2)') 'First Scan Time:',& imon,'/',iday,'/',iyear, & ' ',ihour,':',imin CALL xchsiz(0.05*z_range_max ) ! Set character size CALL xcharc(0.0, z_range_max*1.04, trim(stringtime) ) CALL xchsiz(0.04*z_range_max ) ! Set character size CALL xcpalet(1) ! Plot horizontal color palette CALL xframe END DO ! iazimuth DEALLOCATE( rplt ) DEALLOCATE( rpltm ) DEALLOCATE( zplt ) DEALLOCATE( iwrz ) DEALLOCATE( xrw ) DEALLOCATE( yrw ) DEALLOCATE( datarz ) END IF ! n_azimuth_ref ! ! Velocity data plotting section ! IF (n_elevation_vel > 0 ) THEN WRITE(6,'(a)') ' Allocating x,y,veltilt ...' ALLOCATE( xvplt(maxvelgate,(maxazim+1)) ) ALLOCATE( yvplt(maxvelgate,(maxazim+1)) ) ALLOCATE( veltilt(maxvelgate,(maxazim+1)) ) ALLOCATE( ivwxy(maxvelgate,(maxazim+1)) ) ALLOCATE( xvw(8*maxvelgate) ) ALLOCATE( yvw(8*maxvelgate) ) IF ( overlayrefcontour == 1 ) THEN ALLOCATE( xrplt(maxrefgate,(maxazim+1)) ) ALLOCATE( yrplt(maxrefgate,(maxazim+1)) ) ALLOCATE( reftilt(maxrefgate,(maxazim+1)) ) ALLOCATE( irwxy(maxrefgate,(maxazim+1)) ) ALLOCATE( xrw(8*maxrefgate) ) ALLOCATE( yrw(8*maxrefgate) ) END IF WRITE(6,'(a)') ' Radial Velocity tilt plotting...' DO ktilt = 1,n_elevation_vel elev_to_plot = elevations_vel(ktilt) kvel = 0 DO k = 1,maxelev IF( abs( elvref(k)-elev_to_plot ) & < elev_error_range ) THEN kvel = k EXIT END IF END DO IF(kvel == 0 ) CYCLE ! no elevation angle ! within error range was found velmin = 1000.0 velmax =-1000.0 velmin = 1000.0 velmax =-1000.0 eleva=elvref(kvel) last_azim=-999. jazim=0 DO j=1,maxazim+1 jidx=j IF(jidx > maxazim) jidx=jidx-maxazim data_azim=azmvvol(jidx,kvel) IF(data_azim >= 0. .AND. (abs(data_azim-last_azim) > 1.0E-04)) THEN jazim=jazim+1 last_azim=data_azim xconst=cos( (90.-data_azim)*deg2rad ) yconst=sin( (90.-data_azim)*deg2rad ) DO i=1,maxvelgate rad = 0.001 * rngvvol(i,kvel) xvplt(i,jazim)= rad * xconst yvplt(i,jazim)= rad * yconst veltilt(i,jazim)=max(velvol(i,jidx,kvel),-150.) IF(velvol(i,jidx,kvel) > missing_value ) THEN velmax=max(velvol(i,jidx,kvel),velmax) velmin=min(velvol(i,jidx,kvel),velmin) END IF END DO END IF END DO nazim=jazim WRITE(6,'(a,2f10.2)') ' Velmin, Velmax on tilt= ', velmin, velmax WRITE(6,'(a,i5)') 'N unique azimuths: ',nazim ! !----------------------------------------------------------------------- ! Define plotting space. ! CALL xpspac(0.5-factorx,0.5+factorx,0.5-factory,0.5+factory) CALL xmap(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) CALL xlabmask(1) IF( coltab_vel == -1 ) CALL xstctfn(coltabfn_vel) CALL xsetclrs(coltab_vel) ! Use internally defined color map number 6 CALL xcolor(1) ! Set color to black CALL xcfont(3) ! call xaxfmt('(I4)') ! ! Set tentative color interval and limit on the min/max number of contours ! cl(1)=0.0 cl(2)=5.0 ! Set tentative contour interval CALL xnctrs(2,100) ! Set lower and upper limits ! on the number of contours allowed ! It's used by both xcolfil and xclimit CALL xnctrs(5,20) ! Set lower and upper limits mode=2 call xctrbadv(1) ! Turn on bad value checking ! in contouring routine call xbadval( missing_value ) ! Specify bad value flag CALL xctrclr(ibgncol_vel,iendcol_vel) CALL xctrlim(-50.0, 50.0) CALL xwindw(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) WRITE(6,'(a,f10.2)') & 'Plotting colfil radial velocity at elevation angle=',eleva CALL xcolfil(veltilt,xvplt,yvplt,ivwxy,xvw,yvw, & maxvelgate,maxvelgate,nazim, cl, ncl,mode) IF ( overlayrefcontour == 1 ) THEN kref = 0 DO k = 1,maxelev IF( abs( elvref(k)-elev_to_plot ) & < elev_error_range ) THEN kref = k EXIT END IF END DO IF(kref /= 0 ) THEN refmin = 1000.0 refmax =-1000.0 velmin = 1000.0 refmax =-1000.0 eleva=elvref(kref) last_azim=-999. jazim=0 DO j=1,maxazim+1 jidx=j IF(jidx > maxazim) jidx=jidx-maxazim data_azim=azmrvol(jidx,kref) IF(data_azim >= 0. .AND. (abs(data_azim-last_azim) > 1.0E-04)) THEN jazim=jazim+1 last_azim=data_azim xconst=cos( (90.-data_azim)*deg2rad ) yconst=sin( (90.-data_azim)*deg2rad ) DO i=1,maxrefgate rad = 0.001 * rngrvol(i,kref) xrplt(i,jazim)= rad * xconst yrplt(i,jazim)= rad * yconst reftilt(i,jazim)=max(refvol(i,jidx,kref),-30.) IF(refvol(i,jidx,kref) > missing_value ) THEN refmax=max(refvol(i,jidx,kref),refmax) refmin=min(refvol(i,jidx,kref),refmin) END IF END DO END IF END DO nazim=jazim call XCLTYP(0) CALL xctrclr(14,15) ! Use colors between 10 and 39 inclusive mode=3 ncl=1 cl(1)=15.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=30.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=45.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) cl(1)=55.0 CALL XCONTA(reftilt,xrplt,yrplt,irwxy, & maxrefgate,maxrefgate,nazim, CL, NCL,MODE) END IF END IF CALL xwdwof CALL xaxsca(-h_range_max_w, h_range_max_e, 0.1*h_range_max, & -h_range_max_s, h_range_max_n, 0.1*h_range_max) WRITE(string,'(a,f5.2)') 'Radial Velocity at Elevation Angle=',eleva CALL xchsiz(0.025*2*h_range_max ) ! Set character size CALL xcharc((-h_range_max_w+h_range_max_e)*0.5, & h_range_max_n+2*h_range_max*0.04, trim(string) ) WRITE(stringtime,'(a,i2.2,a,i2.2,a,i4,a,i2.2,a,i2.2)') 'Scan Time:', & imon,'/',iday,'/',iyear, ' ',ihour,':',imin CALL xchsiz(0.025*2*h_range_max ) ! Set character size CALL xcharc((-h_range_max_w+h_range_max_e)*0.5, & h_range_max_n+2*h_range_max*0.01, trim(stringtime) ) CALL xchsiz(0.02* 2*h_range_max ) ! Set character size CALL xcpalet(1) ! Plot horizontal color palette xscale = h_range_max yscale = h_range_max xorig = -(h_range_max_e+h_range_max_w)/2.0 yorig = -(h_range_max_s+h_range_max_n)/2.0 xradar = 0.0 yradar = 0.0 CALL xthick(3) CALL xpenup( xradar-xscale*0.02, yradar ) CALL xpendn( xradar+xscale*0.02, yradar ) CALL xpenup( xradar, yradar-yscale*0.02 ) CALL xpendn( xradar, yradar+yscale*0.02 ) call xthick(1) CALL xchsiz(0.02* 2*h_range_max ) ! Set character size CALL xcharc(xradar, yradar-0.05*yscale, radar_symbol) xscale = h_range_max_e+h_range_max_w yscale = h_range_max_s+h_range_max_n CALL XSTPJGRD(mapproj,trulat1,trulat2,trulon,ctrlat,ctrlon, & xscale,yscale,xorig,yorig) CALL xwindw(-h_range_max_w, h_range_max_e, -h_range_max_s, h_range_max_n) IF( ovrmap /= 0 ) THEN DO i=1,MIN(mxmapfile,nmapfile) lmapfile=len_trim(mapfile(i)) WRITE(6,'(1x,a,a)') 'Input was ',trim(mapfile(i)) INQUIRE(FILE=trim(mapfile(i)), EXIST = fexist ) IF( .NOT.fexist) THEN WRITE(6,'(a)') 'Map file '//mapfile(i)(1:lmapfile)// & ' not found. Corresponding map no plotted.' ELSE CALL xdrawmap(10,mapfile(i)(1:lmapfile),latgrid,longrid) END IF END DO END IF CALL xwdwof CALL xframe END DO ! ktilt DEALLOCATE( ivwxy ) DEALLOCATE( xvw ) DEALLOCATE( yvw ) DEALLOCATE( xvplt ) DEALLOCATE( yvplt ) DEALLOCATE( veltilt ) IF(overlayrefcontour == 1) THEN DEALLOCATE( irwxy ) DEALLOCATE( xrw ) DEALLOCATE( yrw ) DEALLOCATE( xrplt ) DEALLOCATE( yrplt ) DEALLOCATE( reftilt ) END IF END IF ! n_elevation_vel>0 ! !----------------------------------------------------------------------- ! To plot vertical cross-sections of radial velocity !----------------------------------------------------------------------- ! IF(n_azimuth_vel > 0) THEN ALLOCATE( rplt(maxvelgate,nvel_tilts)) ALLOCATE( rpltm(maxvelgate,nvel_tilts)) ALLOCATE( zplt(maxvelgate,nvel_tilts)) ALLOCATE( datarz(maxvelgate,nvel_tilts) ) ALLOCATE( iwrz(maxvelgate,nvel_tilts) ) ALLOCATE( xvw(8*maxvelgate) ) ALLOCATE( yvw(8*maxvelgate) ) DO k=1,nvel_tilts DO i=1,maxvelgate rad = 0.001 * rngvvol(i,k) rplt(i,k)= rad*cos( elvref(k)*deg2rad ) rpltm(i,k)=-rad*cos( elvref(k)*deg2rad ) zplt(i,k)= rad*sin( elvref(k)*deg2rad ) END DO END DO DO iazimuth= 1,n_azimuth_vel WRITE(6,'(a,f10.2)') & ' Plotting vertical cross-section at azimuth angle = ', & azimuth_vel(iazimuth) CALL xpspac(0.1,0.9,0.2,0.6) z_range_max = 12.0 CALL xmap(-h_range_max, h_range_max, 0.0, z_range_max) CALL xlabmask(1) IF( coltab_vel == -1 ) CALL xstctfn(coltabfn_vel) CALL xsetclrs(coltab_vel) ! Use internally defined color map number 6 CALL xcolor(1) ! Set color to black CALL xcfont(3) ! ! Set tentative color interval and limit on the min/max number of contours ! cl(1)=0.0 cl(2)=5.0 ! Set tentative contour interval CALL xnctrs(2,100) ! Set lower and upper limits ! on the number of contours allowed ! It's used by both xcolfil and xclimit CALL xnctrs(5,20) ! Set lower and upper limits mode=2 CALL xctrbadv(1) ! Turn on bad value checking ! in contouring routine CALL xbadval( missing_value ) ! Specify bad value flag ! ! plot a color filled contour field for positive values ! CALL xctrclr(ibgncol_vel,iendcol_vel) ! Use colors between 10 and 39 inclusive CALL xctrlim(50.0, 50.0) ! Only contours above zero are plotted CALL xwindw(-h_range_max, h_range_max, 0.0, z_range_max) DO k=1,nvel_tilts jj=1 dazmin=720. DO j=1,maxazim delazm=azmrvol(j,k)-azimuth_vel(iazimuth) IF(delazm > 180.) delazm=delazm-360. IF(delazm < -180.) delazm=delazm+360. delazm=abs(delazm) IF(delazm < dazmin) THEN jj=j dazmin=delazm END IF END DO DO i=1,maxvelgate datarz(i,k) = max(velvol(i,jj,k),-30.) END DO END DO CALL xcolfil(datarz,rplt,zplt,iwrz,xvw,yvw, & maxvelgate,maxvelgate,nvel_tilts,cl,ncl,mode) tem0 = azimuth_vel(iazimuth)+180.0 IF( tem0 > 360.) tem0=tem0-360. DO k=1,nvel_tilts jj=1 dazmin=720. DO j=1,maxazim delazm=azmrvol(j,k)-tem0 IF(delazm > 180.) delazm=delazm-360. IF(delazm < -180.) delazm=delazm+360. delazm=abs(delazm) IF(delazm < dazmin) THEN jj=j dazmin=delazm END IF END DO DO i=1,maxvelgate datarz(i,k) = max(velvol(i,jj,k),-30.) END DO END DO CALL xcolfil(datarz,rpltm,zplt,iwrz,xvw,yvw, & maxvelgate,maxvelgate,nvel_tilts,cl,ncl,mode) CALL xwdwof CALL xaxsca(-h_range_max, h_range_max, 0.1*h_range_max, & 0.0, z_range_max, 0.5 ) WRITE(string,'(a,f6.2)') & ' Radial Velocity at Azimuth Angle=',azimuth_vel(iazimuth) CALL xchsiz(0.06*z_range_max ) ! Set character size CALL xcharc(0.0, z_range_max*1.10, trim(string) ) WRITE(stringtime,'(a,i2.2,a,i2.2,a,i4,a,i2.2,a,i2.2)') 'First Scan Time:',& imon,'/',iday,'/',iyear, & ' ',ihour,':',imin CALL xchsiz(0.05*z_range_max ) ! Set character size CALL xcharc(0.0, z_range_max*1.04, trim(stringtime) ) CALL xchsiz(0.04*z_range_max ) ! Set character size CALL xcpalet(1) ! Plot horizontal color palette CALL xframe END DO ! iazimuth DEALLOCATE( rplt ) DEALLOCATE( rpltm ) DEALLOCATE( zplt ) DEALLOCATE( iwrz ) DEALLOCATE( xvw ) DEALLOCATE( yvw ) DEALLOCATE( datarz ) END IF ! n_azimuth_vel ! !----------------------------------------------------------------------- ! End of processing for this file !----------------------------------------------------------------------- ! DEALLOCATE( rngrvol) DEALLOCATE( azmrvol) DEALLOCATE( elvrvol) DEALLOCATE( refvol) DEALLOCATE( rngvvol) DEALLOCATE( azmvvol) DEALLOCATE( elvvvol) DEALLOCATE( velvol) IF(wrtuaref /= 0) DEALLOCATE(uarefvol) IF(wrtuavel /= 0) DEALLOCATE(uavelvol) END DO ! ifiles IF( zxplot_initialized == 1 ) CALL xgrend ! End graphics STOP 999 100 CONTINUE Write(6,'(a)') 'Error reading namelist input file.' Write(6,'(a)') 'Program aborted.' STOP 101 END PROGRAM PLT_RADAR_TILT ! !################################################################## !################################################################## !###### ###### !###### Advanced Regional Prediction System (ARPS) ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! SUBROUTINE OUTPUT_REGULAR(regular_file, ntilts_max,ntilts, & radar_symbol,radar_lat,radar_lon,radar_elevation, & ngate,nray,fstgat,gatsp,int_razim,eleva, & volume ) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! Out put 3d regular polar grid data for further retrival. ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Ming Hu ! 03/19/02. ! ! MODIFICATION HISTORY: ! ! INPUT: ! ! regular_file file name that data will be saved into ! ntilts_max maximum number of tilts ! ntilts actual number of tilts in this observation data ! radar_symbol radar Name ! radar_lat radar location: latitude ! radar_lon radar location: longitude ! radar_elevation radar location: elavation ! ! ngate number of gates of observation data per radial ! nray number of radials of data( here = 361 ) ! fstgat distance from radar of first gate of data(meters) ! gatsp gate spacing of data (meters) ! ! volume(ngate,nray,ntilts) Array containing one full volume scan of ! observation(reflectivity or radial velocity ! which has been intepolated to 1-degree ! regular spaced polar angles ! int_razim(nray) azimuth angle for each radial (degrees) ! eleva(ntilts_max) elevation angle for each tilt ! ! ! ! OUTPUT: ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE !----------------------------------------------------------------------- ! ! INPUT RADAR PARAMETER AND OBSERVATION ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: regular_file INTEGER :: ntilts_max INTEGER :: ntilts CHARACTER(LEN=4) :: radar_symbol REAL :: radar_lat REAL :: radar_lon REAL :: radar_elevation INTEGER :: ngate INTEGER :: nray REAL :: fstgat REAL :: gatsp REAL :: volume(ngate,nray,ntilts) ! Array containing one ! full volume scan of ! reflectivity intepolated to ! 1-degree spaced polar angles REAL :: int_razim(nray) REAL :: eleva(ntilts_max) ! !----------------------------------------------------------------------- ! ! WRITE(*,*) regular_file ! WRITE(*,*) ntilts_max,ntilts,ngate,nray,ntilts ! WRITE(*,*) radar_symbol,radar_lat,radar_lon,radar_elevation ! WRITE(*,*) ngate,nray,fstgat,gatsp ! WRITE(*,*) int_razim ! WRITE(*,*) eleva ! WRITE(*,*) volume ! OPEN(27,file=trim(regular_file),status='unknown',form='unformatted') WRITE(27) ntilts_max,ntilts WRITE(27) radar_symbol,radar_lat,radar_lon,radar_elevation WRITE(27) ngate,nray,fstgat,gatsp WRITE(27) int_razim WRITE(27) eleva WRITE(27) volume CLOSE(27) WRITE(*,*) ' The file:',trim(regular_file), & 'has been written out successfully' END SUBROUTINE OUTPUT_REGULAR