!######################################################################## !######################################################################## !######### ######### !######### PROGRAM ncrad2arps ######### !######### ######### !######### Developed by ######### !######### Center for Analysis and Prediction of Storms ######### !######### University of Oklahoma ######### !######### ######### !######################################################################## !######################################################################## PROGRAM ncrad2arps,50 !------------------------------------------------------------------------ ! ! PURPOSE: ! ! Reads radar data from NetCDF data files and remaps data on ARPS ! grid to be used for analysis or display. ! CASA radar data, for example, are stored in netCDF files. ! !------------------------------------------------------------------------ ! ! AUTHOR: ! ! Keith Brewster (May 2006) ! ! MODIFICATIONS: ! ! 1 Feb 2007 Keith Brewster ! Added processing for gzipped files. ! !------------------------------------------------------------------------ IMPLICIT NONE !------------------------------------------------------------------------ ! ! Include files ! !------------------------------------------------------------------------ INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'netcdf.inc' !------------------------------------------------------------------------ ! ! Dimensions and parameters ! !------------------------------------------------------------------------ INTEGER, PARAMETER :: maxfile = 90 INTEGER, PARAMETER :: maxvar = 10 ! INTEGER, PARAMETER :: stdout = 6 INTEGER, PARAMETER :: itime1970 = 315619200 REAL, PARAMETER :: refchek=-35. REAL, PARAMETER :: refmiss=-999. REAL, PARAMETER :: velchek=-190. REAL, PARAMETER :: velmiss=-999. REAL, PARAMETER :: winszazim=10. REAL, PARAMETER :: winszrad=1000. !------------------------------------------------------------------------ ! ! Variable Declarations. ! !------------------------------------------------------------------------ INTEGER :: nx, ny, nz, nzsoil, itime, nstyps REAL :: radelv, radlat, radlon REAL :: radarx, radary !------------------------------------------------------------------------ ! ! CASA Radar Data Variables ! !------------------------------------------------------------------------ ! REAL, ALLOCATABLE :: refl(:,:) REAL, ALLOCATABLE :: attref(:,:) REAL, ALLOCATABLE :: rhohv(:,:) REAL, ALLOCATABLE :: radv(:,:) ! !------------------------------------------------------------------------ ! ! Variables for subroutine remapvol. ! !------------------------------------------------------------------------ ! INTEGER, PARAMETER :: nt=1 INTEGER, PARAMETER :: exbcbufsz=1 INTEGER, PARAMETER :: bkgopt=1 INTEGER, PARAMETER :: shropt=1 INTEGER, PARAMETER :: rfropt = 1 ! 4/3rds earth refraction option INTEGER, PARAMETER :: maxsort = 10000 INTEGER, PARAMETER :: nzsnd = 201 REAL, PARAMETER :: dzsnd=100. INTEGER :: ktsnd(nzsnd) REAL :: zsnd(nzsnd) ! hgt levels for refractivity sounding REAL :: usnd(nzsnd) REAL :: vsnd(nzsnd) REAL :: rfrsnd(nzsnd) ! refractivity sounding REAL, PARAMETER :: rfrconst = (-1./(4.*6371.E03)) ! -1/(4e) REAL, PARAMETER :: envavgr=80000. REAL, ALLOCATABLE :: arfdata(:,:) ! attenuated refl (dBZ) data REAL, ALLOCATABLE :: refdata(:,:) ! refl (dBZ) data REAL, ALLOCATABLE :: veldata(:,:) ! velocity data (m/s) REAL, ALLOCATABLE :: spwdata(:,:) REAL, ALLOCATABLE :: bkgvel(:,:) REAL, ALLOCATABLE :: unfvdata(:,:) ! unfolded velocity data (m/s) REAL, ALLOCATABLE :: rhodata(:,:) ! rho-HV data REAL, ALLOCATABLE :: snstvty(:) ! radar sensitivity REAL, ALLOCATABLE :: rtem(:,:) ! temporary array for radar data processing INTEGER, ALLOCATABLE :: time(:) ! time offset from time1st REAL, ALLOCATABLE :: azim(:) ! azimuth angle for each radial(degree) REAL, ALLOCATABLE :: beamw(:) ! beamwidth for each radial(degree) REAL, ALLOCATABLE :: gtspc(:) ! gate spacing for each radial (mneter) REAL, ALLOCATABLE :: vnyq(:) ! Nyquist velocity (meters/sec) REAL, ALLOCATABLE :: elev(:) ! elevation angle for each radial (degree) INTEGER, ALLOCATABLE :: istrgate(:) INTEGER, ALLOCATABLE :: bgate(:) INTEGER, ALLOCATABLE :: egate(:) INTEGER :: time1st ! time of first radial in volume INTEGER :: irfirstg ! range to first gate (meters) INTEGER :: igatesp ! gate spacing (meters) INTEGER :: irefgatsp ! reflectivity gate spacing (meters) INTEGER :: ivelgatsp ! velocity gate spacing (meters) REAL :: elv ! elevation angle for this tilt REAL :: rfirstg ! range to first gate (m) REAL :: gatesp ! gate spacing (meters) REAL :: vnyquist ! Nyquist velocity for this tilt INTEGER, ALLOCATABLE :: kntrgat(:,:) INTEGER, ALLOCATABLE :: kntrazm(:) INTEGER :: kntrelv INTEGER, ALLOCATABLE :: kntvgat(:,:) INTEGER, ALLOCATABLE :: kntvazm(:) INTEGER :: kntvelv INTEGER, ALLOCATABLE :: timevol(:,:) REAL, ALLOCATABLE :: nyqvvol(:) REAL, ALLOCATABLE :: rngrvol(:,:) REAL, ALLOCATABLE :: azmrvol(:,:) REAL, ALLOCATABLE :: elvrvol(:,:) REAL, ALLOCATABLE :: refvol(:,:,:) REAL, ALLOCATABLE :: rngvvol(:,:) REAL, ALLOCATABLE :: azmvvol(:,:) REAL, ALLOCATABLE :: elvvvol(:,:) REAL, ALLOCATABLE :: velvol(:,:,:) REAL, ALLOCATABLE :: elvmean(:) REAL, ALLOCATABLE :: varsort(:) REAL, ALLOCATABLE :: rxvol(:,:,:) REAL, ALLOCATABLE :: ryvol(:,:,:) REAL, ALLOCATABLE :: rzvol(:,:,:) REAL, ALLOCATABLE :: gridvel(:,:,:) REAL, ALLOCATABLE :: gridref(:,:,:) REAL, ALLOCATABLE :: gridnyq(:,:,:) REAL, ALLOCATABLE :: gridtim(:,:,:) REAL, ALLOCATABLE :: x(:) REAL, ALLOCATABLE :: y(:) REAL, ALLOCATABLE :: z(:) REAL, ALLOCATABLE :: zp(:,:,:) REAL, ALLOCATABLE :: xs(:) REAL, ALLOCATABLE :: ys(:) REAL, ALLOCATABLE :: zps(:,:,:) REAL, ALLOCATABLE :: hterain(:,:) REAL, ALLOCATABLE :: mapfct(:,:,:) ! Map factors at scalar, u and v points REAL, ALLOCATABLE :: j1 (:,:,:) ! Coordinate transformation Jacobian defined ! as - d( zp )/d( x ). REAL, ALLOCATABLE :: j2 (:,:,:) ! Coordinate transformation Jacobian defined ! as - d( zp )/d( y ). REAL, ALLOCATABLE :: j3 (:,:,:) ! Coordinate transformation Jacobian defined ! as d( zp )/d( z ). REAL, ALLOCATABLE :: tem1d1(:) REAL, ALLOCATABLE :: tem1d2(:) REAL, ALLOCATABLE :: tem2d(:,:) REAL, ALLOCATABLE :: tem2dyz(:,:) REAL, ALLOCATABLE :: tem2dxz(:,:) REAL, ALLOCATABLE :: tem2dns(:,:,:) REAL, ALLOCATABLE :: tem3d1(:,:,:) REAL, ALLOCATABLE :: tem3d2(:,:,:) REAL, ALLOCATABLE :: tem3d3(:,:,:) REAL, ALLOCATABLE :: tem3d4(:,:,:) INTEGER, ALLOCATABLE :: tem2dint(:,:) INTEGER, ALLOCATABLE :: soiltyp(:,:,:) REAL, ALLOCATABLE :: stypfrct(:,:,:) REAL, ALLOCATABLE :: tem3dsoil(:,:,:) REAL, ALLOCATABLE :: tem4dsoilns(:,:,:,:) REAL, ALLOCATABLE :: zpsoil(:,:,:) REAL, ALLOCATABLE :: j3soil(:,:,:) REAL, ALLOCATABLE :: j3soilinv(:,:,:) REAL, ALLOCATABLE :: u(:,:,:) REAL, ALLOCATABLE :: v(:,:,:) REAL, ALLOCATABLE :: qv(:,:,:) REAL, ALLOCATABLE :: ubar(:,:,:) REAL, ALLOCATABLE :: vbar(:,:,:) REAL, ALLOCATABLE :: ptbar(:,:,:) REAL, ALLOCATABLE :: uprt(:,:,:) REAL, ALLOCATABLE :: vprt(:,:,:) REAL, ALLOCATABLE :: ptprt(:,:,:) REAL, ALLOCATABLE :: pprt(:,:,:) REAL, ALLOCATABLE :: pbar(:,:,:) REAL, ALLOCATABLE :: qvbar(:,:,:) REAL, ALLOCATABLE :: rhostr(:,:,:) REAL, ALLOCATABLE :: trigs1(:) REAL, ALLOCATABLE :: trigs2(:) INTEGER, ALLOCATABLE :: ifax(:) REAL, ALLOCATABLE :: wsave1(:) REAL, ALLOCATABLE :: wsave2(:) REAL, ALLOCATABLE :: vwork1(:,:) REAL, ALLOCATABLE :: vwork2(:,:) REAL, ALLOCATABLE :: qcumsrc(:,:,:,:) REAL, ALLOCATABLE :: prcrate(:,:,:) REAL, ALLOCATABLE :: exbcbuf(:) REAL, ALLOCATABLE :: grdvel2(:,:,:) REAL, ALLOCATABLE :: grdref2(:,:,:) REAL, ALLOCATABLE :: wgtvel(:,:,:) REAL, ALLOCATABLE :: wgtref(:,:,:) REAL, ALLOCATABLE :: wpotvel(:,:,:) REAL, ALLOCATABLE :: wpotref(:,:,:) REAL, ALLOCATABLE :: gridazm(:,:) REAL, ALLOCATABLE :: gridrng(:,:) REAL(KIND=8), ALLOCATABLE :: tem_double(:) !------------------------------------------------------------------------ ! ! Data dimensions, determined from the data files ! !------------------------------------------------------------------------ INTEGER :: maxgate,maxazim,maxelev !------------------------------------------------------------------------ ! ! Misc. variables ! !------------------------------------------------------------------------ CHARACTER(LEN=256) :: infilname CHARACTER(LEN=80) :: indir CHARACTER(LEN=80) :: dir_name CHARACTER(LEN=256) :: fname CHARACTER(LEN=256) :: full_fname CHARACTER(LEN=6) :: velid = 'radv3d' CHARACTER(LEN=6) :: refid = 'refl3d' CHARACTER(LEN=20) :: refname = 'Reflectivity' CHARACTER(LEN=20) :: velname = 'RawVelocity' CHARACTER(LEN=6) :: refunits = 'dBZ' CHARACTER(LEN=6) :: velunits = 'm/s' CHARACTER(LEN=NF_MAX_NAME), ALLOCATABLE :: ncvarname(:) CHARACTER(LEN=266) :: syscall INTEGER :: fvartype(maxfile) LOGICAL :: gzipped(maxfile) REAL :: felv(maxfile) INTEGER :: fitime(maxfile) INTEGER :: n, i, j, k, istatus INTEGER :: nazim,ngate,nvar INTEGER :: ivcp,itimcdf,initime,mintimcdf INTEGER :: ivar,varfill,dmpfmt,hdf4cmpr INTEGER :: len_dir,len_fname INTEGER :: arbvaropt,velopt LOGICAL :: velproc LOGICAL :: vel2d,vel3d LOGICAL :: ref2d,ref3d LOGICAL :: qc_on LOGICAL :: rho_match INTEGER :: iyear,iyr,imonth,iday,ihr,imin,isec INTEGER :: iiyr,iimonth,iiday,iihr,iimin,iisec INTEGER :: velknt,refknt,rhoknt,kntt2a,kntt2b INTEGER :: isource,icount,iradtype INTEGER :: nyqset,timeset,vardump INTEGER :: xscale REAL :: rdx,dtime,rmin,rmax REAL :: rmisval,rngfval,frtime REAL :: refelvmin,refelvmax REAL :: refrngmin,refrngmax INTEGER :: remapopt REAL :: radius,vconst,bmwidth REAL :: xrad,yrad INTEGER :: maxinfile,infile,lenfn INTEGER :: iarg,jarg,iargc,narg,ifile,jfile,kfile CHARACTER(LEN=4) tsthead CHARACTER(LEN=256) charg CHARACTER(LEN=256) listfile CHARACTER (LEN=6) :: varid CHARACTER (LEN=20) :: varname CHARACTER (LEN=20) :: varunits !------------------------------------------------------------------------ ! ! Some tunable parameters - convert to Namelist in future ! !------------------------------------------------------------------------ INTEGER, PARAMETER :: maxncradfile=90 INTEGER, PARAMETER :: iordref = 2 INTEGER, PARAMETER :: iordvel = 1 REAL, PARAMETER :: refmedl = 15. REAL, PARAMETER :: velmedl = 15. REAL, PARAMETER :: refdazl = 360. REAL, PARAMETER :: veldazl = 20. REAL, PARAMETER :: dazim = 1.0 REAL, PARAMETER :: rngmin = 500. REAL, PARAMETER :: rngmax = 230.E03 REAL, PARAMETER :: rhohvlim = 0.50 CHARACTER(LEN=4) :: radname CHARACTER(LEN=32) :: radarname INTEGER :: nncradfn = 0 CHARACTER(LEN=256) :: ncradfn(maxncradfile) !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ !------------------------------------------------------------------------ ! ! Initializations ! !------------------------------------------------------------------------ isource = 3 ! 1 - WSR-88D raw 2 - WSR-88D NIDS refelvmin = 91.0 refelvmax = -91.0 refrngmin = rngmin refrngmax = rngmax dmpfmt=1 hdf4cmpr=0 kntt2a=0 kntt2b=0 refknt=0 rhoknt=0 velknt=0 velproc = .TRUE. qc_on = .TRUE. ref2d = .FALSE. ref3d = .FALSE. vel2d = .FALSE. vel3d = .FALSE. !------------------------------------------------------------------------ ! ! Initialize ncrad_data NAMELIST variables ! !------------------------------------------------------------------------ radname = 'DMMY' nncradfn = 0 dir_name = './' len_dir = LEN(TRIM(dir_name)) arbvaropt = 0 DO ifile=1,maxncradfile ncradfn(ifile)='dummy' END DO narg=iargc() WRITE(stdout,'(a,i4)') ' Number of command-line arguments: ',narg IF(narg < 2 ) THEN WRITE(stdout,'(a,a,a,a)') & ' Usage: ncrad2arps RADAR_ID radarfile_listfile [-novel] [-hdf 2]', & ' [-binary] [-ref2d] [-ref3d] [-reffile] [-vel2d] [-vel3d] [-velfile]',& ' [-noqc]', & ' < arps.input' STOP END IF !------------------------------------------------------------------------ ! ! Read grid info from the arps.input ! !------------------------------------------------------------------------ CALL initpara(nx, ny, nz, nzsoil, nstyps)! Reads standard ARPS namelists !------------------------------------------------------------------------ ! ! Process command line ! For backward comptability allow input of radar file names via ! files specified on the command line ! !------------------------------------------------------------------------ IF(narg > 1) THEN CALL getarg(1,charg) radname=charg(1:4) WRITE(stdout,'(a,a)') ' radname = ', radname kfile=0 nncradfn=0 iarg=2 DO jarg=2,narg CALL getarg(iarg,charg) IF(charg(1:6) == '-novel') THEN velproc=.FALSE. ELSE IF(charg(1:4) == '-hdf') THEN dmpfmt=2 hdf4cmpr=0 IF(iarg < narg) THEN iarg=iarg+1 CALL getarg(iarg,charg) READ(charg,'(i1)') hdf4cmpr hdf4cmpr=min(max(hdf4cmpr,0),7) END IF WRITE(stdout,'(a,i2)') & ' Output in hdf format with compression level: ',hdf4cmpr ELSE IF(charg(1:7) == '-binary') THEN dmpfmt=1 hdf4cmpr=0 WRITE(stdout,'(a)') ' Output in binary format' ELSE IF(charg(1:6) == '-ref2d') THEN ref2d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 2d reflectivity file of lowest tilt' ELSE IF(charg(1:6) == '-ref3d') THEN ref3d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 3d reflectivity file for plotting' ELSE IF(charg(1:8) == '-reffile') THEN ref3d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 3d reflectivity file for plotting' ELSE IF(charg(1:6) == '-vel2d') THEN vel2d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 2d velocity file of lowest tilt' ELSE IF(charg(1:6) == '-vel3d') THEN vel3d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 3d velocity file for plotting' ELSE IF(charg(1:8) == '-velfile') THEN vel3d=.TRUE. WRITE(stdout,'(a)') & ' Will produce 3d velocity file for plotting' ELSE IF(charg(1:5) == '-noqc') THEN qc_on=.FALSE. WRITE(stdout,'(a)') & ' Will skip qc steps in processing' ELSE listfile=charg WRITE(stdout,'(a,a)') ' Reading file lists ',charg OPEN(31,file=listfile,status='old',form='formatted') DO ifile=(nncradfn+1),maxncradfile READ(31,'(a)',END=101) ncradfn(ifile) kfile=kfile+1 WRITE(stdout,'(a,i4,a,a)') & ' kfile:',kfile,' ncradfn: ',TRIM(ncradfn(ifile)) END DO 101 CONTINUE nncradfn=kfile END IF iarg=iarg+1 IF(iarg > narg) EXIT END DO WRITE(stdout,'(i4,a)') nncradfn,' file names read' END IF istatus = 0 icount = 0 rfirstg=1000. gatesp = 1000. irfirstg = NINT(rfirstg) igatesp = NINT(gatesp) !------------------------------------------------------------------------ ! ! Loop through the NetCDF radar files. ! !------------------------------------------------------------------------ IF (nncradfn > 0 .AND. nncradfn <= maxncradfile) THEN maxinfile = nncradfn ELSE IF (nncradfn > maxncradfile) THEN WRITE(stdout,'(a,i4,a,i4)') & ' WARNING: nncradfn=', nncradfn,' > maxncradfile=',maxncradfile WRITE(stdout,'(a,i4,a)') & ' Will only read in ', maxncradfile,' NIDS files.' maxinfile = maxncradfile ELSE maxinfile = 0 END IF IF (maxinfile > 0) THEN ! !------------------------------------------------------------------------ ! ! First check each file for variable stored and dimensions of data ! !------------------------------------------------------------------------ ! maxgate=1 maxazim=1 maxelev=1 mintimcdf=-1 DO infile = 1,maxinfile fvartype(infile)=0 gzipped(infile)=.FALSE. END DO DO infile = 1,maxinfile nazim = 0 ngate = 0 nvar = 0 infilname = ncradfn(infile) lenfn=len(infilname) CALL strlnth(infilname,lenfn) IF(infilname((lenfn-2):lenfn) == ".gz") THEN gzipped(infile)=.TRUE. WRITE(syscall,'(a,1x,a)') 'gunzip',TRIM(infilname) CALL system(syscall) infilname((lenfn-2):lenfn)=' ' ncradfn(infile)=infilname WRITE(6,'(a,a)') 'Uncompressed ',TRIM(infilname) END IF WRITE(stdout,'(a,a)') 'Checking dims in file:',TRIM(infilname) CALL extdirn(infilname,indir,fname) print *, ' infile: ',TRIM(infilname) print *, ' indir: ',TRIM(indir) print *, ' fname: ',TRIM(fname) CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, & istatus) WRITE(6,'(a,i6,a,i6)') ' istatus: ',istatus,' nvar in file: ',nvar IF(istatus == NF_NOERR .AND. nvar > 0) THEN ALLOCATE(ncvarname(nvar)) IF(iradtype == 21) THEN kntt2a=kntt2a+1 ALLOCATE(tem_double(nazim)) CALL get_ncrad2ainfo(fname,indir,nazim,nvar,tem_double, & radarname,radlat,radlon,radelv,itimcdf,frtime,elv, & ncvarname,istatus) DEALLOCATE(tem_double) fvartype(infile)=100 refknt=refknt+1 velknt=velknt+1 rhoknt=rhoknt+1 ELSE IF(iradtype == 22) THEN kntt2b=kntt2b+1 CALL get_ncrad2binfo(fname,indir,nazim,nvar, & radarname,radlat,radlon,radelv,itimcdf,frtime, & ivcp,elv,ncvarname,istatus) DO k=1,nvar print *, ' variable(',k,') = ',TRIM(ncvarname(k)) IF(ncvarname(k)(1:7) == 'Reflect') THEN fvartype(infile)=1 refknt=refknt+1 EXIT ELSE IF(ncvarname(k)(1:11) == 'RawVelocity' .OR. & ncvarname(k)(1:15) == 'AliasedVelocity' ) THEN fvartype(infile)=2 velknt=velknt+1 EXIT ELSE IF(ncvarname(k)(1:5) == 'RhoHV') THEN fvartype(infile)=3 rhoknt=rhoknt+1 EXIT END IF END DO ELSE WRITE(stdout,'(a,a)') 'Unknown radar file type:',iradtype fvartype(infile)=0 STOP END IF IF(radarname(1:5) == 'cyril') THEN IF(radelv < 100.) radelv=432.9 IF(elv < 0.1) elv=1.0 ELSE IF(radarname(1:5) == 'chick') THEN IF(radelv < 100.) radelv=354.0 IF(elv < 0.1) elv=1.0 ELSE IF(radarname(1:5) == 'lawto') THEN IF(radelv < 100.) radelv=377.4 IF(elv < 0.1) elv=1.0 ELSE IF(radarname(1:5) == 'rushs') THEN IF(radelv < 100.) radelv=414.8 IF(elv < 0.1) elv=1.0 END IF felv(infile)=elv fitime(infile)=itimcdf maxazim=max(nazim,maxazim) maxgate=max(ngate,maxgate) IF(mintimcdf < 0) THEN mintimcdf=itimcdf ELSE mintimcdf=min(mintimcdf,itimcdf) END IF DEALLOCATE(ncvarname) ELSE ! error return from get dims fvartype(infile)=0 END IF END DO maxelev=max(maxelev,refknt) maxelev=max(maxelev,velknt) IF(velknt == 0) velproc=.false. WRITE(6,'(/a,i6)') ' Number of reflectivity files:',refknt WRITE(6,'(a,i6)') ' Number of Rho-HV files:',rhoknt WRITE(6,'(a,i6)') ' Number of velocity files:',velknt WRITE(6,'(/a,i6)') ' Max azimuth dimension of NetCDF data:',maxazim WRITE(6,'(a,i6)') ' Max gate dimension of NetCDF data:',maxgate WRITE(6,'(a,i6/)') ' Max elevations among data:',maxelev IF(kntt2a > 0 .AND. kntt2b > 0) THEN WRITE(6,'(a)') ' Mixed Tier 2a and Tier 2b files in list, stopping' WRITE(6,'(a,i4,a,i4)')' N Tier 2a=',kntt2a,' N Tier 2b=',kntt2b STOP END IF ALLOCATE(x(nx)) ALLOCATE(y(ny)) ALLOCATE(z(nz)) ALLOCATE(zp(nx,ny,nz)) ALLOCATE(hterain(nx,ny)) ALLOCATE(mapfct(nx,ny,8)) ALLOCATE(j1(nx,ny,nz)) ALLOCATE(j2(nx,ny,nz)) ALLOCATE(j3(nx,ny,nz)) ALLOCATE(tem1d1(nz)) ALLOCATE(tem1d2(nz)) ALLOCATE(tem3d1(nx,ny,nz)) IF(qc_on .AND. velknt > 0) THEN !------------------------------------------------------------------------ ! ! Fill refractivity sounding with constant lapse rate ! This is only used when rfropt > 1, but for completeness and ! future upgrade (when an actual sounding made from gridded data ! would be used) this is filled-in. ! !------------------------------------------------------------------------ DO k=1,nzsnd zsnd(k)=(k-1)*dzsnd rfrsnd(k)=325.+rfrconst*zsnd(k) END DO !------------------------------------------------------------------------ ! ! ARPS grid array allocations and variable initializations. ! !------------------------------------------------------------------------ ALLOCATE(u(nx,ny,nz)) ALLOCATE(v(nx,ny,nz)) ALLOCATE(qv(nx,ny,nz)) ALLOCATE(uprt(nx,ny,nz)) ALLOCATE(vprt(nx,ny,nz)) ALLOCATE(ptprt(nx,ny,nz)) ALLOCATE(pprt(nx,ny,nz)) ALLOCATE(ubar(nx,ny,nz)) ALLOCATE(vbar(nx,ny,nz)) ALLOCATE(ptbar(nx,ny,nz)) ALLOCATE(pbar(nx,ny,nz)) ALLOCATE(qvbar(nx,ny,nz)) ALLOCATE(rhostr(nx,ny,nz)) ALLOCATE(tem2d(nx,ny)) ALLOCATE(tem2dyz(ny,nz)) ALLOCATE(tem2dxz(nx,nz)) ALLOCATE(tem2dns(nx,ny,0:nstyps)) ALLOCATE(trigs1(3*(nx-1)/2+1)) ALLOCATE(trigs2(3*(ny-1)/2+1)) ALLOCATE(ifax(13)) ALLOCATE(wsave1(3*(ny-1)+15)) ALLOCATE(wsave2(3*(nx-1)+15)) ALLOCATE(vwork1(nx+1,ny+1)) ALLOCATE(vwork2(ny,nx+1)) ALLOCATE(qcumsrc(nx,ny,nz,5)) ALLOCATE(prcrate(nx,ny,4)) ALLOCATE(exbcbuf(exbcbufsz)) ALLOCATE(soiltyp(nx,ny,nstyps)) ALLOCATE(stypfrct(nx,ny,nstyps)) ALLOCATE(tem2dint(nx,ny)) ALLOCATE(tem3d2(nx,ny,nz)) ALLOCATE(tem3d3(nx,ny,nz)) ALLOCATE(tem3d4(nx,ny,nz)) ALLOCATE(tem3dsoil(nx,ny,nzsoil)) ALLOCATE(tem4dsoilns(nx,ny,nzsoil,0:nstyps)) CALL set_lbcopt(1) CALL initgrdvar(nx,ny,nz,nzsoil,nt,nstyps,exbcbufsz, & x,y,z,zp,tem3dsoil,hterain,mapfct, & j1,j2,j3,tem3dsoil,tem3d1,tem3d1,tem3d1,tem3d1,tem3dsoil, & u,v,tem3d1,tem3d1,ptprt,pprt, & qv,tem3d1,tem3d1,tem3d1,tem3d1,tem3d1,tem3d1, & tem2dyz,tem2dyz,tem2dxz,tem2dxz, & tem2dyz,tem2dyz,tem2dxz,tem2dxz, & trigs1,trigs2,ifax,ifax, & wsave1,wsave2,vwork1,vwork2, & ubar,vbar,ptbar,pbar,tem3d1,tem3d1, & rhostr,tem3d1,qvbar,tem3d1,tem3d1, & soiltyp,stypfrct,tem2dint,tem2d,tem2d,tem2d, & tem4dsoilns,tem4dsoilns,tem2dns,tem2d,tem2dns, & tem3d1,qcumsrc,tem3d1,tem2dint,tem2d, & tem2d,tem2d,tem2d, & tem2d,tem2d,prcrate,exbcbuf, & tem3d1,tem2d,tem2d,tem2d,tem2d, & tem2d,tem2d,tem2d,tem2d, & tem3dsoil,tem3dsoil,tem3dsoil,tem3dsoil,tem3dsoil, & tem3d2,tem3d3,tem3d4,tem3d4,tem3d4, & tem3d4,tem3d4,tem3d4,tem3d4) print *, ' Back from initgrdvar' DEALLOCATE(j1) DEALLOCATE(j2) DEALLOCATE(j3) DEALLOCATE(trigs1) DEALLOCATE(trigs2) DEALLOCATE(ifax) DEALLOCATE(wsave1) DEALLOCATE(wsave2) DEALLOCATE(vwork1) DEALLOCATE(vwork2) DEALLOCATE(qcumsrc) DEALLOCATE(prcrate) DEALLOCATE(soiltyp) DEALLOCATE(stypfrct) DEALLOCATE(exbcbuf) DEALLOCATE(tem2dint) DEALLOCATE(tem2dyz) DEALLOCATE(tem2dxz) DEALLOCATE(tem2dns) DEALLOCATE(tem4dsoilns) DEALLOCATE(tem3dsoil) ALLOCATE(xs(nx)) ALLOCATE(ys(ny)) ALLOCATE(zps(nx,ny,nz)) print *, ' Setting radar coordinate: ' CALL radcoord(nx,ny,nz,x,y,z,zp,xs,ys,zps, & radlat,radlon,radarx,radary) print *, ' Radar x: ',(0.001*radarx),' km' print *, ' Radar y: ',(0.001*radary),' km' print *, ' Environmental wind averaging radius: ',envavgr CALL extenvprf(nx,ny,nz,nzsnd, & x,y,zp,xs,ys,zps, & u,v,ptprt,pprt,qv,ptbar,pbar,tem3d1,tem3d2,tem3d3, & radarx,radary,radelv,envavgr, & zsnd,ktsnd,usnd,vsnd,rfrsnd) print *, ' Back from extenvprf' DEALLOCATE(ubar) DEALLOCATE(vbar) DEALLOCATE(pbar) DEALLOCATE(ptbar) DEALLOCATE(qvbar) DEALLOCATE(rhostr) DEALLOCATE(uprt) DEALLOCATE(vprt) DEALLOCATE(ptprt) DEALLOCATE(pprt) DEALLOCATE(qv) DEALLOCATE(tem3d2) DEALLOCATE(tem3d3) DEALLOCATE(tem3d4) ELSE ALLOCATE(xs(nx)) ALLOCATE(ys(ny)) ALLOCATE(zps(nx,ny,nz)) ALLOCATE(zpsoil(nx,ny,nzsoil)) ALLOCATE(j3soil(nx,ny,nzsoil)) ALLOCATE(j3soilinv(nx,ny,nzsoil)) print *, ' Setting ARPS grid coordinates ...' CALL inigrd(nx,ny,nz,nzsoil,x,y,z,zp,zpsoil, & hterain,mapfct,j1,j2,j3,j3soil,j3soilinv, & tem1d1,tem1d2,tem3d1) print *, ' Setting radar coordinate...' CALL radcoord(nx,ny,nz,x,y,z,zp,xs,ys,zps, & radlat,radlon,radarx,radary) DEALLOCATE(j1) DEALLOCATE(j2) DEALLOCATE(j3) DEALLOCATE(zpsoil) DEALLOCATE(j3soil) DEALLOCATE(j3soilinv) END IF DEALLOCATE(hterain) DEALLOCATE(mapfct) DEALLOCATE(tem3d1) ALLOCATE(gridvel(nx,ny,nz)) ALLOCATE(gridref(nx,ny,nz)) ALLOCATE(gridnyq(nx,ny,nz)) ALLOCATE(gridtim(nx,ny,nz)) ALLOCATE(grdvel2(nx,ny,nz)) ALLOCATE(grdref2(nx,ny,nz)) ALLOCATE(wgtvel(nx,ny,nz)) ALLOCATE(wgtref(nx,ny,nz)) ALLOCATE(wpotvel(nx,ny,nz)) ALLOCATE(wpotref(nx,ny,nz)) ALLOCATE(gridazm(nx,ny)) ALLOCATE(gridrng(nx,ny)) ALLOCATE(varsort(maxsort)) !------------------------------------------------------------------------ ! ! Radar data array allocations ! !------------------------------------------------------------------------ IF(refknt > 0) THEN ALLOCATE(refdata(maxgate,maxazim)) ALLOCATE(arfdata(maxgate,maxazim)) END IF IF(velknt > 0) THEN ALLOCATE(veldata(maxgate,maxazim)) IF(qc_on) THEN ALLOCATE(spwdata(maxgate,maxazim)) ALLOCATE(unfvdata(maxgate,maxazim)) ALLOCATE(bkgvel(maxgate,maxazim)) ALLOCATE(bgate(maxazim)) ALLOCATE(egate(maxazim)) spwdata=0. END IF END IF IF(rhoknt > 0) ALLOCATE(rhodata(maxgate,maxazim)) ALLOCATE(snstvty(maxgate)) ALLOCATE(rtem(maxgate,maxazim)) ALLOCATE(time(maxazim)) ALLOCATE(azim(maxazim)) ALLOCATE(beamw(maxazim)) ALLOCATE(gtspc(maxazim)) ALLOCATE(vnyq(maxazim)) ALLOCATE(elev(maxazim)) ALLOCATE(istrgate(maxazim)) ALLOCATE(kntrgat(maxazim,maxelev)) ALLOCATE(kntrazm(maxelev)) ALLOCATE(kntvgat(maxazim,maxelev)) ALLOCATE(kntvazm(maxelev)) ALLOCATE(nyqvvol(maxelev)) ALLOCATE(timevol(maxazim,maxelev)) ALLOCATE(rngrvol(maxgate,maxelev)) ALLOCATE(azmrvol(maxazim,maxelev)) ALLOCATE(elvrvol(maxazim,maxelev)) ALLOCATE(refvol(maxgate,maxazim,maxelev)) ALLOCATE(rngvvol(maxgate,maxelev)) ALLOCATE(azmvvol(maxazim,maxelev)) ALLOCATE(elvvvol(maxazim,maxelev)) ALLOCATE(velvol(maxgate,maxazim,maxelev)) ALLOCATE(elvmean(maxelev)) ALLOCATE(rxvol(maxgate,maxazim,maxelev)) ALLOCATE(ryvol(maxgate,maxazim,maxelev)) ALLOCATE(rzvol(maxgate,maxazim,maxelev)) !------------------------------------------------------------------------ ! ! Initialize 3-d radar data arrays ! !------------------------------------------------------------------------ itimcdf=mintimcdf+itime1970 time1st=itimcdf CALL abss2ctim(itimcdf, iyear, imonth, iday, ihr, imin, isec) iyr = MOD(iyear, 100) CALL rmpinit(nx,ny,nz,maxgate,maxgate,maxazim,maxelev, & kntrgat,kntrazm,kntrelv, & kntvgat,kntvazm,kntvelv, & nyqvvol,timevol, & rngrvol,azmrvol,elvrvol,refvol, & rngvvol,azmvvol,elvvvol,velvol, & gridvel,gridref,gridnyq,gridtim) iiyr = iyr iimonth = imonth iiday = iday iihr = ihr iimin = imin iisec = isec IF(kntt2a > 0) THEN DO infile = 1,maxinfile infilname = ncradfn(infile) WRITE(6,'(a,i4,a,i4)') ' infile= ',infile,' type=',fvartype(infile) CALL extdirn(infilname,indir,fname) CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, & istatus) ALLOCATE(tem_double(nazim)) ALLOCATE(ncvarname(nvar)) CALL get_ncrad2ainfo(fname,indir,nazim,nvar,tem_double, & radarname,radlat,radlon,radelv,itimcdf,frtime,elv, & ncvarname,istatus) DEALLOCATE(ncvarname) ALLOCATE(attref(ngate,nazim)) ALLOCATE(refl(ngate,nazim)) ALLOCATE(radv(ngate,nazim)) ALLOCATE(rhohv(ngate,nazim)) CALL rd2atiltcdf(nazim,ngate,fname,indir, & rmisval,rngfval,itimcdf,frtime,initime, & vnyquist,rfirstg,istrgate, & azim,beamw,gtspc,attref,refl,radv,rhohv, & tem_double) DEALLOCATE(tem_double) bmwidth=beamw(1) arfdata=-999.0 refdata=-999.0 DO j=1,nazim DO i=istrgate(j),ngate IF(abs(refl(i,j)) < 99. ) THEN arfdata(i,j)=attref(i,j) refdata(i,j)=refl(i,j) END IF END DO END DO veldata=-999.0 DO j=1,nazim DO i=istrgate(j),ngate IF(abs(radv(i,j)) < 199. ) veldata(i,j)=radv(i,j) END DO END DO rhodata=-999.0 DO j=1,nazim DO i=istrgate(j),ngate IF(abs(rhohv(i,j)) < 1.1 ) rhodata(i,j)=rhohv(i,j) END DO END DO DEALLOCATE(attref) DEALLOCATE(refl) DEALLOCATE(radv) DEALLOCATE(rhohv) !------------------------------------------------------------------------ ! ! Process reflectivity data. ! !------------------------------------------------------------------------ WRITE(stdout,'(a,i4)') & ' Processing base reflectivity data... ', icount WRITE(stdout,*)'Transferring ',nazim,' reflectivity radials.' dtime=float(itime-time1st) refelvmin=min(elv,refelvmin) refelvmax=max(elv,refelvmax) DO j = 1, maxazim time(j) = dtime elev(j) = elv END DO WRITE(stdout,'(a,i5,a,4f9.1)') & ' Ref j azim elev refmin refmax' DO j = 1, nazim rmin=999. rmax=-999. DO i=1,ngate IF(refdata(i,j) > refchek) THEN rmin=min(refdata(i,j),rmin) rmax=max(refdata(i,j),rmax) END IF END DO IF(mod(j,5) == 0) & WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax END DO gatesp = gtspc(1) igatesp = NINT(gatesp) irfirstg = NINT(rfirstg) IF( qc_on) THEN print *, ' calling SNR sensitivity screen' CALL snrchk(maxgate,maxazim,ngate,nazim,refchek,rfirstg,gatesp, & snstvty,refdata,veldata,arfdata) print *, ' calling RhoHV screen' CALL rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, & refdata,veldata,rhodata) print *, ' calling despekl ' CALL despekl(maxgate,maxazim,maxgate,maxazim,refchek, & refdata,rtem) END IF print *, ' calling volbuild reflectivity' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,ngate,nazim, & nyqset,timeset, & kntrgat,kntrazm,kntrelv, & igatesp,irfirstg,refchek, & vnyquist,time, & azim,elev,refdata, & nyqvvol,timevol, & rngrvol,azmrvol,elvrvol,refvol) !------------------------------------------------------------------------ ! ! Process velocity data. ! !------------------------------------------------------------------------ icount = icount + 1 WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') & ' DATE: ', iyear, '/', imonth, '/', iday WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') & ' TIME: ', ihr, ':', imin, ':', isec WRITE(stdout,'(a,i4)')' Processing base velocity data... ', icount WRITE(stdout,'(a,i6,a)')'Transferring',nazim,' velocity radials.' WRITE(stdout,'(a,f10.2)') ' Nyquist velocity: ',vnyquist dtime=float(itime-time1st) DO j = 1, maxazim time(j) = dtime elev(j) = elv END DO WRITE(stdout,'(a,i5,a,4f9.1)') & ' Ref j azim elev vr min vr max' DO j = 1, nazim rmin=999. rmax=-999. DO i=1,ngate IF(veldata(i,j) > velchek) THEN rmin=min(veldata(i,j),rmin) rmax=max(veldata(i,j),rmax) END IF END DO IF(mod(j,5) == 0) & WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax END DO gatesp = gtspc(1) igatesp = NINT(gatesp) irfirstg = NINT(rfirstg) !------------------------------------------------------------------------ ! ! Call radar volume builder. ! !------------------------------------------------------------------------ ivelgatsp=igatesp nyqset=1 timeset=1 IF( qc_on ) THEN print *, ' calling despekl ' CALL despekl(maxgate,maxazim,maxgate,maxazim,velchek, & veldata,rtem) print *,' calling unfnqc' print *,' vnyquist=',vnyquist CALL unfnqc(maxgate,maxazim,maxgate,nazim, & nzsnd,zsnd,usnd,vsnd,rfrsnd, & bkgopt,shropt,rfropt,igatesp,irfirstg, & radlat,radlon,radelv, & veldata,spwdata,elev,azim,vnyquist, & unfvdata,bkgvel,bgate,egate,rtem) print *, ' calling volbuild velocity ' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, & nyqset,timeset, & kntvgat,kntvazm,kntvelv, & igatesp,irfirstg,velchek, & vnyquist,time, & azim,elev,unfvdata, & nyqvvol,timevol, & rngvvol,azmvvol,elvvvol,velvol) ELSE print *, ' calling volbuild velocity ' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, & nyqset,timeset, & kntvgat,kntvazm,kntvelv, & igatesp,irfirstg,velchek, & vnyquist,time, & azim,elev,veldata, & nyqvvol,timevol, & rngvvol,azmvvol,elvvvol,velvol) END IF END DO ! infile ELSE IF(kntt2b > 0) THEN ! !------------------------------------------------------------------------ ! ! Process reflectivity files ! !------------------------------------------------------------------------ ! DO infile = 1,maxinfile WRITE(6,'(a,i4,a,i4)') ' infile= ',infile, & ' type=',fvartype(infile) IF(fvartype(infile) == 1) THEN ! type is reflectivity infilname = ncradfn(infile) WRITE(stdout,'(a,/a)') & 'Reading reflectivity file:',TRIM(infilname) CALL extdirn(infilname,indir,fname) CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, & istatus) ALLOCATE(ncvarname(nvar)) CALL get_ncrad2binfo(fname,indir,nazim,nvar, & radarname,radlat,radlon,radelv,itimcdf,frtime, & ivcp,elv,ncvarname,istatus) DEALLOCATE(ncvarname) ALLOCATE(refl(ngate,nazim)) CALL rdrftiltcdf(nazim,ngate,fname,indir, & rmisval,rngfval,itimcdf,frtime,initime,rfirstg, & azim,beamw,gtspc,refl) bmwidth=beamw(1) refdata=rmisval DO j=1,nazim DO i=21,ngate refdata(i,j)=refl(i,j) END DO END DO DEALLOCATE(refl) IF(radarname(1:5) == 'cyril') THEN IF(radelv < 0.1) radelv=445.3 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'chick') THEN IF(radelv < 0.1) radelv=355.4 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'lawto') THEN IF(radelv < 0.1) radelv=295.0 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'rushs') THEN IF(radelv < 0.1) radelv=365.8 IF(elv < 0.1) elv=1.5 END IF !------------------------------------------------------------------------ ! ! Now begin processing the data. ! !------------------------------------------------------------------------ icount = icount + 1 WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') & ' DATE: ', iyear, '/', imonth, '/', iday WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') & ' TIME: ', ihr, ':', imin, ':', isec !------------------------------------------------------------------------ ! ! Process reflectivity data. ! !------------------------------------------------------------------------ WRITE(stdout,'(a,i4)') & ' Processing base reflectivity data... ', icount WRITE(stdout,*)'Transferring ',nazim,' reflectivity radials.' dtime=float(itime-time1st) refelvmin=min(elv,refelvmin) refelvmax=max(elv,refelvmax) DO j = 1, maxazim time(j) = dtime elev(j) = elv END DO WRITE(stdout,'(a,i5,a,4f9.1)') & ' Ref j azim elev refmin refmax' DO j = 1, nazim, 5 rmin=999. rmax=-999. DO i=1,ngate IF(refdata(i,j) > refchek) THEN rmin=min(refdata(i,j),rmin) rmax=max(refdata(i,j),rmax) END IF END DO WRITE(stdout,'(2x,i5,4f9.1)') j,azim(j),elev(j),rmin,rmax END DO gatesp = gtspc(1) igatesp = NINT(gatesp) rfirstg = NINT(rfirstg) !------------------------------------------------------------------------ ! ! Call radar volume builder. ! !------------------------------------------------------------------------ irefgatsp=igatesp print *, ' gatesp,rfirstg,refchek=',gatesp,rfirstg,refchek ! !------------------------------------------------------------------------ ! ! Get matching rho-HV data ! !------------------------------------------------------------------------ ! rho_match=.FALSE. DO jfile=1,maxinfile IF(fvartype(jfile) == 3 .AND. fitime(jfile) == itimcdf .AND. & abs(felv(jfile)-elv) < 0.01 ) THEN rho_match=.TRUE. infilname = ncradfn(jfile) WRITE(stdout,'(a,a)') 'Reading Rho-HV file:',TRIM(infilname) CALL extdirn(infilname,indir,fname) ALLOCATE(rhohv(ngate,nazim)) CALL rdrhotiltcdf(nazim,ngate,fname,indir,radarname, & radlat,radlon,radelv,elv, & rmisval,rngfval,itimcdf,frtime,initime, & rfirstg,azim,rhohv) rhodata=rmisval print *, ' ngate =',ngate print *, ' max range =',(ngate*0.024) DO j=1,nazim DO i=21,ngate rhodata(i,j)=rhohv(i,j) END DO END DO DEALLOCATE(rhohv) EXIT END IF END DO ! jfile loop IF(radarname(1:5) == 'cyril') THEN IF(radelv < 0.1) radelv=445.3 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'chick') THEN IF(radelv < 0.1) radelv=355.4 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'lawto') THEN IF(radelv < 0.1) radelv=295.0 IF(elv < 0.1) elv=1.5 ELSE IF(radarname(1:5) == 'rushs') THEN IF(radelv < 0.1) radelv=365.8 IF(elv < 0.1) elv=1.5 END IF IF( qc_on) THEN IF(rho_match) THEN print *, ' calling RhoHV screen' CALL rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, & refdata,veldata,rhodata) END IF print *, ' calling despekl ' CALL despekl(maxgate,maxazim,maxgate,maxazim,refchek, & refdata,rtem) END IF nyqset=0 IF ( velproc ) THEN timeset=0 ELSE timeset=1 END IF vnyquist=0. print *, ' calling volbuild reflectivity' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,ngate,nazim, & nyqset,timeset, & kntrgat,kntrazm,kntrelv, & igatesp,irfirstg,refchek, & vnyquist,time, & azim,elev,refdata, & nyqvvol,timevol, & rngrvol,azmrvol,elvrvol,refvol) END IF ! type is reflectivity END DO ! file loop !------------------------------------------------------------------------ ! ! Process radial velocity files ! !------------------------------------------------------------------------S ! DO infile = 1,maxinfile IF(fvartype(infile) == 2) THEN ! type is velocity infilname = ncradfn(infile) WRITE(stdout,'(a,a)') 'Reading velocity file:',TRIM(infilname) CALL extdirn(infilname,indir,fname) CALL get_ncraddims(fname,indir,iradtype,ngate,nazim,nvar, & istatus) ALLOCATE(ncvarname(nvar)) CALL get_ncrad2binfo(fname,indir,nazim,nvar, & radarname,radlat,radlon,radelv,itimcdf,frtime, & ivcp,elv,ncvarname,istatus) DEALLOCATE(ncvarname) ALLOCATE(radv(ngate,nazim)) CALL rdvrtiltcdf(nazim,ngate,fname,indir, & rmisval,rngfval,itimcdf,frtime,initime, & vnyquist,rfirstg, & azim,beamw,gtspc,vnyq,radv) bmwidth=beamw(1) veldata=-999. DO j=1,nazim DO i=1,ngate veldata(i,j)=radv(i,j) END DO END DO DEALLOCATE(radv) !------------------------------------------------------------------------ ! ! Correction for early experimental CASA data ! !------------------------------------------------------------------------ IF(radarname(1:5) == 'cyril' .OR. & radarname(1:5) == 'chick' .OR. & radarname(1:5) == 'lawto' .OR. & radarname(1:5) == 'rushs' ) THEN vconst=15./acos(-1.) DO j=1,nazim DO i=1,ngate IF(veldata(i,j) > -4.0) THEN veldata(i,j)=veldata(i,j)*vconst END IF END DO END DO vnyquist=15.0 vnyq=15.0 END IF !------------------------------------------------------------------------ ! ! Process velocity data. ! !------------------------------------------------------------------------ icount = icount + 1 WRITE(stdout,'(a,i4)') 'VCP number for this file:', ivcp WRITE(stdout,'(1x,a,i4.4,a1,i2.2,a1,i2.2)', ADVANCE='NO') & ' DATE: ', iyear, '/', imonth, '/', iday WRITE(stdout,'(1X,a,i2.2,a1,i2.2,a1,i2.2)') & ' TIME: ', ihr, ':', imin, ':', isec WRITE(stdout,'(a,i4)')' Processing base velocity data... ', icount WRITE(stdout,'(a,i6,a)')'Transferring',nazim,' velocity radials.' WRITE(stdout,'(a,f10.2)') ' Nyquist velocity: ',vnyquist dtime=float(itime-time1st) DO j = 1, maxazim time(j) = dtime elev(j) = elv END DO DO j = 1, nazim, 20 WRITE(stdout,'(a,i5,a,f8.1)') & ' Vel j =',j,' azim =',azim(j) END DO gatesp = gtspc(1) igatesp = NINT(gtspc(1)) irfirstg = NINT(rfirstg) !------------------------------------------------------------------------ ! ! Call radar volume builder. ! !------------------------------------------------------------------------ ivelgatsp=igatesp nyqset=1 timeset=1 IF( qc_on ) THEN print *, ' calling despekl ' CALL despekl(maxgate,maxazim,maxgate,maxazim,velchek, & veldata,rtem) print *, ' calling unfnqc' print *, ' Nyquist velocity: ',vnyquist CALL unfnqc(maxgate,maxazim,maxgate,nazim, & nzsnd,zsnd,usnd,vsnd,rfrsnd, & bkgopt,shropt,rfropt,igatesp,irfirstg, & radlat,radlon,radelv, & veldata,spwdata,elev,azim,vnyquist, & unfvdata,bkgvel,bgate,egate,rtem) print *, ' calling volbuild velocity ' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, & nyqset,timeset, & kntvgat,kntvazm,kntvelv, & igatesp,irfirstg,velchek, & vnyquist,time, & azim,elev,unfvdata, & nyqvvol,timevol, & rngvvol,azmvvol,elvvvol,velvol) ELSE print *, ' calling volbuild velocity ' print *, ' elv=',elev(1),' nazim =',nazim CALL volbuild(maxgate,maxazim,maxelev,maxgate,nazim, & nyqset,timeset, & kntvgat,kntvazm,kntvelv, & igatesp,irfirstg,velchek, & vnyquist,time, & azim,elev,veldata, & nyqvvol,timevol, & rngvvol,azmvvol,elvvvol,velvol) END IF END IF ! velocity END DO ! DO ifile = 1,maxinfile END IF ! knt2b > 0 ! ! Restore any gzipped files ! DO infile=1,maxinfile IF(gzipped(infile)) THEN WRITE(syscall,'(a,1x,a)') 'gzip',TRIM(ncradfn(infile)) CALL system(syscall) END IF END DO END IF ! IF (maxinfile > 0) !------------------------------------------------------------------------ ! ! Create filename for output of remapped reflectivity and velocity. ! !------------------------------------------------------------------------ IF (icount > 0) THEN !------------------------------------------------------------------------ ! ! Call remapping routines ! !------------------------------------------------------------------------ IF( qc_on ) THEN print *, ' Calling apdetect ' CALL apdetect(maxgate,maxgate,maxazim,maxelev, & kntrgat,kntrazm,kntrelv, & kntvgat,kntvazm,kntvelv, & refchek,velchek, & irefgatsp,ivelgatsp, & winszrad,winszazim,ivcp, & rngrvol,azmrvol,elvrvol, & rngvvol,azmvvol,elvvvol, & refvol,velvol,rtem, & istatus) END IF nyqset=0 IF( velproc ) THEN timeset=0 ELSE timeset=1 END IF IF(ref3d) THEN vardump=1 ELSE vardump=0 END IF varfill=0 ivar=1 print *, ' Calling remapvol for reflectivity ' CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz,nzsnd,maxsort, & vardump,varfill,ivar,nyqset,timeset,rfropt, & refid,refname,refunits, & refchek,refmiss,bmwidth,refmedl,refdazl,iordref, & kntrgat,kntrazm,kntrelv, & radlat,radlon,radarx,radary,radelv,dazim, & rngmin,rngmax,time1st, & rngrvol,azmrvol,elvrvol, & refvol,timevol,nyqvvol,rxvol,ryvol,rzvol, & xs,ys,zps,zsnd,rfrsnd,varsort,elvmean, & gridref,gridtim,gridnyq,istatus) IF( ref2d ) THEN vardump=1 ivar=1 varid='refl2d' varname='Low-level reflect' varunits='dBZ' print *, ' Calling remap2d for reflectivity' CALL remap2d(maxgate,maxazim,maxelev,nx,ny,nzsnd,maxsort, & vardump,ivar,rfropt,varid,varname,varunits,dmpfmt,hdf4cmpr, & refchek,refmiss,refmedl,refdazl,iordref, & kntrgat,kntrazm,kntrelv, & radlat,radlon,radarx,radary,radelv,dazim, & rngmin,rngmax,rngrvol,azmrvol,elvrvol, & refvol,rxvol,ryvol,xs,ys,zsnd,rfrsnd,varsort, & tem2d,istatus) END IF IF( velproc ) THEN nyqset=1 timeset=1 IF(vel3d) THEN vardump=1 ELSE vardump=0 END IF vardump=0 varfill=0 ivar=2 WRITE(stdout,'(a)') ' Calling remapvol for velocity ' CALL remapvol(maxgate,maxazim,maxelev,nx,ny,nz,nzsnd,maxsort, & vardump,varfill,ivar,nyqset,timeset,rfropt, & velid,velname,velunits, & velchek,velmiss,bmwidth,velmedl,veldazl,iordvel, & kntvgat,kntvazm,kntvelv, & radlat,radlon,radarx,radary,radelv,dazim, & rngmin,rngmax,time1st, & rngvvol,azmvvol,elvvvol, & velvol,timevol,nyqvvol,rxvol,ryvol,rzvol, & xs,ys,zps,zsnd,rfrsnd,varsort,elvmean, & gridvel,gridtim,gridnyq,istatus) END IF IF( vel2d ) THEN vardump=1 ivar=1 varid='radv2d' varname='Low-level Velocity' varunits='m/s' print *, ' Calling remap2d for velocity ' CALL remap2d(maxgate,maxazim,maxelev,nx,ny,nzsnd,maxsort, & vardump,ivar,rfropt,varid,varname,varunits,dmpfmt,hdf4cmpr, & velchek,velmiss,velmedl,veldazl,iordvel, & kntvgat,kntvazm,kntvelv, & radlat,radlon,radarx,radary,radelv,dazim, & rngmin,rngmax,rngvvol,azmvvol,elvvvol, & velvol,rxvol,ryvol,xs,ys,zsnd,rfrsnd,varsort, & tem2d,istatus) END IF !------------------------------------------------------------------------ ! ! Create filename and write file. ! !------------------------------------------------------------------------ full_fname=' ' CALL mkradfnm(dmpfmt,dir_name,len_dir,radname,iiyr,iimonth,iiday, & iihr, iimin, iisec, full_fname, len_fname) WRITE(stdout,'(a,a)') ' Filename for this volume: ',TRIM(full_fname) ! CALL wtradcol(nx, ny, nz, dmpfmt, 1, hdf4cmpr, full_fname, & ! radname, radlat, radlon, radelv, & ! iyr,imonth,iday,iihr,iimin,iisec,ivcp,isource, & ! refelvmin,refelvmax,refrngmin,refrngmax, & ! xs, ys, zps, gridvel, gridref, gridnyq, gridtim, & ! tem1d1,tem1d2,tem1d3,tem1d4,tem1d5,tem1d6) CALL wtradcol(nx,ny,nz,dmpfmt,1,hdf4cmpr,full_fname, & radname,radlat,radlon,radelv, & iyr,imonth,iday,iihr,iimin,iisec,ivcp,isource, & refelvmin,refelvmax,refrngmin,refrngmax, & xs,ys,zps,gridvel,gridref,gridnyq,gridtim,tem2d); END IF !------------------------------------------------------------------------ ! ! The End. ! !------------------------------------------------------------------------ WRITE(stdout,'(/a)') ' Normal termination of ncrad2arps.' STOP END PROGRAM ncrad2arps SUBROUTINE rhohvchk(maxgate,maxazim,ngate,nazim,refchek,rhohvlim, & 2 refl,radv,rhohv) IMPLICIT NONE INTEGER, INTENT(IN) :: maxgate INTEGER, INTENT(IN) :: maxazim INTEGER, INTENT(IN) :: ngate INTEGER, INTENT(IN) :: nazim REAL, INTENT(IN) :: refchek REAL, INTENT(IN) :: rhohvlim REAL, INTENT(INOUT) :: refl(maxgate,maxazim) REAL, INTENT(INOUT) :: radv(maxgate,maxazim) REAL, INTENT(IN) :: rhohv(maxgate,maxazim) ! REAL, PARAMETER :: refmiss=-901.0 REAL, PARAMETER :: velmiss=-901.0 INTEGER, PARAMETER :: nhist=20 INTEGER :: histknt(nhist) INTEGER :: igate,jazim,ihist INTEGER :: kntvalid,kntflag REAL :: histcst,histinv,histbin,pctflg,pcthist histknt=0 histcst=1.0/float(nhist) histinv=float(nhist) kntvalid=0 kntflag=0 DO jazim=1,nazim DO igate=1,ngate IF(refl(igate,jazim) > refchek ) THEN kntvalid=kntvalid+1 IF( rhohv(igate,jazim) < rhohvlim) THEN kntflag=kntflag+1 refl(igate,jazim) = refmiss radv(igate,jazim) = velmiss END IF ihist=1+INT(rhohv(igate,jazim)*histinv) ihist=min(max(ihist,1),nhist) histknt(ihist)=histknt(ihist)+1 END IF END DO END DO ! pctflg=0. IF(kntvalid > 0) THEN pctflg=100.*float(kntflag)/float(kntvalid) END IF WRITE(6,'(a,i9,a,i9,a,f5.1,a)') & ' Rho-HV screened ',kntflag,' of ',kntvalid,' = ',pctflg, & '% reflectivity gates' ! WRITE(6,'(/a)') ' Histogram of Rho-HV (count, percent):' DO ihist=1,nhist histbin=histcst*(ihist-1) pcthist=100.*float(histknt(ihist))/float(kntvalid) WRITE(6,'(f7.3,i9,f10.1)') histbin,histknt(ihist),pcthist END DO RETURN END SUBROUTINE rhohvchk SUBROUTINE snrchk(maxgate,maxazim,ngate,nazim,refchek,rfirstg,gatesp, & 1,1 snstvty,refl,radv,attref) IMPLICIT NONE INTEGER, INTENT(IN) :: maxgate INTEGER, INTENT(IN) :: maxazim INTEGER, INTENT(IN) :: ngate INTEGER, INTENT(IN) :: nazim REAL, INTENT(IN) :: refchek REAL, INTENT(IN) :: rfirstg REAL, INTENT(IN) :: gatesp REAL, INTENT(OUT) :: snstvty(maxgate) REAL, INTENT(INOUT) :: refl(maxgate,maxazim) REAL, INTENT(INOUT) :: radv(maxgate,maxazim) REAL, INTENT(IN) :: attref(maxgate,maxazim) ! REAL, PARAMETER :: refmiss=-911.0 REAL, PARAMETER :: velmiss=-911.0 INTEGER :: igate,jazim INTEGER :: kntvalid,kntflag REAL :: snr,pctflg CALL xbsnstvty(maxgate,rfirstg,gatesp,snstvty) DO igate=1,maxgate,10 print *, ' igate: ',igate,' range=',(0.001*(rfirstg+(igate-1)*gatesp)), & ' sens= ',snstvty(igate) END DO kntvalid=0 kntflag=0 DO jazim=1,nazim DO igate=1,ngate IF(refl(igate,jazim) > refchek ) THEN kntvalid=kntvalid+1 IF( attref(igate,jazim) < snstvty(igate)) THEN kntflag=kntflag+1 refl(igate,jazim) = refmiss radv(igate,jazim) = velmiss END IF END IF END DO END DO ! pctflg=0. IF(kntvalid > 0) THEN pctflg=100.*float(kntflag)/float(kntvalid) END IF WRITE(6,'(a,i9,a,i9,a,f5.1,a)') & ' SNR screened ',kntflag,' of ',kntvalid,' = ',pctflg, & '% reflectivity gates' RETURN END SUBROUTINE snrchk SUBROUTINE xbsnstvty(ngate,rfirstg,gatesp,sens) 1 ! !----------------------------------------------------------------------- ! Calculate X-band radar sensitivity (dBZ) as a function of range. ! ! Based on information provided by Francesc Joyvount ! ! Keith Brewster and Erin Fay, CAPS !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: ngate REAL, INTENT(IN) :: rfirstg REAL, INTENT(IN) :: gatesp REAL, INTENT (OUT):: sens(ngate) ! Sensitivity values for all dist from radar. ! ! Parameters for CASA radar ! ! REAL, PARAMETER :: lmbda = 3.0E-02 ! Wavelength (m) REAL, PARAMETER :: freqkhz=11505000.0 ! kHz REAL, PARAMETER :: ptxdbm=66.68 ! Peak power (dBm) ! REAL, PARAMETER :: pt = 12500 ! Peak power (Watts) REAL, PARAMETER :: G = 38.0 ! Antenna gain (dB) REAL, PARAMETER :: F = 5.5 ! Noise figure (dB) ! REAL, PARAMETER :: tau = 0.67E-06 ! Radar pulse length (s) REAL, PARAMETER :: tau = 4.06E-07 ! Radar pulse length (s) REAL, PARAMETER :: theta = 1.8 ! Antenna half-power beamwidth (deg) REAL, PARAMETER :: B = 2.5 ! Bandwidth (MHz) REAL, PARAMETER :: lm = 1.0 ! Receiver mis-match loss(dB) REAL, PARAMETER :: Kw2 = 0.91 ! Refractive Index of water squarred REAL, PARAMETER :: T0 = 300. ! Rx temperature (K) REAL, PARAMETER :: Ta = 200. ! Antenna temperature (K) ! ! Physical constants ! REAL, PARAMETER :: K = 1.38E-23 ! Boltsmann's Constant (J/K) REAL, PARAMETER :: c = 2.99792E08 ! Speed of light (m/s) REAL, PARAMETER :: rconst =1.0E18 ! m6/m3 to mm6/m3 ! ! Misc internal variables ! INTEGER :: igate REAL :: ln2,pi,pi3,four3,bwrad,rnoise,BHz,Ni,sconst,rlm,rG REAL :: pt,lmbda REAL :: range ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! pt=0.001*(10.0**(0.1*ptxdbm)) print *, ' pt= ',pt lmbda=c/(freqkhz*1000.) print *, ' lmbda (m)= ',lmbda ln2=log(2.0) pi=acos(-1.) pi3=pi*pi*pi four3=4.*4.*4. bwrad=theta*pi/180. rnoise=10.**(0.1*F) rlm=10.**(0.1*lm) rG=10.**(0.1*G) BHz=1.0E06*b Ni=K*(Ta+(rnoise-1.0)*T0)*BHz sconst=rconst *(Ni/(pi3*Kw2)) * ((8.*ln2)/(bwrad*bwrad)) * & (2./(c*tau)) * ((lmbda*lmbda*four3*rlm)/(Pt*rG*rG)) DO igate = 1, ngate range=rfirstg+((igate-1)*gatesp) sens(igate) = 10.*LOG10(range*range*sconst) END DO END SUBROUTINE xbsnstvty SUBROUTINE extdirn(filename,dirname,fname) 5 ! ! Extract directory name from complete filename and output ! directory name and filename in two character string variables ! IMPLICIT NONE CHARACTER (LEN=*) :: filename CHARACTER (LEN=*) :: dirname CHARACTER (LEN=*) :: fname INTEGER :: lenf INTEGER :: i lenf=LEN(filename) dirname='./' fname=filename DO i=(lenf-1),1,-1 IF(filename(i:i) == '/') THEN dirname=filename(1:(i-1)) fname=filename((i+1):lenf) EXIT END IF END DO RETURN END SUBROUTINE extdirn