PROGRAM fakerad,18
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM FAKERAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Generate fake radar data from an ARPS history file.
!
! It shares with the model the include file 'globcst.inc'
! for storage parameters.
!
! It reads in a history file produced by ARPS 4.0 in a user
! specified format.
!
! Parameters grdin,basin,mstin,icein,trbin are read in from the
! data file itself, therefore are determined internally.
! Arrays that are not read in retain their initial zero values.
! These parameters are passed among subroutines through
! a common block defined in 'indtflg.inc'.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS, December, 1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (km)
! zp z coordinate of grid points in computational space (m)
! zpsoil z coordinate of soil levels (m)
! uprt x component of perturbation velocity (m/s)
! vprt y component of perturbation velocity (m/s)
! wprt Vertical component of perturbation velocity in Cartesian
! coordinates (m/s).
!
! ptprt Perturbation potential temperature (K)
! pprt Perturbation pressure (Pascal)
!
! qvprt Perturbation water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state air density (kg/m**3)
! qvbar Base state water vapor mixing ratio (kg/kg)
!
! CALCULATED DATA ARRAYS:
!
! u x component of velocity (m/s)
! v y component of velocity (m/s)
! w z component of velocity (m/s)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'indtflg.inc'
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! Soil level depth.
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 :: j3soil(:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z )
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: km (:,:,:) ! The turbulent mixing coefficient for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: rhobar(:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific humidity
! (kg/kg)
INTEGER, ALLOCATABLE :: soiltyp(:,:) ! Soil type
INTEGER, ALLOCATABLE :: vegtyp(:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! soil moisture (g/kg)
REAL, ALLOCATABLE :: wetcanp(:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing (:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc (:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: radswnet (:,:) ! Net solar radiation, SWin - SWout
REAL, ALLOCATABLE :: radlwin (:,:) ! Incoming longwave radiation
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: rhoprt(:,:,:) ! Perturbation air density (kg/m**3)
REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation water vapor specific
! humidity (kg/kg)
!
!-----------------------------------------------------------------------
!
! Other data variables
!
!-----------------------------------------------------------------------
!
REAL :: time
!
!-----------------------------------------------------------------------
!
! Temporary working arrays:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
!
!-----------------------------------------------------------------------
!
! ARPS dimensions:
!
! nx, ny, nz: Dimensions of computatioal grid. When run on MPP
! with PVM or MPI, they represent of the size of the
! subdomains. See below.
!
! Given nx, ny and nz, the physical domain size will be
! xl=(nx-3)*dx by yl=(ny-3)*dy by zh=(nz-3)*dz. The
! variables nx, ny, nz, dx, dy and dz are read in from
! the input file by subroutine INITPARA.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nz ! Number of grid points in the z-direction
INTEGER :: nzsoil ! Number of soil levels
!
!-----------------------------------------------------------------------
!
! User request stuff
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=4) :: radid
CHARACTER (LEN=80) :: grdbasfn,filename
CHARACTER (LEN=100) :: rfname
REAL :: latrad,lonrad,elvrad
REAL :: thref,elvmin,elvmax,dstmin,dstmax
!
!-----------------------------------------------------------------------
!
! Scalar grid location variables
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: lat(:,:),lon(:,:)
!-----------------------------------------------------------------------
!
! Radar variables
!
!-----------------------------------------------------------------------
!
INCLUDE 'remap.inc'
REAL, ALLOCATABLE :: gridvel(:,:,:)
REAL, ALLOCATABLE :: gridref(:,:,:)
REAL, ALLOCATABLE :: gridnyq(:,:,:)
REAL, ALLOCATABLE :: gridtim(:,:,:)
REAL, ALLOCATABLE :: xsc(:),ysc(:)
REAL, ALLOCATABLE :: zpsc(:,:,:)
!
!-----------------------------------------------------------------------
!
! Radar angles
!
!-----------------------------------------------------------------------
!
INTEGER :: maxtilt
PARAMETER (maxtilt=14)
REAL :: hlfwid
PARAMETER (hlfwid=0.45)
INTEGER :: ntilt
REAL :: ang(maxtilt)
REAL :: ang11(maxtilt),ang12(maxtilt)
REAL :: ang31(maxtilt),ang32(maxtilt)
DATA ang11 / 0.5,1.5,2.4,3.4,4.3,5.3,6.2,7.5,8.7, &
10.0,12.0,14.0,16.7,19.5/
DATA ang12 / 0.5,1.5,2.4,3.4,4.3,6.0,9.9,14.6,19.5, &
19.5,19.5,19.5,19.5,19.5/
DATA ang31 / 0.5,1.5,2.5,3.5,4.5,4.5,4.5,4.5,4.5, &
4.5, 4.5, 4.5, 4.5, 4.5/
DATA ang32 / 0.5,1.5,2.5,3.5,4.5,4.5,4.5,4.5,4.5, &
4.5, 4.5, 4.5, 4.5, 4.5/
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=60) :: vcptxt
INTEGER :: i,j,k,itilt,ireturn,mode
INTEGER :: ngchan,nchanl,lenfil,lengbf,nchin
INTEGER :: iradfmt,vcpnum,myr
INTEGER :: knt,knttry,kntref,kntdst,kntang
INTEGER :: abstsec,iyr,imon,idy,ihr,imin,isec
REAL :: latnot(2)
REAL :: xctr,yctr,x0sc,y0sc
REAL :: azim,dist,delz,eleva,range,dhdr,dsdr
REAL :: ddrot,uazmrad,vazmrad
REAL :: perc,prctry,prcref,prcdst,prcang
INTEGER :: istatus
!-----------------------------------------------------------------------
!
! NAMELIST DECLARATION
!
!-----------------------------------------------------------------------
! NAMELIST /grid_dims/ nx,ny,nz
NAMELIST /file_name/ hdmpfmt,grdbasfn,filename
NAMELIST /radar_info/ radid, latrad, lonrad, elvrad, vcpnum, dstmin, &
dstmax, thref
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(/11(/5x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to FAKERAD ###', &
'### This program reads in the history dump data ###', &
'### sets generated by ARPS, and generates a fake ###', &
'### radar data file. ###', &
'### ###', &
'###############################################################', &
'###############################################################'
!-----------------------------------------------------------------------
!
! Read in the dimensions
!
!-----------------------------------------------------------------------
! READ(5,grid_dims,END=100)
! WRITE(6,'(/a)')' Namelist grid_dims sucessfully read in.'
!
!-----------------------------------------------------------------------
!
! Get the name of the input data set.
!
!-----------------------------------------------------------------------
!
READ(5,file_name,END=100)
WRITE(6,'(/a)')' Namelist file_name sucessfully read in.'
lengbf=LEN_trim(grdbasfn)
lenfil=LEN_trim(filename)
WRITE(6,'(/a,a)')' The data set name is ', filename(1:lenfil)
! CALL get_dims_from_data(hinfmt,grdbasfn, nx,ny,nz,nstyps, ireturn)
CALL get_dims_from_data
(hinfmt,grdbasfn, nx,ny,nz,nzsoil,nstyps, ireturn)
!-----------------------------------------------------------------------
!
! Initialize variables.
!
!-----------------------------------------------------------------------
!
istatus=0
ALLOCATE( x=0
ALLOCATE( y=0
ALLOCATE( z=0
ALLOCATE(zp(nx,ny,nz),STAT=istatus)
zp=0
ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus)
zpsoil=0
ALLOCATE(j1(nx,ny,nz),STAT=istatus)
j1=0
ALLOCATE(j2(nx,ny,nz),STAT=istatus)
j2=0
ALLOCATE(j3(nx,ny,nz),STAT=istatus)
j3=0
ALLOCATE(j3(nx,ny,nzsoil),STAT=istatus)
j3soil=0
ALLOCATE( u=0
ALLOCATE( v=0
ALLOCATE( w=0
ALLOCATE(qc(nx,ny,nz),STAT=istatus)
qc=0
ALLOCATE(qr(nx,ny,nz),STAT=istatus)
qr=0
ALLOCATE(qi(nx,ny,nz),STAT=istatus)
qi=0
ALLOCATE(qs(nx,ny,nz),STAT=istatus)
qs=0
ALLOCATE(qh(nx,ny,nz),STAT=istatus)
qh=0
ALLOCATE(km(nx,ny,nz),STAT=istatus)
km=0
ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
ubar=0
ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
vbar=0
ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
wbar=0
ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
ptbar=0
ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
rhobar=0
ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
pbar=0
ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
qvbar=0
ALLOCATE(soiltyp(nx,ny),STAT=istatus)
soiltyp=0
ALLOCATE(vegtyp(nx,ny),STAT=istatus)
vegtyp=0
ALLOCATE(lai(nx,ny),STAT=istatus)
lai=0
ALLOCATE(roufns(nx,ny),STAT=istatus)
roufns=0
ALLOCATE(veg(nx,ny),STAT=istatus)
veg=0
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tsoil")
tsoil = 0
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qsoil")
qsoil = 0
ALLOCATE(wetcanp(nx,ny),STAT=istatus)
wetcanp=0
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
snowdpth=0
ALLOCATE(raing(nx,ny),STAT=istatus)
raing=0
ALLOCATE(rainc(nx,ny),STAT=istatus)
rainc=0
ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
prcrate=0
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
radfrc=0
ALLOCATE(radsw(nx,ny),STAT=istatus)
radsw=0
ALLOCATE(rnflx(nx,ny),STAT=istatus)
rnflx=0
ALLOCATE(radswnet (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radswnet")
radswnet = 0
ALLOCATE(radlwin (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radlwin")
radlwin = 0
ALLOCATE(usflx(nx,ny),STAT=istatus)
usflx=0
ALLOCATE(vsflx(nx,ny),STAT=istatus)
vsflx=0
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
ptsflx=0
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
qvsflx=0
ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
uprt=0
ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
vprt=0
ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
wprt=0
ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
ptprt=0
ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
pprt=0
ALLOCATE(rhoprt(nx,ny,nz),STAT=istatus)
rhoprt=0
ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
qvprt=0
ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
tem1=0
ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
tem2=0
ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
tem3=0
ALLOCATE(lat(nx,ny),STAT=istatus)
lat=0
ALLOCATE(lon(nx,ny),STAT=istatus)
lon=0
!-------------------------------------------------------------
! Allocate for variables from remap_d.inc
!-------------------------------------------------------------
ALLOCATE (gridvel(nx,ny,nz),STAT=istatus)
gridvel=velmis
ALLOCATE (gridref(nx,ny,nz),STAT=istatus)
gridref=refmis
ALLOCATE (gridnyq(nx,ny,nz),STAT=istatus)
gridnyq=100.
ALLOCATE (gridtim(nx,ny,nz),STAT=istatus)
gridtim=0.
ALLOCATE (xsc(nx),ysc(ny),STAT=istatus)
ALLOCATE (zpsc(nx,ny,nz),STAT=istatus)
!-----------------------------------------------------------------------
!
! Get the name of the input data set.
!
!-----------------------------------------------------------------------
!
! READ(5,file_name,END=100)
! WRITE(6,'(/a)')' Namelist file_name sucessfully read in.'
!
! lengbf=LEN_trim(grdbasfn)
! lenfil=LEN_trim(filename)
!
! WRITE(6,'(/a,a)')' The data set name is ', filename(1:lenfil)
!
!------------------------------------------------------------------------
!
! Get the rada information
!
!------------------------------------------------------------------------
READ(5,radar_info,END=100)
WRITE(6,'(/a)')' Namelist radar_info sucessfully read in.'
dstmin=dstmin*1000.
dstmax=dstmax*1000.
!-----------------------------------------------------------------------
!
! Set elevation angles
!
!-----------------------------------------------------------------------
!
IF(vcpnum == 11) THEN
ntilt=14
vcptxt='Storm mode 14 tilts 0.5-19.5 deg'
DO itilt=1,ntilt
ang(itilt)=ang11(itilt)
END DO
ELSE IF (vcpnum == 12) THEN
ntilt=9
vcptxt='Storm mode 9 tilts 0.5-19.5 deg'
DO itilt=1,ntilt
ang(itilt)=ang11(itilt)
END DO
ELSE IF (vcpnum == 31) THEN
ntilt=5
vcptxt='Clear-air 5 tilts 0.5- 4.5 deg'
DO itilt=1,ntilt
ang(itilt)=ang11(itilt)
END DO
ELSE IF (vcpnum == 32) THEN
ntilt=5
vcptxt='Clear-air 5 tilts 0.5- 4.5 deg'
DO itilt=1,ntilt
ang(itilt)=ang11(itilt)
END DO
ELSE
WRITE(6,*) vcpnum,' is a bogus vcpnum'
STOP
END IF
elvmin=ang(1)-hlfwid
elvmax=ang(ntilt)+hlfwid
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hdmpfmt,nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, &
x,y,z,zp,zpsoil,uprt,vprt,wprt,ptprt,pprt, &
qvprt, qc, qr, qi, qs, qh, km, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1,tem2,tem3)
curtim = time
!
!-----------------------------------------------------------------------
!
! ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
IF( ireturn == 0 ) THEN ! successful read
!
!-----------------------------------------------------------------------
!
! Set map projection parameters
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr(mapproj,sclfct,latnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
!
gdx=x(2)-x(1)
x0sc=xctr - 0.5*(nx-3)*gdx
gdy=y(2)-y(1)
y0sc=yctr - 0.5*(ny-3)*gdy
CALL setorig(1,x0sc,y0sc)
!
!-----------------------------------------------------------------------
!
! Calculate lat,lon locations of the scalar grid points
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
xsc(i)=0.5*(x(i)+x(i+1))
END DO
xsc(nx)=2.*xsc(nx-1)-xsc(nx-2)
DO j=1,ny-1
ysc(j)=0.5*(y(j)+y(j+1))
END DO
ysc(ny)=2.*ysc(ny-1)-ysc(ny-2)
CALL xytoll(nx,ny,xsc,ysc,lat,lon)
!
!-----------------------------------------------------------------------
!
! Move z field onto the scalar grid.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zpsc(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
zpsc(i,j,nz)=(2.*zpsc(i,j,nz-1)) &
-zpsc(i,j,nz-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Combine wind perturbations and mean fields and put them
! on scalar grid.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
u(i,j,k)=0.5*(ubar(i,j,k) +uprt(i,j,k) + &
ubar(i+1,j,k)+uprt(i+1,j,k))
END DO
END DO
END DO
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
v(i,j,k)=0.5*(vbar(i,j,k) +vprt(i,j,k) + &
vbar(i,j+1,k)+vprt(i,j+1,k))
END DO
END DO
END DO
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
w(i,j,k)=0.5*(wbar(i,j,k) +wprt(i,j,k) + &
wbar(i,j,k+1)+wprt(i,j,k+1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Get reflectivity from qr.
!
!-----------------------------------------------------------------------
!
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, gridref)
!
!-----------------------------------------------------------------------
!
! Get radial velocity from 3-d velocity components
!
!-----------------------------------------------------------------------
!
knt=0
kntref=0
kntang=0
kntdst=0
DO j=2,ny-2
DO i=2,nx-2
CALL disthead(lat(i,j),lon(i,j), &
latrad,lonrad,azim,dist)
IF(dist > dstmin .AND. dist < dstmax)THEN
DO k=2,nz-2
IF(gridref(i,j,k) > thref) THEN
delz=zpsc(i,j,k)-elvrad
CALL beamelv(delz,dist,eleva,range)
IF(eleva > elvmin .AND. eleva < elvmax ) THEN
DO itilt=1,ntilt
IF(ABS(eleva-ang(itilt)) < hlfwid) GO TO 145
END DO
kntang=kntang+1
CYCLE
145 CONTINUE
knt=knt+1
CALL dhdrange(eleva,range,dhdr)
dsdr=SQRT(AMAX1(0.,(1.-dhdr*dhdr)))
CALL ddrotuv(1,lon(i,j),azim,1., &
ddrot,uazmrad,vazmrad)
gridvel(i,j,k)=u(i,j,k)*uazmrad*dsdr &
+v(i,j,k)*vazmrad*dsdr &
+w(i,j,k)*dhdr
ELSE
kntang=kntang+1
END IF
ELSE
kntref=kntref+1
END IF
END DO
ELSE
kntdst=kntdst+(nz-3)
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Write the data counters.
!
!-----------------------------------------------------------------------
!
knttry=(nx-3)*(ny-3)*(nz-3)
perc=100.*FLOAT(knt)/FLOAT(knttry)
prcdst=100.*FLOAT(kntdst)/FLOAT(knttry)
prcang=100.*FLOAT(kntang)/FLOAT(knttry)
prcref=100.*FLOAT(kntref)/FLOAT(knttry)
!
WRITE(6,'(a,i4/,a/)') ' Used VCP Number',vcpnum,vcptxt
WRITE(6,'(a,i12)') ' Total points: ',knttry
WRITE(6,'(a,i12,f6.1,a)') ' Generated data: ',knt,perc,' %'
WRITE(6,'(a,i12,f6.1,a)') ' Dist rejected: ', &
kntdst,prcdst,' %'
WRITE(6,'(a,i12,f6.1,a)') ' Angle samp rej: ', &
kntang,prcang,' %'
WRITE(6,'(a,i12,f6.1,a//)') ' Ref thresh rej: ', &
kntref,prcref,' %'
!
!-----------------------------------------------------------------------
!
! Write the radar data out.
!
!-----------------------------------------------------------------------
!
iradfmt=1
CALL ctim2abss(year,month,day,hour,minute,second, abstsec)
abstsec=abstsec+curtim
CALL abss2ctim( abstsec, iyr, imon, idy, ihr, imin, isec )
!
myr=MOD(iyr,100)
WRITE(rfname,'(a,a,i2.2,i2.2,i2.2,a,i2.2,i2.2)') &
radid,'.',myr,imon,idy,'.',ihr,imin
PRINT *, ' Writing data into ',rfname
CALL wtradcol(iradfmt,rfname,radid,latrad,lonrad,elvrad, &
myr,imon,idy,ihr,imin,isec,vcpnum, &
nx,ny,nz,gridvel,gridref,gridnyq,gridtim,xsc,ysc,zpsc)
END IF ! successful read
GOTO 110
100 CONTINUE
Write(6,'(1x,a)')'Namelist read end encountered. Program stopped.'
110 CONTINUE
STOP
END PROGRAM fakerad