PROGRAM rademul,146
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads an ARPS history file and creates data emulating that observed
! by a radar. Output is in polar (az-ran) coordinates.
!
! The initial version uses the "nearest neighbor" in the horizontal and
! linear vertical interpolation without Gaussian beam smoothing, a
! feature to be added later.
!
! Output are NetCDF files which are readable by the NSSL WDSS-II system.
!
! Radar parameters (tilt angles, beamwidth, etc) are set in the input file.
!
! AUTHOR:
! Keith Brewster, Center for Analysis and Prediction of Storms
! 02/10/2004
!
! MODIFICATIONS
! 10/30/2004 Version with beam-weight filtering. KB, CAPS
!
! 11/15/2004 Added processing for effective horizontal beam width as a
! function of rotation rate, PRT and number of samples.
! Also added an option for X-band attenuation.
! Keith Brewster, CAPS
!
! 12/05/2004 Added Gaussian random noise for velocity and reflectivity
! Sometime this could be updated to parameterize the noise
! Gaussian width as a function of SNR.
! Keith Brewster, CAPS
!
! 02/06/2005 Added processing for reading individual variables from
! ARPS hdf history files in order to save memory when reading
! very large datasets.
!
! 03/18/2005 Added loop to process mulitple files (model output times)
! in a single run.
!
! 05/10/2005 Keith Brewster, CAPS
! Added vertical vorticity variable for
! verification of vortex detection algorithms.
!
! 05/24/2005 Keith Brewster, CAPS
! Added unattenuated variables and variable write flags.
! Added smoothing option for vorticity variable.
!
! 06/14/2005 Keith Brewster, CAPS
! Modified logic to compute weighted average of data re-
! interpolated to points in a regular 7x7x7 az-ran-elv array.
!
! 12/14/2005 Keith Brewster, CAPS
! Various updates to make computing on array more flexible, add
! dual-pulse mode, add forced consistency between PRT, number
! of samples, rotation rate and output radial spacing.
!
! 04/20/2006 Keith Brewster, CAPS
! Added beam blockage.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'grid.inc'
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INTEGER :: nx, ny, nz
INTEGER :: nzsoil
INTEGER :: nstyps
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x(:) ! The x-coord. of the physical and
! computational grid.
REAL, ALLOCATABLE :: y(:) ! The y-coord. of the physical and
! computational grid.
REAL, ALLOCATABLE :: z(:) ! The z-coord. of the computational
! grid. Defined at w-point.
REAL, ALLOCATABLE :: zp(:,:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The height for soil model
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 :: qvprt (:,:,:) ! Perturbation water vapor specific
! humidity (kg/kg)
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 :: 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 :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific
! humidity (kg/kg)
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Fraction of soil types
INTEGER, ALLOCATABLE :: vegtyp(:,:) ! Vegetation type
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
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 ::tem2d(:,:) ! 2D work array
REAL, ALLOCATABLE ::tem1(:,:,:) ! Work arrays
REAL, ALLOCATABLE ::tem2(:,:,:) ! Work arrays
REAL, ALLOCATABLE ::tem3(:,:,:) ! Work arrays
REAL, ALLOCATABLE ::tem4(:,:,:) ! Work arrays
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
REAL, ALLOCATABLE :: hmax(:)
REAL, ALLOCATABLE :: hmin(:)
REAL, ALLOCATABLE :: xsc(:)
REAL, ALLOCATABLE :: ysc(:)
REAL, ALLOCATABLE :: latsc(:,:)
REAL, ALLOCATABLE :: lonsc(:,:)
REAL, ALLOCATABLE :: azmsc(:,:)
REAL, ALLOCATABLE :: sfcr(:,:)
REAL, ALLOCATABLE :: cmpref(:,:)
REAL, ALLOCATABLE :: usc(:,:,:)
REAL, ALLOCATABLE :: vsc(:,:,:)
REAL, ALLOCATABLE :: usm(:,:,:)
REAL, ALLOCATABLE :: vsm(:,:,:)
REAL, ALLOCATABLE :: wsc(:,:,:)
REAL, ALLOCATABLE :: zps(:,:,:)
REAL, ALLOCATABLE :: ref(:,:,:)
REAL, ALLOCATABLE :: refz(:,:,:)
REAL, ALLOCATABLE :: elvsc(:,:,:)
REAL, ALLOCATABLE :: rngsc(:,:,:)
REAL, ALLOCATABLE :: vort(:,:,:)
REAL, ALLOCATABLE :: vrsc(:,:,:)
!
!-----------------------------------------------------------------------
!
! NEXRAD Radar variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: maxtilt=14
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/
!
!-----------------------------------------------------------------------
!
! Grid adjustment variables -- adjustments to model grid specifications
!
!-----------------------------------------------------------------------
!
INTEGER :: adjhgt,adjctr,adjmove
REAL :: ctrlatem,ctrlonem
REAL :: umovein,vmovein
REAL :: hgtoffset
NAMELIST /grid_adj/ adjhgt,adjctr,adjmove,ctrlatem,ctrlonem, &
umovein,vmovein,hgtoffset
!
!-----------------------------------------------------------------------
!
! I/O Options
!
!-----------------------------------------------------------------------
!
INTEGER :: ifmt,creidx,ipktyp,nbits,nsmvort
INTEGER :: wrtuaref,wrtuavel,wrtvort
CHARACTER (LEN=80) :: outdir
NAMELIST /output_opts/ ifmt,creidx,ipktyp,nbits,outdir, &
nsmvort,wrtuaref,wrtuavel,wrtvort
!-----------------------------------------------------------------------
!
! Radar variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=4) :: radname
INTEGER :: locopt
REAL :: xrad,yrad,radlat,radlon,radelv
REAL :: wavelen,beamwid
REAL :: pwrxmt,gaindb,lossdb,noisedbm
NAMELIST /radar_specs/ radname,locopt, &
xrad,yrad,radlat,radlon,radelv,wavelen,beamwid, &
pwrxmt,gaindb,lossdb,noisedbm
!
!-----------------------------------------------------------------------
!
! Radar Operation Mode Variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: maxelv=30
INTEGER :: nelev,vcp,ngate
INTEGER :: unifpuls,dualprf,pulseopt,rotropt,rotdir
REAL :: gatesp,rotrate,delaz,azim1,azim2
REAL :: elev(maxelv)
REAL :: pulselen1(maxelv),pulselen2(maxelv)
REAL :: pulsewid1(maxelv),pulsewid2(maxelv)
REAL :: prt1(maxelv),prt2(maxelv)
INTEGER :: npulse1(maxelv),npulse2(maxelv)
NAMELIST /radar_opmode/ nelev,vcp,elev,ngate,gatesp, &
rotdir,rotropt,rotrate,delaz,azim1,azim2, &
unifpuls,dualprf,pulseopt, &
pulselen1,pulselen2,pulsewid1,pulsewid2, &
prt1,prt2,npulse1,npulse2
!
!-----------------------------------------------------------------------
!
! Radar emulator options
!
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: refchk = -99.
REAL, PARAMETER :: refmax = 80.
INTEGER :: nptsazm,nptselv,nptsrng
REAL :: evalwid,refmin,rngmin
REAL :: rmisval,rngfval
REAL :: hblkmin,hblkmax,gndrefl
INTEGER :: int_method,blockopt,kntrmin,kntvmin
INTEGER :: attenopt,senstvopt
INTEGER :: rferropt,vrerropt,nyqstopt
INTEGER :: tmadvopt
REAL :: sigmarf,sigmavr
REAL :: samploi,stdvrmul
REAL :: vnyquist(maxelv)
REAL :: timeincr
NAMELIST /emul_specs/ int_method,nptsazm,nptselv,nptsrng, &
evalwid,refmin,rngmin,kntrmin,kntvmin, &
blockopt,hblkmin,hblkmax,gndrefl, &
tmadvopt,attenopt,senstvopt, &
rferropt,sigmarf,vrerropt,sigmavr,samploi,stdvrmul, &
timeincr,rmisval,rngfval,nyqstopt,vnyquist
!-----------------------------------------------------------------------
!
! Radar tilt variables
!
!-----------------------------------------------------------------------
!
INTEGER :: maxazim
INTEGER :: itilt,jazim,igate,iigate
INTEGER :: iloc,jloc,iloc1,jloc1,ilc,jlc
REAL :: deg2rad,rad2deg
REAL :: irloc,jrloc
REAL :: elvang(maxelv)
REAL, ALLOCATABLE :: azim(:)
REAL, ALLOCATABLE :: beamw(:)
REAL, ALLOCATABLE :: gtspc(:)
REAL, ALLOCATABLE :: vnyq(:)
REAL, ALLOCATABLE :: radv(:,:)
REAL, ALLOCATABLE :: attradv(:,:)
REAL, ALLOCATABLE :: refl(:,:)
REAL, ALLOCATABLE :: attrefl(:,:)
REAL, ALLOCATABLE :: stdvr(:,:)
REAL, ALLOCATABLE :: vvort(:,:)
REAL, ALLOCATABLE :: sens(:)
!
!-----------------------------------------------------------------------
!
! Radar volume variables
!
!-----------------------------------------------------------------------
!
INTEGER, ALLOCATABLE :: itimvol(:)
REAL, ALLOCATABLE :: vnyqvol(:)
REAL, ALLOCATABLE :: rngvol(:,:)
REAL, ALLOCATABLE :: azmvol(:,:)
REAL, ALLOCATABLE :: elvvol(:,:)
REAL, ALLOCATABLE :: refvol(:,:,:)
REAL, ALLOCATABLE :: uarefvol(:,:,:)
REAL, ALLOCATABLE :: velvol(:,:,:)
REAL, ALLOCATABLE :: uavelvol(:,:,:)
REAL, ALLOCATABLE :: vorvol(:,:,:)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
REAL :: alatnot(2)
REAL, PARAMETER :: c = 2.9979e08 ! speed of light m/s
REAL, ALLOCATABLE :: delrng(:)
REAL, ALLOCATABLE :: delazm(:)
REAL, ALLOCATABLE :: delelv(:)
REAL, ALLOCATABLE :: wgtpt(:,:,:)
REAL, ALLOCATABLE :: blockfct(:,:)
INTEGER, PARAMETER :: nhisfile_max=1000
CHARACTER (LEN=256) :: fname
CHARACTER (LEN=256) :: fnamevol
CHARACTER (LEN=256) :: idxfname
CHARACTER (LEN=256) :: grdbasfn
CHARACTER (LEN=256) :: hisfile(nhisfile_max)
CHARACTER (LEN=80) :: vcptxt
CHARACTER (LEN=40) :: varname
CHARACTER (LEN=5) :: cmprext
INTEGER, PARAMETER :: itim1970=315619200
INTEGER :: hinfmt,nchin,nhisfile,idxunit,istatus
INTEGER :: lengbf,lenfil,itime,itimcdf,npulsmn
INTEGER :: itimtilt,iyear,imon,iday,ihour,imin,isec
INTEGER :: ifile,initime,ifsecs,ireturn
INTEGER :: i,j,k,ibgn,iend,jbgn,jend,ipt,jpt,kpt,iptmid
INTEGER :: ipmid,jpmid,kpmid,nzm2
INTEGER :: ntilt,nazim,ielv,felv,ifold,irot,ismth
INTEGER :: iradius,jradius,kntr,kntv
INTEGER :: kntgate,kntintr,kntintv,kntfewr,kntfewv
INTEGER :: kntavgr,kntavgv
INTEGER :: kntclr,knthob,kntvob,kntbob
REAL :: twdx2inv,twdxinv,twdy2inv,twdyinv
REAL :: dr,de,da,p0inv,azmsect,azmbgn,denom,ttime
REAL :: time,time0,rfrgate,frtime,radlata,radlona,hlfgtsp
REAL :: elv,azimuth,azrad,srange,hgt,hgtmsl,sfcrng
REAL :: ctrx,ctry,xsw,ysw,gtlat,gtlon,gtx,gty,xpt,ypt
REAL :: delx,dely,xrat,yrat,dhdr,dsdr,twovnyq,twovninv
REAL :: ucmp,vcmp,urad,vrad,wrad,vt,vr,vvmin,vvmax
REAL :: xrada,yrada,ipos,jpos
REAL :: w1i,w2i,w1j,w2j,w11,w12,w21,w22
REAL :: r11h,r12h,r21h,r22h,rmax
REAL :: r11l,r12l,r21l,r22l
REAL :: delv,dazm,drng,depth,elradius
REAL :: whigh,wlow,a,b,fjp1,fj0,fjm1
REAL :: flow,fhigh,reflz,refdbz,reflzgnd
REAL :: pi,cinv,b6,hlftau,hlfctau,aconst,rngg
REAL :: efbmwid,bmwrad,bmwm,sradius,vnyqstl
REAL :: fourln4,sigr2,sr2inv,wesqinv,wasqinv
REAL :: wgtg,wgtr,wgtv,sumwr,sumwv,sumv,sumr
REAL :: sumv1,sumv2,flkntv,varvr
REAL :: wgti,wgtj,wgtk
REAL :: prtmn,pulsewmn,pulselmn
REAL :: rngwid,srgpt,azmpt,azmrpt,elvpt,sfcrpt,sfcrngt,sfcrngb
REAL :: elvtop,elvbot,trnpt,trngt,topgt,tophgt,topmsl,bothgt,botmsl
REAL :: hgtpt,hgtagl,xmit,dhblkinv,blkcst,gndatten
REAL :: reflmin,reflmax
LOGICAL :: echo
!
!-----------------------------------------------------------------------
!
! External Functions
!
!-----------------------------------------------------------------------
!
REAL :: randnor
REAL :: effbmw
REAL :: erf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Initializations of constants and namelist variables
!
!-----------------------------------------------------------------------
!
itime=0
frtime=0.
time0=0.
pi=acos(-1.)
deg2rad=pi/180.
rad2deg=180./pi
cinv=1./c
p0inv=1./p0
radname='NULL'
locopt=2
xrad=0.
yrad=0.
radlat=0.
radlon=0.
radelv=0.
wavelen=10.5
beamwid=0.89
pwrxmt=475000.
gaindb=44.5
lossdb=5.
noisedbm=-113.
int_method=2
nptsazm=7
nptselv=7
nptsrng=7
evalwid=2.0
refmin=-20.
rngmin=100.
kntrmin=9
kntvmin=9
blockopt=2
hblkmin=1.0
hblkmax=10.0
gndrefl=50.0
tmadvopt=1
attenopt=0
senstvopt=0
rferropt=1
sigmarf=5.0
vrerropt=1
sigmavr=1.0
samploi=0.5
stdvrmul=1.5
timeincr=1.
rmisval=-99900.
rngfval=-99901.
nyqstopt=1
vnyquist=100.
nelev=1
vcp=0
elev=0.5
ngate=100
gatesp=1000.
rotdir=1
rotropt=1
rotrate=18.0
delaz=1.0
azim1=0.
azim2=360.
unifpuls=1
dualprf=0
pulseopt=1
pulselen1=1.57E-06
pulselen2=1.57E-06
pulsewid1=441.
pulsewid2=441.
prt1=1.06e-03
prt2=1.06e-03
npulse1=54
npulse2=54
ifmt=1
creidx=0
ipktyp=2
nbits=16
nsmvort=2
wrtuaref=0
wrtuavel=0
wrtvort =0
outdir='./'
adjhgt=0
adjctr=0
adjmove=0
ctrlatem=45.
ctrlonem=-90.
umovein=0.
vmovein=0.
hgtoffset=0.
irloc=45.
jrloc=45.
cmprext='.gz'
!
!-----------------------------------------------------------------------
!
! Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
WRITE(6,'(a,i4,a/)') ' Processing ',nhisfile,' model files.'
lengbf = len_trim(grdbasfn)
!
!-----------------------------------------------------------------------
!
! Read other input namelists
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a)') ' Reading grid_adj namelist'
READ(5,grid_adj)
WRITE(6,'(a)') ' Reading output_opts namelist'
READ(5,output_opts)
WRITE(6,'(a)') ' Reading radar_specs namelist'
READ(5,radar_specs)
WRITE(6,'(a)') ' Reading radar_opmode namelist'
READ(5,radar_opmode)
WRITE(6,'(a)') ' Reading emul_specs namelist'
READ(5,emul_specs)
IF(vcp == 0) THEN
ntilt=nelev
write(vcptxt,'(a,i2,a,f4.1,a,f4.1,a)') &
'Custom vcp',ntilt,' tilts',elev(1),'-',elev(ntilt),' deg'
vcptxt='Custom vcp 14 tilts 0.5-19.5 deg'
DO itilt=1,ntilt
elvang(itilt)=elev(itilt)
END DO
ELSE IF(vcp == 11) THEN
ntilt=14
vcptxt='Storm mode 14 tilts 0.5-19.5 deg'
DO itilt=1,ntilt
elvang(itilt)=ang11(itilt)
END DO
ELSE IF (vcp == 12) THEN
ntilt=9
vcptxt='Storm mode 9 tilts 0.5-19.5 deg'
DO itilt=1,ntilt
elvang(itilt)=ang11(itilt)
END DO
ELSE IF (vcp == 31) THEN
ntilt=5
vcptxt='Clear-air 5 tilts 0.5- 4.5 deg'
DO itilt=1,ntilt
elvang(itilt)=ang11(itilt)
END DO
ELSE IF (vcp == 32) THEN
ntilt=5
vcptxt='Clear-air 5 tilts 0.5- 4.5 deg'
DO itilt=1,ntilt
elvang(itilt)=ang11(itilt)
END DO
ELSE
WRITE(6,*) vcp,' is not a recognized VCP number'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Print a warning about blocking/clutter using integration method 1.
!
!-----------------------------------------------------------------------
!
IF(int_method == 1 .AND. blockopt > 0) THEN
WRITE(6,'(//a,/a,i4,/a,i4,/a//)') &
' ***********WARNING*********************', &
' Terrain blocking calculation not possible using integration method ', &
int_method, ' Select integration method 2 with blockopt=',blockopt,&
'Continuing using integration method 2...'
int_method=2
END IF
!
!-----------------------------------------------------------------------
!
! Calculate constants for antenna gain weighting.
!
!-----------------------------------------------------------------------
!
kntvmin=max(kntvmin,2)
fourln4=4.0*alog(4.0)
wesqinv=fourln4/(beamwid*beamwid)
reflzgnd=10.**(gndrefl*0.10)
IF(azim2 > azim1) THEN
azmsect=azim2-azim1
ELSE
azmsect=(360.+azim2)-azim1
END IF
IF(unifpuls == 1) THEN
DO itilt=2,ntilt
prt1(itilt)=prt1(1)
prt2(itilt)=prt2(1)
npulse1(itilt)=npulse1(1)
npulse2(itilt)=npulse2(1)
pulselen1(itilt)=pulselen1(1)
pulselen2(itilt)=pulselen2(1)
pulsewid1(itilt)=pulsewid1(1)
pulsewid2(itilt)=pulsewid2(1)
END DO
END IF
IF(rotropt > 1) THEN
nazim=NINT(azmsect/delaz)+1
maxazim=max(nazim,1)
ELSE
maxazim=1
DO itilt=1,ntilt
IF(dualprf > 0) THEN
delaz=rotrate*((prt1(itilt)*npulse1(itilt))+ &
(prt2(itilt)*npulse2(itilt)))
ELSE
delaz=rotrate*prt1(itilt)*npulse1(itilt)
END IF
nazim=NINT(azmsect/delaz)+1
maxazim=max(nazim,maxazim)
END DO
END IF
WRITE(6,'(a,i6)') ' Maximum azimuths: ',maxazim
rfrgate=gatesp
hlfgtsp=0.5*gatesp
IF(blockopt == 0) THEN
hblkmin=0.
hblkmax=1.0
dhblkinv=1.0
ELSE
hblkmin=max(hblkmin,0.0)
hblkmax=max(hblkmax,(hblkmin+1.0))
dhblkinv=1./(hblkmax-hblkmin)
END IF
nptsrng=max(nptsrng,5)
ALLOCATE (delrng(nptsrng))
dr=gatesp/float(nptsrng-1)
WRITE(6,'(a,f10.2,a)') ' Range evaluation points (gatesp=',gatesp,'):'
DO ipt=1,nptsrng
delrng(ipt)=(-0.5*gatesp)+((ipt-1)*dr)
WRITE(6,'(a,i4,a,f12.4)') ' delrng(',ipt,') = ',delrng(ipt)
END DO
iptmid=(nptsrng/2)+1
blkcst=(0.001*gatesp)/float(nptsrng)
nptsazm=max(nptsazm,5)
ALLOCATE (delazm(nptsazm))
nptselv=max(nptselv,5)
ALLOCATE (delelv(nptselv))
elradius=evalwid*beamwid
de=elradius/float(nptselv-1)
WRITE(6,'(a,f10.2,a,f10.2)') &
' Beam width ',beamwid,' Eval region width:',(evalwid*beamwid)
WRITE(6,'(a)') &
' Vertical cross-beam evaluation pts (elev deg above center):'
DO kpt=1,nptselv
delelv(kpt)=(-0.5*evalwid*beamwid)+((kpt-1)*de)
WRITE(6,'(a,i4,a,f12.4)') ' delelv(',kpt,') = ',delelv(kpt)
END DO
ipmid=(nptsrng/2)+1
jpmid=(nptsazm/2)+1
kpmid=(nptselv/2)+1
ALLOCATE (wgtpt(nptsrng,nptsazm,nptselv))
ALLOCATE (blockfct(nptsazm,nptselv))
!
!-----------------------------------------------------------------------
!
! Obtain the grid dimensions from input data.
!
!-----------------------------------------------------------------------
!
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps,ireturn)
nstyp = nstyps ! Copy to global variable
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,', nzsoil',nzsoil
nzm2=nz-2
!-----------------------------------------------------------------------
!
! Allocate grid position arrays.
!
!-----------------------------------------------------------------------
!
allocate(x(nx),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:x")
x=0.
allocate(y(ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:y")
y=0.
allocate(z(nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:z")
z=0.
allocate(zp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:zp")
zp=0.
DO ifile=1,nhisfile
WRITE(6,'(a,a)') ' Processing file: ',TRIM(hisfile(ifile))
IF( hinfmt .NE. 3) THEN
!-----------------------------------------------------------------------
!
! Allocate arrays needed for dtaread
! These are deallocated after being read and used so they need
! to be reallocated with each new file.
!
!-----------------------------------------------------------------------
!
allocate(zpsoil (nx,ny,nzsoil),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:zpsoil")
zpsoil=0.
allocate(uprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:uprt")
uprt=0.
allocate(vprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vprt")
vprt=0.
allocate(wprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:wprt")
wprt=0.
allocate(ptprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:ptprt")
ptprt=0.
allocate(pprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:pprt")
pprt=0.
allocate(qvprt (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qvprt")
qvprt=0.
allocate(ubar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:ubar")
ubar=0.
allocate(vbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vbar")
vbar=0.
allocate(wbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:wbar")
wbar=0.
allocate(ptbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:ptbar")
ptbar=0.
allocate(pbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:pbar")
pbar=0.
allocate(rhobar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:rhobar")
rhobar=0.
allocate(qvbar (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qvbar")
qvbar=0.
allocate(qc (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qc")
qc=0.
allocate(qr (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qr")
qr=0.
allocate(qi (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qi")
qi=0.
allocate(qs (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qs")
qs=0.
allocate(qh (nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qh")
qh=0.
allocate(soiltyp (nx,ny,nstyps),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:soiltyp")
soiltyp=0
allocate(stypfrct(nx,ny,nstyps),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:stypfrct")
stypfrct=0.
allocate(vegtyp(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vegtyp")
vegtyp=0
allocate(tsoil (nx,ny,nzsoil,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tsoil")
tsoil=0.
allocate(qsoil (nx,ny,nzsoil,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:qsoil")
qsoil=0.
allocate(wetcanp(nx,ny,0:nstyps),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:wetcanp")
wetcanp=0.
allocate(prcrate(nx,ny,4),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:prcrate")
prcrate=0.
allocate(tem2d(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tem2d")
tem2d=0.
allocate(tem1(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tem1")
tem1=0.
allocate(tem2(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tem2")
tem2=0.
allocate(tem3(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tem3")
tem3=0.
allocate(tem4(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:tem4")
tem4=0.
lenfil = len_trim(hisfile(ifile))
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
hisfile(ifile),lenfil,time, &
x,y,z,zp,zpsoil,uprt,vprt,wprt,ptprt,pprt, &
qvprt, qc, qr, qi, qs, qh, tem4,tem4,tem4, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,tem2d,tem2d,tem2d, &
tsoil,qsoil,wetcanp,tem2d, &
tem2d,tem2d,prcrate, &
tem4,tem2d,tem2d,tem2d,tem2d, &
tem2d,tem2d,tem2d,tem2d, &
ireturn, tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! Deallocate the no-longer-needed dtaread arrays.
!
!-----------------------------------------------------------------------
!
deallocate(qvprt)
deallocate(qvbar)
deallocate(wbar)
deallocate(qc)
deallocate(qi)
deallocate(zpsoil)
deallocate(soiltyp)
deallocate(stypfrct)
deallocate(vegtyp)
deallocate(tsoil)
deallocate(qsoil)
deallocate(wetcanp)
deallocate(prcrate)
deallocate(tem2d)
deallocate(tem2)
deallocate(tem3)
deallocate(tem4)
allocate(usc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:usc")
usc=0.
allocate(vsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vsc")
vsc=0.
allocate(wsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:wsc")
wsc=0.
dxinv=1./dx
dyinv=1./dx
twdx2inv=1./(2.*dx*dx)
twdxinv=1./(2.*dx)
twdy2inv=1./(2.*dy*dy)
twdyinv=1./(2.*dy)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
usc(i,j,k)=0.5*((uprt( i,j,k)+ubar( i,j,k))+ &
(uprt(i+1,j,k)+ubar(i+1,j,k)))+umove
vsc(i,j,k)=0.5*((vprt(i, j,k)+vbar(i, j,k))+ &
(vprt(i,j+1,k)+vbar(i,j+1,k)))+vmove
wsc(i,j,k)=0.5*(wprt(i,j,k)+wprt(i,j,k+1))
END DO
END DO
END DO
deallocate(ubar)
deallocate(uprt)
deallocate(vbar)
deallocate(vprt)
deallocate(wprt)
!
!-----------------------------------------------------------------------
!
! Calculate reflectivity at all scalar grid points.
! First need Temperature in K, which is stored in tem1
!
!-----------------------------------------------------------------------
!
tem1=0.
tem2=0.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=(ptprt(i,j,k)+ptbar(i,j,k))* &
(((pprt(i,j,k) + pbar(i,j,k))*p0inv) ** rddcp)
END DO
END DO
END DO
deallocate(pbar)
deallocate(pprt)
deallocate(ptbar)
deallocate(ptprt)
allocate(ref(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:ref")
ref=0.
! CALL reflec_ferrier(nx,ny,nz, tem2, qr, qs, qh, tem1, ref)
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, ref)
!
deallocate(qr)
deallocate(qs)
deallocate(qh)
deallocate(tem1)
deallocate(rhobar)
ELSE ! hdf format
allocate(itmp(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:itmp")
itmp=0
allocate(hmax(nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:hmax")
hmax=0.
allocate(hmin(nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:hmin")
hmin=0.
IF(ifile == 1) THEN
CALL get_gridxyzp_hdf(nx,ny,nz,grdbasfn(1:lengbf), &
x,y,zp, &
itmp,hmax,hmin,ireturn)
IF( ireturn /= 0 ) THEN
WRITE(6,'(a,a)') ' Problem getting x,y,zp from hdf file ', &
grdbasfn(1:lengbf)
WRITE(6,'(a)') ' Program stopped.'
STOP
END IF
print *,' HDF info: x,y,z of input data read in.'
print *,' x(1 )=',x(1)
print *,' x(nx)=',x(nx)
print *,' y(1 )=',y(1)
print *,' y(ny)=',y(ny)
END IF
!
!-----------------------------------------------------------------------
!
! Read in the wind data and average to the scalar points
! Here the tem arrays are used to store the total wind components
! Note also: umove and vmove are returned via common
!
!-----------------------------------------------------------------------
!
allocate(tem1(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:tem1")
tem1=0.
allocate(tem2(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:tem2")
tem2=0.
allocate(tem3(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:tem3")
tem3=0.
allocate(tem4(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:tem4")
tem4=0.
allocate(usc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:usc")
usc=0.
allocate(vsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:vsc")
vsc=0.
allocate(wsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:wsc")
wsc=0.
allocate(rhobar(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:rhobar")
rhobar=0.
CALL hdfreaduvw(nx,ny,nz,TRIM(hisfile(ifile)), &
time,tem1,tem2,tem3, &
itmp,hmax,hmin,ireturn)
IF(ireturn < 0) EXIT
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
usc(i,j,k)=(0.5*(tem1(i,j,k)+tem1(i+1,j,k)))+umove
vsc(i,j,k)=(0.5*(tem2(i,j,k)+tem2(i+1,j,k)))+vmove
wsc(i,j,k)=0.5*(tem3(i,j,k)+tem3(i,j,k+1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Get pressure and potential temperature to find temperature in K
! Here the tem1 = total pressure and tem2 is total potential temp
! Use the rhobar array to store rho.
!
!-----------------------------------------------------------------------
!
tem1=0.
tem2=0.
CALL hdfreadppt(nx,ny,nz,TRIM(hisfile(ifile)), &
time,tem1,tem2, &
itmp,hmax,hmin,ireturn)
IF(ireturn < 0) EXIT
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k)=(tem2(i,j,k))*((tem1(i,j,k)*p0inv) ** rddcp)
rhobar(i,j,k)=tem1(i,j,k)/(rd*tem4(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Get the hydrometeors qr,qs, and qh in tem1,tem2, and tem3
!
!-----------------------------------------------------------------------
!
tem1=0.
tem2=0.
tem3=0.
CALL hdfreadqrsh(nx,ny,nz,TRIM(hisfile(ifile)), &
time,tem1,tem2,tem3, &
itmp,hmax,hmin,ireturn)
IF(ireturn < 0) EXIT
!
!-----------------------------------------------------------------------
!
! Free the hdf io arrays
!
!-----------------------------------------------------------------------
!
deallocate(itmp)
deallocate(hmax)
deallocate(hmin)
!
!-----------------------------------------------------------------------
!
! Compute the reflectivity
!
!-----------------------------------------------------------------------
!
allocate(ref(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:hdf:ref")
ref=0.
! CALL reflec_ferrier(nx,ny,nz, rhobar, tem1, tem2, tem3, tem4, ref)
CALL reflec(nx,ny,nz, rhobar, tem1, tem2, tem3, ref)
!
!-----------------------------------------------------------------------
!
! Free the temporary i/o arrays
!
!-----------------------------------------------------------------------
!
deallocate(tem1)
deallocate(tem2)
deallocate(tem3)
deallocate(tem4)
deallocate(rhobar)
END IF ! hdf format
!
!-----------------------------------------------------------------------
!
! Set up fake reflectivity for testing terrain blocking.
!
!-----------------------------------------------------------------------
! CALL fake_reflec(nx,ny,nz,ref)
! CALL fake_vel(nx,ny,nz,usc,vsc)
!
!-----------------------------------------------------------------------
!
! Make a few calculations now that we have read the datafile.
!
!-----------------------------------------------------------------------
!
dx=x(2)-x(1)
dy=y(2)-y(1)
dxinv=1./dx
dyinv=1./dy
curtim = time
IF(ifile == 1 .AND. adjctr == 2) time0=curtim
CALL ctim2abss(year,month,day,hour,minute,second,itime)
ifsecs=NINT(curtim)
itime=itime+ifsecs
itimtilt=itime
initime=itimcdf-ifsecs
CALL gtlfnkey(runname,lfnkey)
print *, ' yr,mo,day,hr,min,sec:',year,month,day
print *, ' hr,min,sec:',hour,minute,second
print *, ' ifsecs: ',ifsecs
print *, ' curtim: ',curtim
print *, ' mapproj,ctrlat,ctrlon: ',mapproj,ctrlat,ctrlon
CALL abss2ctim(itime,year,month,day,hour,minute,second)
!
!-----------------------------------------------------------------------
!
! Average spatial locations to scalar points
! Need only do this once, for the first file.
!
!-----------------------------------------------------------------------
!
IF(ifile == 1) THEN
allocate(xsc(nx),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:xsc")
xsc=0.
allocate(ysc(ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:ysc")
ysc=0.
allocate(zps(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:zps")
zps=0.
DO i=1,nx-1
xsc(i)=0.5*(x(i)+x(i+1))
END DO
xsc(nx)=xsc(nx-1)+dx
DO j=1,ny-1
ysc(j)=0.5*(y(j)+y(j+1))
END DO
ysc(ny)=ysc(ny-1)+dy
write(6,'(a,f10.2,a,f10.2,a)') ' Coordinate of x(2): ',(x(2)*0.001), &
' km y(2): ',(y(2)*0.001),' km'
IF( adjhgt > 0) THEN
write(6,'(a,f10.2,a)') &
' Adding height offset to grid:',hgtoffset,' m.'
ELSE
hgtoffset=0.
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zps(i,j,k)=(0.5*(zp(i,j,k)+zp(i,j,k+1)))+hgtoffset
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
zps(i,j,nz)=(2.0*zps(i,j,nz-1))-zps(i,j,nz-2)
END DO
END DO
END IF ! first file
!
!-----------------------------------------------------------------------
!
! Allocate additional grid location and geometry arrays
! These speed calculations down the road.
!
!-----------------------------------------------------------------------
!
allocate(latsc(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:latsc")
latsc=0.
allocate(lonsc(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:lonsc")
lonsc=0.
allocate(azmsc(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:azmsc")
azmsc=0.
allocate(sfcr(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:sfcr")
sfcr=0.
allocate(elvsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:elvsc")
elvsc=0.
allocate(rngsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:rngsc")
rngsc=0.
allocate(vrsc(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vrsc")
vrsc=0.
allocate(usm(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:usm")
usm=0.
allocate(vsm(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vsm")
vsm=0.
allocate(vort(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vort")
vort=0.
allocate(refz(nx,ny,nz),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:refz")
refz=0.
allocate(cmpref(nx,ny),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:cmpref")
cmpref=0.
!
!-----------------------------------------------------------------------
!
! Define output data size and find location of radar in the grid.
! Account for umove and vmove by moving the radar relative to the grid.
!
! If no map projection has been defined, for the purposes of emulation,
! set the projection to Lambert conformal and center the grid on
! ctrlat, ctrlon.
!
!-----------------------------------------------------------------------
!
IF ( adjmove > 0 ) THEN
umove=umovein
vmove=vmovein
END IF
IF ( ifile == 1 ) THEN
IF ( mapproj == 0 ) THEN
IF ( adjctr > 0 ) THEN
ctrlat = ctrlatem
ctrlon = ctrlonem
END IF
print *, ' ctrlat, ctrlon: ',ctrlat,ctrlon
trulat1 = ctrlat
trulat2 = ctrlat
trulon = ctrlon
alatnot(1)=trulat1
alatnot(2)=trulat2
CALL setmapr(2,1.0,alatnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,ctrx,ctry)
xsw=ctrx-0.5*((nx-3)*dx)
ysw=ctry-0.5*((ny-3)*dy)
CALL setorig(1,xsw,ysw)
ELSE IF ( adjctr > 0 ) THEN
ctrlat = ctrlatem
ctrlon = ctrlonem
alatnot(1)=trulat1
alatnot(2)=trulat2
CALL setmapr(mapproj,1.0,alatnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,ctrx,ctry)
xsw=ctrx-0.5*((nx-3)*dx)
ysw=ctry-0.5*((ny-3)*dy)
CALL setorig(1,xsw,ysw)
ELSE
alatnot(1)=trulat1
alatnot(2)=trulat2
CALL setmapr(mapproj,1.0,alatnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,ctrx,ctry)
xsw=ctrx-0.5*((nx-3)*dx)
ysw=ctry-0.5*((ny-3)*dy)
CALL setorig(1,xsw,ysw)
END IF
END IF
IF(locopt == 1) THEN
print *, ' radlat,radlon=',radlat,radlon
CALL lltoxy(1,1,radlat,radlon,xrad,yrad)
xrada=xrad-umove*(curtim-time0)
yrada=yrad-vmove*(curtim-time0)
CALL xytoll(1,1,xrada,yrada,radlata,radlona)
irloc=((xrada-xsc(1))*dxinv)+1.0
jrloc=((yrada-ysc(1))*dyinv)+1.0
WRITE(6,'(a,/a,f10.4,a,f10.4,/a,f10.4,a,f10.4,/a,f10.4,a,f10.4)') &
' Radar position (adjusted for umove, vmove)', &
' Radar latitude:',radlata,' longitude: ',radlona, &
' x-coord (km):',(0.001*xrada),' y-coord (km):',(0.001*yrada), &
' i-location: ',irloc,' j-location:',jrloc
ELSE
xrada=xrad*1000.
yrada=yrad*1000.
CALL xytoll(1,1,xrada,yrada,radlat,radlon)
radlata=radlat
radlona=radlon
irloc=((xrada-xsc(1))*dxinv)+1.0
jrloc=((yrada-ysc(1))*dyinv)+1.0
WRITE(6,'(a,f10.4,a,f10.4,/a,f10.4,a,f10.4,/a,f10.4,a,f10.4)') &
' Radar latitude:',radlata,' longitude: ',radlona, &
' x-coord (km):',(0.001*xrada),' y-coord (km):',(0.001*yrada),&
' i-location: ',irloc,' j-location:',jrloc
END IF
iloc=INT(irloc)
jloc=INT(jrloc)
trngt=0.
IF(iloc > 0 .AND. iloc < nx-2 .AND. &
jloc > 0 .AND. jloc < ny-2 ) THEN
w2i=irloc-iloc
w1i=1.0-w2i
w2j=jrloc-jloc
w1j=1.0-w2j
w11=w1i*w1j
w21=w2i*w1j
w12=w1i*w2j
w22=w2i*w2j
trngt=w11*zp(iloc,jloc,2) + &
w21*zp(iloc+1,jloc,2) + &
w12*zp(iloc,jloc+1,2) + &
w22*zp(iloc+1,jloc+1,2)
END IF
WRITE(6,'(a,f8.1,a,f8.1//)') &
' Radar Elevation: ',radelv,' Terrain Height= ',trngt
!
!-----------------------------------------------------------------------
!
! Calculate radar parameters at each scalar grid point.
!
!-----------------------------------------------------------------------
print *, ' calculating radar parameters'
!
IF( mapproj > 0 ) THEN
CALL xytoll(nx,ny,xsc,ysc,latsc,lonsc)
DO j=1,ny-1
DO i=1,nx-1
CALL disthead(radlata,radlona,latsc(i,j),lonsc(i,j), &
azmsc(i,j),sfcr(i,j))
IF(sfcr(i,j) > rngmin) THEN
DO k=1,nz-1
CALL beamelv((zps(i,j,k)-radelv),sfcr(i,j), &
elvsc(i,j,k),rngsc(i,j,k))
END DO
ELSE
rngsc(i,j,k)=sfcr(i,j)
elvsc(i,j,k)=-99.
END IF
END DO
END DO
ELSE
DO j=1,ny-1
DO i=1,nx-1
delx=xsc(i)-xrada
dely=ysc(j)-yrada
sfcr(i,j)=sqrt(delx*delx + dely*dely)
IF(dely > 0.) THEN
azmsc(i,j)=rad2deg*atan(delx/dely)
ELSE IF (dely < 0.) THEN
azmsc(i,j)=180.+rad2deg*atan(delx/dely)
ELSE IF (delx > 0.) THEN
azmsc(i,j)=90.
ELSE
azmsc(i,j)=270.
END IF
IF(azmsc(i,j) < 0.) azmsc(i,j)=azmsc(i,j)+360.
IF(sfcr(i,j) > rngmin) THEN
DO k=1,nz-1
CALL beamelv((zps(i,j,k)-radelv),sfcr(i,j), &
elvsc(i,j,k),rngsc(i,j,k))
END DO
ELSE
rngsc(i,j,k)=sfcr(i,j)
elvsc(i,j,k)=-99.
END IF
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Calculate radial velocity, including terminal velocity, and reflectivity
! factor (in mm6/m3) at all scalar grid points.
!
! The termvel subroutine can be updated sometime to use actual freezing
! level and/or the model species variables to calculate vt.
!
!-----------------------------------------------------------------------
!
print *, ' Calculating radial velocity'
DO j=1,ny-1
DO i=1,nx-1
ucmp=sin(deg2rad*azmsc(i,j))
vcmp=cos(deg2rad*azmsc(i,j))
cmpref(i,j)=-99.
DO k=2,nz-1
CALL termvel(ref(i,j,k),zps(i,j,k),vt)
IF(ref(i,j,k) > 0.) THEN
refz(i,j,k)=10.**(0.1*ref(i,j,k))
ELSE
refz(i,j,k)=0.
END IF
CALL dhdrange(elvsc(i,j,k),rngsc(i,j,k),dhdr)
dsdr=sqrt(1.-dhdr*dhdr)
vrsc(i,j,k)=dsdr*(ucmp*usc(i,j,k) + vcmp*vsc(i,j,k)) + &
dhdr*(wsc(i,j,k)-vt)
cmpref(i,j)=max(cmpref(i,j),ref(i,j,k))
END DO
END DO
END DO
usm=usc
vsm=usc
print *, ' Smoothing velocities, nsmvort=',nsmvort
IF(nsmvort > 0) THEN
DO ismth=1,nsmvort
CALL smooth9p( usm, nx,ny, 1, nx-1, 1, ny-1, 0, vort )
CALL smooth9p( vsm, nx,ny, 1, nx-1, 1, ny-1, 0, vort )
END DO
vort=0.
END IF
!
! Calculate vertical vorticity at scalar points
!
WRITE(6,'(a,f10.8,a,f10.8)') ' Calculating vorticity: dxinv=', &
dxinv,' dyinv: ',dyinv
DO k=1,nz
DO j=2,ny-2
DO i=2,nx-2
vort(i,j,k)=0.5*((vsm(i+1,j,k)-vsm(i-1,j,k))*dxinv - &
(usm(i,j+1,k)-usm(i,j-1,k))*dyinv)
END DO
END DO
END DO
DO k=1,nz
DO j=2,ny-2
vort(1,j,k)=vort(2,j,k)
vort(nx-1,j,k)=vort(nx-2,j,k)
END DO
END DO
DO k=1,nz
DO i=1,nx-1
vort(i,1,k)=vort(1,2,k)
vort(i,ny-1,k)=vort(i,ny-2,k)
END DO
END DO
vvmax=maxval(vort)
vvmin=minval(vort)
WRITE(6,'(a,g16.9,a,g16.9)') ' Grid min vert vort=',vvmin, &
' Max vert vort = ',vvmax
!
deallocate(usm)
deallocate(vsm)
!
! Allocate output radar arrays
!
allocate(azim(maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:azim")
azim=0.
allocate(beamw(maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:beamw")
beamw=0.
allocate(gtspc(maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:gtspc")
gtspc=0.
allocate(vnyq(maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vnyq")
vnyq=0.
allocate(refl(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:refl")
refl=0.
allocate(attrefl(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:attrefl")
attrefl=0.
allocate(radv(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:radv")
radv=0.
allocate(attradv(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:attradv")
attradv=0.
allocate(stdvr(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:stdvr")
stdvr=0.
allocate(vvort(ngate,maxazim),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vvort")
vvort=0.
allocate(sens(ngate),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:sens")
sens=0.
IF( ifmt == 1 ) THEN
allocate(itimvol(nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:itimvol")
itimvol=0
allocate(vnyqvol(nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vnyqvol")
vnyqvol=0.
allocate(rngvol(ngate,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:rngvol")
rngvol=0.
allocate(azmvol(maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:azmvol")
azmvol=0.
allocate(elvvol(maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:elvvol")
elvvol=0.
allocate(refvol(ngate,maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:refvol")
refvol=0.
allocate(uarefvol(ngate,maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:uarefvol")
uarefvol=0.
allocate(velvol(ngate,maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:velvol")
velvol=0.
allocate(uavelvol(ngate,maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:uavelvol")
uavelvol=0.
allocate(vorvol(ngate,maxazim,nelev),stat=istatus)
CALL check_alloc_status(istatus, "radaremul:vorvol")
vorvol=0.
END IF
!
!-----------------------------------------------------------------------
!
! Open log file for index records
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 .AND. ifmt == 2) THEN
write(idxfname,'(a,a,a,a,a,a,i6.6,a)') &
TRIM(outdir),'/',radname,'.',runname(1:lfnkey),'.',ifsecs,'.xml'
CALL getunit(idxunit)
write(6,'(a,a)') ' Writing index to:',TRIM(idxfname)
OPEN(idxunit,file=idxfname,form='formatted',status='unknown')
write(idxunit,'(a)') &
'<?xml version="1.0" encoding="iso-8859-1" ?>'
write(idxunit,'(a,a,a)') '<codeindex type="netcdf" dataset="', &
TRIM(outdir),'">'
END IF
!
!-----------------------------------------------------------------------
!
! Create data from Cartesian grid.
!
!-----------------------------------------------------------------------
!
irot=1
IF(rotdir == 2) irot=-1
print *, ' entering tilt loop'
DO itilt=1,ntilt
kntgate=0
kntintr=0
kntintv=0
kntfewr=0
kntfewv=0
kntclr=0
knthob=0
kntvob=0
kntbob=0
write(6,'(/a,i3,a,f6.2,a)') &
' Processing tilt ',itilt,' elev=',elvang(itilt),' deg'
azim=rmisval
beamw=rmisval
gtspc=rmisval
vnyq=rmisval
refl=rmisval
attrefl=rmisval
radv=rmisval
attradv=rmisval
vvort=0.
vvmax=-999.
vvmin=999.
elv=elvang(itilt)
ielv=INT(elv)
felv=NINT(100.*(elv-ielv))
elvtop=elv+delelv(nptselv)
elvbot=elv+delelv(1)
write(6,'(a,f6.2,a,f6.2)') &
' Elv ang beam bot=',elvbot,' Elv ang beam top=',elvtop
IF(tmadvopt == 1) THEN
itimtilt=itime+NINT((itilt-1)*timeincr)
END IF
print *, ' itilt, itimtilt, tmadvopt, timeincr: ', &
itilt,itimtilt,tmadvopt,timeincr
IF(rotropt > 1) THEN
nazim=NINT(azmsect/delaz)+1
IF(dualprf > 0) THEN
ttime=nazim*((prt1(itilt)*npulse1(itilt))+ &
(prt2(itilt)*npulse2(itilt)))
ELSE
ttime=nazim*prt1(itilt)*npulse1(itilt)
END IF
rotrate=azmsect/ttime
ELSE
IF(dualprf > 0) THEN
delaz=rotrate*((prt1(itilt)*npulse1(itilt))+ &
(prt2(itilt)*npulse2(itilt)))
ELSE
delaz=rotrate*prt1(itilt)*npulse1(itilt)
END IF
nazim=NINT(azmsect/delaz)+1
END IF
WRITE(6,'(a,f10.2,a,/a,f10.2,a,/a,i6)') &
' Rotation rate:',rotrate,' deg/s', &
' Azimuth sampling increment:',delaz,' deg', &
' Number of azimuth samples:',nazim
!
!-----------------------------------------------------------------------
!
! Calculate constants for range weighting.
! Assumed here is a square wave transmitted with a matched filter
! in the receiver resulting in a received wave envelope that is
! approximately Gaussian.
! See Eqs 11.118 and 5.76 in Doviak and Zrnic', 2nd Ed, 1993.
!
! Note that pulselen has units of time, pulsewid units of meters
!
!-----------------------------------------------------------------------
!
IF(dualprf == 0) THEN
IF(pulseopt == 1) THEN
pulselmn=pulselen1(itilt)
pulsewmn=c*pulselmn
ELSE
pulsewmn=pulsewid1(itilt)
pulselmn=cinv*pulsewmn
END IF
prtmn=prt1(itilt)
npulsmn=npulse1(itilt)
IF(nyqstopt == 1) THEN
vnyqstl=vnyquist(itilt)
ELSE
vnyqstl=0.25*(0.01*wavelen)/prt1(itilt)
vnyquist(itilt)=vnyqstl
END IF
ELSE
IF(pulseopt == 1) THEN
pulselmn=0.5*(pulselen1(itilt)*pulselen2(itilt))
pulsewmn=c*pulselmn
ELSE
pulsewmn=0.5*(pulsewid1(itilt)*pulsewid2(itilt))
pulselmn=cinv*pulsewmn
END IF
prtmn=0.5*(prt1(itilt)+prt2(itilt))
npulsmn=NINT(0.5*(npulse1(itilt)+npulse2(itilt)))
IF(nyqstopt == 1) THEN
vnyqstl=vnyquist(itilt)
ELSE
vnyqstl=abs(0.25*(0.01*wavelen)/(prt1(itilt)-prt2(itilt)))
vnyquist(itilt)=vnyqstl
END IF
END IF
hlftau=0.5*pulselmn
aconst=pi/(2.*sqrt(log(2.0)))
b6=1.04/pulselmn
hlfctau=0.5*pulsewmn
sigr2=(0.35*hlfctau)*(0.35*hlfctau)
sr2inv=2.0/(4.0*sigr2)
rngwid=sqrt(sigr2)
twovnyq=2.*vnyqstl
twovninv=1./twovnyq
WRITE(6,'(a,f10.2,a)') ' Mean PRT: ',(prtmn*1000.),' millisec'
WRITE(6,'(a,f10.2,a)') ' Mean pulse length: ',(pulselmn*1.0E06), &
' microsec'
WRITE(6,'(a,f10.2,a)') ' Mean pulse width: ',pulsewmn,' m'
WRITE(6,'(a,f10.2,a)') &
' Gaussian filtered pulse half-width: ',rngwid,' m'
WRITE(6,'(a,f10.2,a)') ' Unambiguous velocity: ',vnyqstl,' m/s'
print *, 'rotrate,prtmn,npulsmn: ',rotrate,prtmn,npulsmn
IF( abs(rotrate) > 0.) THEN
efbmwid=EFFBMW(beamwid,rotrate,prtmn,npulsmn)
ELSE
efbmwid=beamwid
END IF
WRITE(6,'(a,f8.2,a)') ' Static beamwidth:',beamwid,' degrees.'
WRITE(6,'(a,f8.2,a)') ' Effective beamwidth:',efbmwid,' degrees.'
da=(evalwid*efbmwid)/float(nptsazm-1)
WRITE(6,'(a)') &
' Horiz cross-beam evaluation pts (azimuth degrees from center):'
DO jpt=1,nptsazm
delazm(jpt)=(-0.5*evalwid*efbmwid)+((jpt-1)*da)
WRITE(6,'(a,i4,a,f12.4)') ' delazm(',jpt,') = ',delazm(jpt)
END DO
!
!-----------------------------------------------------------------------
!
! Calculate sensitivity as a function of range.
!
!-----------------------------------------------------------------------
!
IF(senstvopt == 1) CALL xbsens(ngate,gatesp,pulselmn,beamwid,sens)
wasqinv=fourln4/(efbmwid*efbmwid)
bmwrad=deg2rad*efbmwid
DO kpt=1,nptselv
DO jpt=1,nptsazm
DO ipt=1,nptsrng
wgtpt(ipt,jpt,kpt)=exp(-((delrng(ipt)*delrng(ipt))*sr2inv + &
(delazm(jpt)*delazm(jpt))*wasqinv + &
(delelv(kpt)*delelv(kpt))*wesqinv))
END DO
END DO
END DO
IF(rotdir == 1) THEN
irot=1
ELSE IF(rotdir == -1) THEN
irot=-1
ELSE IF (rotdir == 2 .OR. rotdir == -2) THEN
irot=-1*irot
END IF
IF(irot > 0) THEN
azmbgn=azim1
ELSE
azmbgn=azim2
END IF
WRITE(6,'(a,f10.2,a,i4,a,f10.2)') &
' Elev:',elv,' irot:',irot,' azmbgn:',azmbgn
IF(int_method == 1) THEN
WRITE(6,'(a,i3)') 'Using integration method ',int_method
DO jazim=1,nazim
azimuth=azmbgn+irot*((jazim-1)*delaz)
IF(azimuth > 360.) azimuth = azimuth-360.
IF(azimuth < 0.) azimuth = azimuth+360.
! IF(MOD((jazim-1),10) == 0) &
! WRITE(6,'(a,f10.2)') ' Processing azimuth: ', azimuth
azrad=deg2rad*azimuth
azim(jazim)=azimuth
beamw(jazim)=beamwid
gtspc(jazim)=gatesp
vnyq(jazim)=vnyquist(itilt)
DO igate=1, ngate
kntgate=kntgate+1
srange=igate*gatesp
CALL beamhgt(elv,srange,hgt,sfcrng)
IF(mapproj > 0 ) THEN
CALL gcircle(radlata,radlona,azimuth,sfcrng,gtlat,gtlon)
CALL lltoxy(1,1,gtlat,gtlon,gtx,gty)
delx=gtx-xrada
dely=gty-yrada
xrat=delx/srange
yrat=dely/srange
ipos=(gtx*dxinv)+1.5
jpos=(gty*dyinv)+1.5
ELSE
delx=sin(azrad)*sfcrng
dely=cos(azrad)*sfcrng
xrat=delx/srange
yrat=dely/srange
ipos=irloc+(delx*dxinv)
jpos=jrloc+(dely*dyinv)
END IF
iloc1=NINT(ipos)
jloc1=NINT(jpos)
iloc=INT(ipos)
jloc=INT(jpos)
IF(iloc1 > 0 .AND. iloc1 < nx-1 .AND. &
jloc1 > 0 .AND. jloc1 < ny-1 ) THEN
hgtmsl=hgt+radelv
IF( hgtmsl > zps(iloc1,jloc1,2) .AND. &
hgtmsl < zps(iloc1,jloc1,nz-1) ) THEN
DO k=3,nz-2
IF(zps(iloc1,jloc1,k) > hgtmsl) EXIT
END DO
bmwm=bmwrad*evalwid*srange
sradius=max(bmwm,pulsewmn)
iradius=1+INT(sradius/dx)
jradius=1+INT(sradius/dy)
! print *, 'bmwm,pulsewid = ',bmwm,pulsewmn
! print *, 'iradius,jradius = ',iradius,jradius
!
! First calculate provisional values using tri-linear interpolation
! These will be used if there is insufficient density of data for a
! weighted-average calculation
!
rmax=0.
r11l=max(refz( iloc, jloc,k-1),0.)
rmax=max(rmax,r11l)
r21l=max(refz(iloc+1, jloc,k-1),0.)
rmax=max(rmax,r21l)
r12l=max(refz( iloc,jloc+1,k-1),0.)
rmax=max(rmax,r12l)
r22l=max(refz(iloc+1,jloc+1,k-1),0.)
rmax=max(rmax,r22l)
r11h=max(refz( iloc, jloc, k),0.)
rmax=max(rmax,r11h)
r21h=max(refz(iloc+1, jloc, k),0.)
rmax=max(rmax,r21h)
r12h=max(refz( iloc,jloc+1, k),0.)
rmax=max(rmax,r12h)
r22h=max(refz(iloc+1,jloc+1, k),0.)
rmax=max(rmax,r22h)
whigh=(hgtmsl-zps(iloc1,jloc1,k-1))/ &
(zps(iloc1,jloc1,k)-zps(iloc1,jloc1,k-1))
wlow=1.-whigh
w2i=ipos-iloc
w1i=1.0-w2i
w2j=jpos-jloc
w1j=1.0-w2j
w11=w1i*w1j
w21=w2i*w1j
w12=w1i*w2j
w22=w2i*w2j
vvort(igate,jazim)=whigh*(w11*vort(iloc,jloc,k) + &
w21*vort(iloc+1,jloc,k) + &
w12*vort(iloc,jloc+1,k) + &
w22*vort(iloc+1,jloc+1,k)) + &
wlow*(w11*vort(iloc,jloc,k-1) + &
w21*vort(iloc+1,jloc,k-1) + &
w12*vort(iloc,jloc+1,k-1) + &
w22*vort(iloc+1,jloc+1,k-1))
vvmin=min(vvmin,vvort(igate,jazim))
vvmax=max(vvmax,vvort(igate,jazim))
IF( rmax > 0. ) THEN
CALL dhdrange(elv,srange,dhdr)
dsdr=sqrt(1.-dhdr*dhdr)
kntintv=kntintv+1
kntintr=kntintr+1
reflz=whigh*(w11*r11h+w21*r21h + &
w12*r12h+w22*r22h) + &
wlow*(w11*r11l+w21*r21l + &
w12*r12l+w22*r22l)
refdbz=10.*alog10(reflz)
refl(igate,jazim)=refdbz+sigmarf*randnor()
refl(igate,jazim)=max(refl(igate,jazim),refmin)
IF(refdbz > refmin ) THEN
urad=whigh*(w11*usc( iloc, jloc,k) + &
w21*usc(iloc+1, jloc,k) + &
w12*usc( iloc,jloc+1,k) + &
w22*usc(iloc+1,jloc+1,k))+ &
wlow*(w11*usc( iloc, jloc,k-1) + &
w21*usc(iloc+1, jloc,k-1) + &
w12*usc( iloc,jloc+1,k-1) + &
w22*usc(iloc+1,jloc+1,k-1))
vrad=whigh*(w11*vsc( iloc, jloc,k) + &
w21*vsc(iloc+1, jloc,k) + &
w12*vsc( iloc,jloc+1,k) + &
w22*vsc(iloc+1,jloc+1,k))+ &
wlow*(w11*vsc( iloc, jloc,k-1) + &
w21*vsc(iloc+1, jloc,k-1) + &
w12*vsc( iloc,jloc+1,k-1) + &
w22*vsc(iloc+1,jloc+1,k-1))
wrad=whigh*(w11*wsc( iloc, jloc,k) + &
w21*wsc(iloc+1, jloc,k) + &
w12*wsc( iloc,jloc+1,k) + &
w22*wsc(iloc+1,jloc+1,k))+ &
wlow*(w11*wsc( iloc, jloc,k-1) + &
w21*wsc(iloc+1, jloc,k-1) + &
w12*wsc( iloc,jloc+1,k-1) + &
w22*wsc(iloc+1,jloc+1,k-1))
!
! Estimate terminal velocity (vt is positive down)
!
CALL termvel(refl(igate,jazim),hgtmsl,vt)
!
! Calculate radial component
!
vr=dsdr*(xrat*urad+yrat*vrad)+dhdr*(wrad-vt)
!
! Add Gaussian random noise
!
vr=vr+sigmavr*randnor()
!
! Apply Doppler aliasing, if any
!
ifold=NINT(vr*twovninv)
radv(igate,jazim)=vr-ifold*twovnyq
END IF
END IF
!
! Calculate weighted average values using antenna pattern gain.
!
ibgn=max(1,(iloc1-iradius))
iend=min((nx-1),(iloc1+iradius))
jbgn=max(1,(jloc1-jradius))
jend=min((ny-1),(jloc1+jradius))
!
! First check to see of the entire observation volume area is clear.
!
echo=.false.
DO jpt=jbgn,jend
DO ipt=ibgn,iend
IF(cmpref(ipt,jpt) > 0.) echo=.true.
END DO
END DO
!
! If echo was found in 2d composite, proceed with processing 3d data.
!
IF(echo) THEN
sumwr=0.
sumwv=0.
sumv=0.
sumr=0.
kntr=0
kntv=0
DO jpt=jbgn,jend
DO ipt=ibgn,iend
DO kpt=2,nz-2
delv=elvsc(ipt,jpt,kpt)-elv
IF( abs(delv) < elradius .AND. &
refz(ipt,jpt,kpt) > 0. ) THEN
depth=0.5*(zps(ipt,jpt,kpt+1)-zps(ipt,jpt,kpt-1))
drng=rngsc(ipt,jpt,kpt)-srange
dazm=azmsc(ipt,jpt)-azimuth
IF(dazm > 180.) THEN
dazm=dazm-360.
ELSE IF( dazm < -180.) THEN
dazm=dazm+360.
END IF
wgtg=(drng*drng)*sr2inv + &
(dazm*dazm)*wasqinv + &
(delv*delv)*wesqinv
wgtr=depth*exp(-wgtg)
kntr=kntr+1
sumwr=sumwr+wgtr
sumr=sumr+wgtr*refz(ipt,jpt,kpt)
IF( ref(ipt,jpt,kpt) > refmin ) THEN
wgtv=refz(ipt,jpt,kpt)*wgtr
kntv=kntv+1
sumwv=sumwv+wgtv
sumv=sumv+wgtv*vrsc(ipt,jpt,kpt)
END IF
END IF
END DO
END DO
END DO
!
! Solve for mean reflectivity
!
! print *, ' igate,jazim,sumwr, kntr = ', &
! igate,jazim,sumwr,kntr
IF(sumwr > 0. .AND. kntr >= kntrmin ) THEN
kntavgr=kntavgr+1
reflz=sumr/sumwr
refdbz=10.*alog10(reflz)
refl(igate,jazim)=refdbz+sigmarf*randnor()
refl(igate,jazim)=max(refl(igate,jazim),refmin)
END IF
!
! Solve for mean radial velocity
!
! print *, ' igate,jazim,sumwv, kntv = ', &
! igate,jazim,sumwv,kntv
IF(sumwv > 0. .AND. kntv >= kntvmin ) THEN
kntavgv=kntavgv+1
vr=sumv/sumwv
!
! Add Gaussian random noise
!
vr=vr+sigmavr*randnor()
!
! Apply Doppler aliasing, if any
!
ifold=NINT(vr*twovninv)
radv(igate,jazim)=vr-ifold*twovnyq
END IF
ELSE ! no echo found in area
kntclr=kntclr+1
refl(igate,jazim)=0.
radv(igate,jazim)=rmisval
END IF
ELSE ! outside vertical limits of grid
kntvob=kntvob+1
END IF
ELSE ! outside grid domain
knthob=knthob+1
END IF
END DO ! igate
END DO ! jazim
!
ELSE ! Integration method 2
WRITE(6,'(a,i3)') 'Using integration method ',int_method
DO jazim=1,nazim
azimuth=azmbgn+irot*((jazim-1)*delaz)
IF(azimuth > 360.) azimuth = azimuth-360.
IF(azimuth < 0.) azimuth = azimuth+360.
! IF(MOD((jazim-1),20) == 0) &
! WRITE(6,'(a,f10.2)') ' Processing azimuth: ', azimuth
azrad=deg2rad*azimuth
azim(jazim)=azimuth
beamw(jazim)=beamwid
gtspc(jazim)=gatesp
vnyq(jazim)=vnyqstl
DO kpt=1,nptselv
DO jpt=1,nptsazm
blockfct(jpt,kpt)=1.0
END DO
END DO
DO igate=1, ngate
kntgate=kntgate+1
srange=igate*gatesp
CALL beamhgt(elv,srange,hgt,sfcrng)
IF(mapproj > 0 ) THEN
CALL gcircle(radlata,radlona,azimuth,sfcrng,gtlat,gtlon)
CALL lltoxy(1,1,gtlat,gtlon,gtx,gty)
delx=gtx-xrada
dely=gty-yrada
xrat=delx/srange
yrat=dely/srange
ipos=(gtx*dxinv)+1.5
jpos=(gty*dyinv)+1.5
ELSE
delx=sin(azrad)*sfcrng
dely=cos(azrad)*sfcrng
xrat=delx/srange
yrat=dely/srange
ipos=irloc+(delx*dxinv)
jpos=jrloc+(dely*dyinv)
END IF
iloc1=NINT(ipos)
jloc1=NINT(jpos)
iloc=INT(ipos)
jloc=INT(jpos)
IF(iloc1 > 0 .AND. iloc1 < nx-1 .AND. &
jloc1 > 0 .AND. jloc1 < ny-1 ) THEN
hgtmsl=hgt+radelv
w2i=ipos-iloc
w1i=1.0-w2i
w2j=jpos-jloc
w1j=1.0-w2j
w11=w1i*w1j
w21=w2i*w1j
w12=w1i*w2j
w22=w2i*w2j
trngt=w11*zp(iloc,jloc,2) + &
w21*zp(iloc+1,jloc,2) + &
w12*zp(iloc,jloc+1,2) + &
w22*zp(iloc+1,jloc+1,2)
topgt=w11*zps(iloc,jloc,nzm2) + &
w21*zps(iloc+1,jloc,nzm2) + &
w12*zps(iloc,jloc+1,nzm2) + &
w22*zps(iloc+1,jloc+1,nzm2)
CALL beamhgt(elvtop,srange,tophgt,sfcrngt)
CALL beamhgt(elvbot,srange,bothgt,sfcrngb)
topmsl=tophgt+radelv
botmsl=bothgt+radelv
IF(botmsl < topgt .AND. topmsl > (trngt+hblkmin)) THEN
DO k=3,nz-2
IF(zps(iloc1,jloc1,k) > hgtmsl) EXIT
END DO
!
! Calculate vorticity by doing a linear interpolation directly to
! the gate center location.
!
whigh=(hgtmsl-zps(iloc1,jloc1,k-1))/ &
(zps(iloc1,jloc1,k)-zps(iloc1,jloc1,k-1))
wlow=1.-whigh
vvort(igate,jazim)=whigh*(w11*vort(iloc,jloc,k) + &
w21*vort(iloc+1,jloc,k) + &
w12*vort(iloc,jloc+1,k) + &
w22*vort(iloc+1,jloc+1,k)) + &
wlow*(w11*vort(iloc,jloc,k-1) + &
w21*vort(iloc+1,jloc,k-1) + &
w12*vort(iloc,jloc+1,k-1) + &
w22*vort(iloc+1,jloc+1,k-1))
vvmin=min(vvmin,vvort(igate,jazim))
vvmax=max(vvmax,vvort(igate,jazim))
!
! Check to see of the entire observation volume area is clear.
!
bmwm=bmwrad*evalwid*srange
sradius=max(bmwm,pulsewmn)
iradius=1+INT(sradius/dx)
jradius=1+INT(sradius/dy)
ibgn=max(1,(iloc1-iradius))
iend=min((nx-1),(iloc1+iradius))
jbgn=max(1,(jloc1-jradius))
jend=min((ny-1),(jloc1+jradius))
echo=.false.
DO jpt=jbgn,jend
DO ipt=ibgn,iend
IF(cmpref(ipt,jpt) > 0.) THEN
echo=.true.
EXIT
END IF
END DO
END DO
IF(echo .OR. (botmsl < trngt)) THEN
!
! Calculate beam-weighted sums for an uniform array of points around this
! range-gate using bi-quadratic horizontal and linear vertical interpolation
!
kntr=0
sumr=0.
sumwr=0.
kntv=0.
sumv=0.
sumv1=0.
sumv2=0.
sumwv=0.
DO kpt=1,nptselv
elvpt=elv+delelv(kpt)
! IF(igate == 100) print *, 'kpt, hgtmsl: ',kpt,hgtmsl
DO jpt=1,nptsazm
azmpt=azimuth+delazm(jpt)
azmrpt=deg2rad*azmpt
DO ipt=1,nptsrng
srgpt=(igate*gatesp)+delrng(ipt)
IF(srgpt <= 0.0) CYCLE
CALL beamhgt(elvpt,srgpt,hgtpt,sfcrpt)
hgtmsl=hgtpt+radelv
IF(mapproj > 0 ) THEN
CALL gcircle(radlata,radlona,azmpt,sfcrpt,gtlat,gtlon)
CALL lltoxy(1,1,gtlat,gtlon,gtx,gty)
delx=gtx-xrada
dely=gty-yrada
xrat=delx/sfcrpt
yrat=dely/sfcrpt
ipos=(gtx*dxinv)+1.5
jpos=(gty*dyinv)+1.5
ELSE
delx=sin(azmrpt)*sfcrpt
dely=cos(azmrpt)*sfcrpt
xrat=delx/sfcrpt
yrat=dely/sfcrpt
ipos=irloc+(delx*dxinv)
jpos=jrloc+(dely*dyinv)
END IF
iloc1=NINT(ipos)
jloc1=NINT(jpos)
iloc=INT(ipos)
jloc=INT(jpos)
IF(iloc1 > 0 .AND. iloc1 < nx-1 .AND. &
jloc1 > 0 .AND. jloc1 < ny-1 ) THEN
ilc=max(min(iloc1,nx-2),2)
jlc=max(min(jloc1,ny-2),2)
xpt=gtx-xsc(ilc)
ypt=gty-ysc(jlc)
w2i=ipos-iloc
w1i=1.0-w2i
w2j=jpos-jloc
w1j=1.0-w2j
w11=w1i*w1j
w21=w2i*w1j
w12=w1i*w2j
w22=w2i*w2j
trnpt=w11*zp(iloc,jloc,2) + &
w21*zp(iloc+1,jloc,2) + &
w12*zp(iloc,jloc+1,2) + &
w22*zp(iloc+1,jloc+1,2)
hgtagl=hgtmsl-trnpt
IF(blockopt > 0 .AND. hgtagl < hblkmax) THEN
xmit=max(0.0,min(1.0,((hgtagl-hblkmin)*dhblkinv)))
gndatten=max(0.0,min(1.0,(blkcst*(1.-xmit))))
blockfct(jpt,kpt)= &
max(0.0,min(1.0,(blockfct(jpt,kpt)-gndatten)))
END IF
IF( hgtagl > hblkmin .AND. &
hgtmsl < zps(iloc1,jloc1,nz-1) ) THEN
DO k=3,nz-2
IF(zps(iloc1,jloc1,k) > hgtmsl) EXIT
END DO
whigh=(hgtmsl-zps(iloc1,jloc1,k-1))/ &
(zps(iloc1,jloc1,k)-zps(iloc1,jloc1,k-1))
wlow=1.-whigh
a=(refz(ilc-1,jlc-1,k-1)-2.*refz(ilc,jlc-1,k-1)+ &
refz(ilc+1,jlc-1,k-1))*twdx2inv
b=(refz(ilc+1,jlc-1,k-1)-refz(ilc-1,jlc-1,k-1)) &
*twdxinv
fjm1=a*xpt*xpt+b*xpt+refz(ilc,jlc-1,k-1)
a=(refz(ilc-1,jlc,k-1)-2.*refz(ilc,jlc,k-1)+ &
refz(ilc+1,jlc,k-1))*twdx2inv
b=(refz(ilc+1,jlc,k-1)-refz(ilc-1,jlc,k-1)) &
*twdxinv
fj0=a*xpt*xpt+b*xpt+refz(ilc,jlc,k-1)
a=(refz(ilc-1,jlc+1,k-1)-2.*refz(ilc,jlc+1,k-1)+ &
refz(ilc+1,jlc+1,k-1))*twdx2inv
b=(refz(ilc+1,jlc+1,k-1)-refz(ilc-1,jlc+1,k-1)) &
*twdxinv
fjp1=a*xpt*xpt+b*xpt+refz(ilc,jlc+1,k-1)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
flow=a*ypt*ypt + b*ypt + fj0
a=(refz(ilc-1,jlc-1,k)-2.*refz(ilc,jlc-1,k)+ &
refz(ilc+1,jlc-1,k))*twdx2inv
b=(refz(ilc+1,jlc-1,k)-refz(ilc-1,jlc-1,k)) &
*twdxinv
fjm1=a*xpt*xpt+b*xpt+refz(ilc,jlc-1,k)
a=(refz(ilc-1,jlc,k)-2.*refz(ilc,jlc,k)+ &
refz(ilc+1,jlc,k))*twdx2inv
b=(refz(ilc+1,jlc,k)-refz(ilc-1,jlc,k))*twdxinv
fj0=a*xpt*xpt+b*xpt+refz(ilc,jlc,k)
a=(refz(ilc-1,jlc+1,k)-2.*refz(ilc,jlc+1,k)+ &
refz(ilc+1,jlc+1,k))*twdx2inv
b=(refz(ilc+1,jlc+1,k)-refz(ilc-1,jlc+1,k)) &
*twdxinv
fjp1=a*xpt*xpt+b*xpt+refz(ilc,jlc+1,k)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
fhigh=a*ypt*ypt + b*ypt + fj0
reflz=whigh*fhigh+wlow*flow
kntr=kntr+1
sumr=sumr+reflz*wgtpt(ipt,jpt,kpt)*blockfct(jpt,kpt)
sumwr=sumwr+wgtpt(ipt,jpt,kpt)
IF( reflz > 0. ) THEN
CALL dhdrange(elvpt,srgpt,dhdr)
dsdr=sqrt(1.-dhdr*dhdr)
refdbz=10.*alog10(reflz)
IF(refdbz > refmin ) THEN
a=(usc(ilc-1,jlc-1,k-1)-2.*usc(ilc,jlc-1,k-1)&
+usc(ilc+1,jlc-1,k-1))*twdx2inv
b=(usc(ilc+1,jlc-1,k-1)-usc(ilc-1,jlc-1,k-1))&
*twdxinv
fjm1=a*xpt*xpt+b*xpt+usc(ilc,jlc-1,k-1)
a=(usc(ilc-1,jlc,k-1)-2.*usc(ilc,jlc,k-1)+ &
usc(ilc+1,jlc,k-1))*twdx2inv
b=(usc(ilc+1,jlc,k-1)-usc(ilc-1,jlc,k-1)) &
*twdxinv
fj0=a*xpt*xpt+b*xpt+usc(ilc,jlc,k-1)
a=(usc(ilc-1,jlc+1,k-1)-2.*usc(ilc,jlc+1,k-1)&
+usc(ilc+1,jlc+1,k-1))*twdx2inv
b=(usc(ilc+1,jlc+1,k-1)-usc(ilc-1,jlc+1,k-1))&
*twdxinv
fjp1=a*xpt*xpt+b*xpt+usc(ilc,jlc+1,k-1)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
flow=a*ypt*ypt + b*ypt + fj0
a=(usc(ilc-1,jlc-1,k)-2.*usc(ilc,jlc-1,k)+ &
usc(ilc+1,jlc-1,k))*twdx2inv
b=(usc(ilc+1,jlc-1,k)-usc(ilc-1,jlc-1,k)) &
*twdxinv
fjm1=a*xpt*xpt+b*xpt+usc(ilc,jlc-1,k)
a=(usc(ilc-1,jlc,k)-2.*usc(ilc,jlc,k)+ &
usc(ilc+1,jlc,k))*twdx2inv
b=(usc(ilc+1,jlc,k)-usc(ilc-1,jlc,k))*twdxinv
fj0=a*xpt*xpt+b*xpt+usc(ilc,jlc,k)
a=(usc(ilc-1,jlc+1,k)-2.*usc(ilc,jlc+1,k)+ &
usc(ilc+1,jlc+1,k))*twdx2inv
b=(usc(ilc+1,jlc+1,k)-usc(ilc-1,jlc+1,k)) &
*twdxinv
fjp1=a*xpt*xpt+b*xpt+usc(ilc,jlc+1,k)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
fhigh=a*ypt*ypt + b*ypt + fj0
urad=whigh*fhigh+wlow*flow
a=(vsc(ilc-1,jlc-1,k-1)-2.*vsc(ilc,jlc-1,k-1)&
+vsc(ilc+1,jlc-1,k-1))*twdx2inv
b=(vsc(ilc+1,jlc-1,k-1)-vsc(ilc-1,jlc-1,k-1))&
*twdxinv
fjm1=a*xpt*xpt+b*xpt+vsc(ilc,jlc-1,k-1)
a=(vsc(ilc-1,jlc,k-1)-2.*vsc(ilc,jlc,k-1)+ &
vsc(ilc+1,jlc,k-1))*twdx2inv
b=(vsc(ilc+1,jlc,k-1)-vsc(ilc-1,jlc,k-1)) &
*twdxinv
fj0=a*xpt*xpt+b*xpt+vsc(ilc,jlc,k-1)
a=(vsc(ilc-1,jlc+1,k-1)-2.*vsc(ilc,jlc+1,k-1)&
+vsc(ilc+1,jlc+1,k-1))*twdx2inv
b=(vsc(ilc+1,jlc+1,k-1)-vsc(ilc-1,jlc+1,k-1))&
*twdxinv
fjp1=a*xpt*xpt+b*xpt+vsc(ilc,jlc+1,k-1)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
flow=a*ypt*ypt + b*ypt + fj0
a=(vsc(ilc-1,jlc-1,k)-2.*vsc(ilc,jlc-1,k)+ &
vsc(ilc+1,jlc-1,k))*twdx2inv
b=(vsc(ilc+1,jlc-1,k)-vsc(ilc-1,jlc-1,k)) &
*twdxinv
fjm1=a*xpt*xpt+b*xpt+vsc(ilc,jlc-1,k)
a=(vsc(ilc-1,jlc,k)-2.*vsc(ilc,jlc,k)+ &
vsc(ilc+1,jlc,k))*twdx2inv
b=(vsc(ilc+1,jlc,k)-vsc(ilc-1,jlc,k))*twdxinv
fj0=a*xpt*xpt+b*xpt+vsc(ilc,jlc,k)
a=(vsc(ilc-1,jlc+1,k)-2.*vsc(ilc,jlc+1,k)+ &
vsc(ilc+1,jlc+1,k))*twdx2inv
b=(vsc(ilc+1,jlc+1,k)-vsc(ilc-1,jlc+1,k)) &
*twdxinv
fjp1=a*xpt*xpt+b*xpt+vsc(ilc,jlc+1,k)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
fhigh=a*ypt*ypt + b*ypt + fj0
vrad=whigh*fhigh+wlow*flow
a=(wsc(ilc-1,jlc-1,k-1)-2.*wsc(ilc,jlc-1,k-1)&
+wsc(ilc+1,jlc-1,k-1))*twdx2inv
b=(wsc(ilc+1,jlc-1,k-1)-wsc(ilc-1,jlc-1,k-1))&
*twdxinv
fjm1=a*xpt*xpt+b*xpt+wsc(ilc,jlc-1,k-1)
a=(wsc(ilc-1,jlc,k-1)-2.*wsc(ilc,jlc,k-1)+ &
wsc(ilc+1,jlc,k-1))*twdx2inv
b=(wsc(ilc+1,jlc,k-1)-wsc(ilc-1,jlc,k-1)) &
*twdxinv
fj0=a*xpt*xpt+b*xpt+wsc(ilc,jlc,k-1)
a=(wsc(ilc-1,jlc+1,k-1)-2.*wsc(ilc,jlc+1,k-1)&
+wsc(ilc+1,jlc+1,k-1))*twdx2inv
b=(wsc(ilc+1,jlc+1,k-1)-wsc(ilc-1,jlc+1,k-1))&
*twdxinv
fjp1=a*xpt*xpt+b*xpt+wsc(ilc,jlc+1,k-1)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
flow=a*ypt*ypt + b*ypt + fj0
a=(wsc(ilc-1,jlc-1,k)-2.*wsc(ilc,jlc-1,k)+ &
wsc(ilc+1,jlc-1,k))*twdx2inv
b=(wsc(ilc+1,jlc-1,k)-wsc(ilc-1,jlc-1,k)) &
*twdxinv
fjm1=a*xpt*xpt+b*xpt+wsc(ilc,jlc-1,k)
a=(wsc(ilc-1,jlc,k)-2.*wsc(ilc,jlc,k)+ &
wsc(ilc+1,jlc,k))*twdx2inv
b=(wsc(ilc+1,jlc,k)-wsc(ilc-1,jlc,k))*twdxinv
fj0=a*xpt*xpt+b*xpt+wsc(ilc,jlc,k)
a=(wsc(ilc-1,jlc+1,k)-2.*wsc(ilc,jlc+1,k)+ &
wsc(ilc+1,jlc+1,k))*twdx2inv
b=(wsc(ilc+1,jlc+1,k)-wsc(ilc-1,jlc+1,k)) &
*twdxinv
fjp1=a*xpt*xpt+b*xpt+wsc(ilc,jlc+1,k)
a=(fjm1-2.*fj0+fjp1)*twdy2inv
b=(fjp1-fjm1)*twdyinv
fhigh=a*ypt*ypt + b*ypt + fj0
wrad=whigh*fhigh+wlow*flow
!
! Estimate terminal velocity (vt is positive down)
!
CALL termvel(refdbz,hgtmsl,vt)
!
! Calculate radial component
!
vr=dsdr*(xrat*urad+yrat*vrad)+dhdr*(wrad-vt)
kntv=kntv+1
sumv=sumv+vr*wgtpt(ipt,jpt,kpt)* &
blockfct(jpt,kpt)*reflz
sumwv=sumwv+wgtpt(ipt,jpt,kpt)*reflz
sumv1=sumv1+vr
sumv2=sumv2+vr*vr
END IF ! refdbz > refmin
END IF ! reflz > 0
ELSE IF(hgtagl <= hblkmin .AND. blockopt > 1) THEN
vr=0.
reflz=reflzgnd
kntr=kntr+1
sumr=sumr+reflz*wgtpt(ipt,jpt,kpt)*blockfct(jpt,kpt)
sumwr=sumwr+wgtpt(ipt,jpt,kpt)
kntv=kntv+1
sumv=sumv+vr*wgtpt(ipt,jpt,kpt)* &
blockfct(jpt,kpt)*reflz
sumwv=sumwv+wgtpt(ipt,jpt,kpt)*reflz
sumv1=sumv1+vr
sumv2=sumv2+vr*vr
END IF ! hgt in grid
END IF ! iloc in grid
END DO ! ipt
END DO ! jpt
END DO ! kpt
!
! Calculate reflectivity and velocity from the uniform
! array center on this range gate.
!
IF(sumwr > 0. .AND. kntr > kntrmin) THEN
kntintr=kntintr+1
reflz=sumr/sumwr
refdbz=10.*alog10(reflz)
refl(igate,jazim)=min(max(refdbz,refmin),refmax)
ELSE
kntfewr=kntfewr+1
refl(igate,jazim)=0.
END IF
IF(sumwv > 0. .AND. kntv > kntvmin) THEN
kntintv=kntintv+1
vr=sumv/sumwv
ifold=NINT(vr*twovninv)
radv(igate,jazim)=vr-ifold*twovnyq
flkntv=float(kntv)
varvr=max(0.01, &
((sumv2-(sumv1*sumv1/flkntv))/(flkntv-1.)))
stdvr(igate,jazim)=stdvrmul*SQRT(varvr)
ELSE
kntfewv=kntfewv+1
radv(igate,jazim)=rmisval
END IF
ELSE ! no echo
kntclr=kntclr+1
refl(igate,jazim)=0.
radv(igate,jazim)=rmisval
END IF
ELSE IF ( botmsl > topgt) THEN ! beam completely above grid
kntvob=kntvob+(ngate-igate)+1
kntgate=kntgate+(ngate-igate)
EXIT
ELSE IF ( topmsl < (trnpt+hblkmin) ) THEN ! beam completely below ground
print *, ' Entire beam below ground: ',azimuth,igate
IF(blockopt > 1) THEN
sumr=0.
sumwr=0.
DO kpt=1,nptselv
DO jpt=1,nptsazm
reflz=reflzgnd
sumr=sumr+reflz*wgtpt(ipt,jpt,kpt)*blockfct(jpt,kpt)
sumwr=sumwr+wgtpt(ipt,jpt,kpt)
END DO
END DO
refl(igate,jazim)=sumr/sumwr
ELSE
refl(igate,jazim)=0.
END IF
radv(igate,jazim)=0.
kntbob=kntbob+(ngate-igate)+1
kntgate=kntgate+(ngate-igate)
DO iigate=igate+1,ngate
radv(iigate,jazim)=0.
refl(iigate,jazim)=0.
END DO
EXIT
END IF
ELSE ! outside grid domain
knthob=knthob+1
END IF
END DO ! igate
reflmin=999.
reflmax=-999.
DO igate=1,ngate
IF(refl(igate,jazim) > refchk) THEN
reflmin=min(reflmin,refl(igate,jazim))
reflmax=max(reflmax,refl(igate,jazim))
END IF
END DO
IF(mod(jazim,30) == 0) &
print *, ' Azimuth,reflmin,reflmax=',azimuth,reflmin,reflmax
END DO ! jazim
END IF ! integration method
!
!-----------------------------------------------------------------------
!
! Calculate attenuation for this tilt
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a,i6)') ' Attenopt block, attenopt= ',attenopt
IF (attenopt > 0) THEN
CALL xbatten(ngate,maxazim,nazim, &
refl,gatesp,rmisval,refmin,attrefl)
ELSE
DO jazim=1,nazim
DO igate=1,ngate
attrefl(igate,jazim)=refl(igate,jazim)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Apply noise to the reflectivity and velocity according the
! indicated option.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a,i4)') ' radial vel error, vrerropt= ',vrerropt
IF(vrerropt > 0 ) &
CALL radverr(ngate,maxazim,nazim,radv,attrefl,stdvr,gatesp, &
vrerropt,sigmavr,samploi,vnyqstl, &
wavelen,pwrxmt,gaindb,lossdb,noisedbm,beamwid, &
pulselmn,npulsmn,rmisval)
WRITE(6,'(a,i4)') ' reflectivity error, rferropt= ',rferropt
IF(rferropt > 0 ) &
CALL reflerr(ngate,maxazim,nazim,refl,attrefl,gatesp, &
rferropt,sigmarf,refmin,samploi, &
wavelen,pwrxmt,gaindb,lossdb,noisedbm,beamwid, &
pulselmn,npulsmn,rmisval)
!
!-----------------------------------------------------------------------
!
! Apply X-band radar sensitivity thresholding to the
! attenuated velocity using the attenuated reflectivity.
!
!-----------------------------------------------------------------------
!
IF(senstvopt > 0) THEN
DO jazim=1,nazim
DO igate=1,ngate
IF(attrefl(igate,jazim) > sens(igate)) THEN
attradv(igate,jazim)=radv(igate,jazim)
ELSE
attrefl(igate,jazim)=rmisval
attradv(igate,jazim)=rmisval
END IF
END DO
END DO
ELSE
DO jazim=1,nazim
DO igate=1,ngate
attradv(igate,jazim)=radv(igate,jazim)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Report statistics
!
!-----------------------------------------------------------------------
!
denom=1./float(kntgate)
write(6,'(a,i9)') ' Total gates: ',kntgate
write(6,'(a,i9,f7.1,a)') ' Out-of-bounds Horiz: ',knthob, &
(100.*float(knthob)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Above grid points: ',kntvob, &
(100.*float(kntvob)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Ground-blocked points: ',kntbob, &
(100.*float(kntbob)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Gates no echo: ',kntclr, &
(100.*float(kntclr)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Gates interp refl: ',kntintr, &
(100.*float(kntintr)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Gates interp rvel: ',kntintv, &
(100.*float(kntintv)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Gates too few pts refl: ',kntfewr, &
(100.*float(kntfewr)*denom),' percent'
write(6,'(a,i9,f7.1,a)') ' Gates too few pts rvel: ',kntfewv, &
(100.*float(kntfewv)*denom),' percent'
write(6,'(a,f16.8,a,f16.8)') &
' Tilt min vert vort (*10**3)=',vvmin,' Max vert vort = ',vvmax
print *, ' Stats by azimuth, refl:'
DO jazim=1,nazim,30
reflmin=999.
reflmax=-999.
DO igate=1,ngate
IF(refl(igate,jazim) > refchk) THEN
reflmin=min(reflmin,refl(igate,jazim))
reflmax=max(reflmax,refl(igate,jazim))
END IF
END DO
print *, ' Azimuth,reflmin,reflmax=',azim(jazim),reflmin,reflmax
END DO
print *, ' Stats by azimuth, Attenuted refl:'
DO jazim=1,nazim,30
reflmin=999.
reflmax=-999.
DO igate=1,ngate
reflmin=min(reflmin,attrefl(igate,jazim))
reflmax=max(reflmax,attrefl(igate,jazim))
END DO
print *, ' Azimuth,reflmin,reflmax=',azim(jazim),reflmin,reflmax
END DO
CALL abss2ctim(itimtilt,iyear,imon,iday,ihour,imin,isec)
itimcdf=itimtilt-itim1970
!
!-----------------------------------------------------------------------
!
! Write tilt netCDF file - Reflectivity
!
!-----------------------------------------------------------------------
!
IF( ifmt == 1) THEN
itimvol(itilt)=itimtilt
vnyqvol(itilt)=vnyqstl
DO igate=1,ngate
rngvol(igate,itilt)=igate*gatesp
END DO
DO jazim=1,nazim
elvvol(jazim,itilt)=elvang(itilt)
azmvol(jazim,itilt)=azim(jazim)
END DO
DO jazim=1,nazim
DO igate=1,ngate
refvol(igate,jazim,itilt)=attrefl(igate,jazim)
velvol(igate,jazim,itilt)=attradv(igate,jazim)
uarefvol(igate,jazim,itilt)=refl(igate,jazim)
uavelvol(igate,jazim,itilt)=radv(igate,jazim)
vorvol(igate,jazim,itilt)=vvort(igate,jazim)
END DO
END DO
ELSE IF( ifmt == 2 ) THEN
write(fname,'(a,i2.2,a,i2.2,a,i4.4,2(i2.2),a,3(i2.2),a)') &
'Reflectivity_',ielv,'.',felv,'_',iyear,imon,iday,'-', &
ihour,imin,isec,'.netcdf'
varname='Reflectivity'
CALL wtrftiltcdf(ngate,maxazim,nazim,fname,outdir,varname, &
radname,radlata,radlona,radelv,vcp,elv, &
rmisval,rngfval,itimtilt,frtime,initime, &
vnyqstl,rfrgate, &
azim,beamw,gtspc,attrefl)
!
!-----------------------------------------------------------------------
!
! Write index record
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 ) THEN
write(idxunit,'(a)') '<item>'
write(idxunit,'(a,f8.6,a,i10,a)') &
'<time fractional="',frtime,'">',itimcdf,'</time>'
write(idxunit,'(a,a,a,a)') '<params>netcdf {indexlocation} ', &
TRIM(fname),TRIM(cmprext),'</params>'
write(idxunit,'(a,i4.4,2i2.2,a,3i2.2,1x,a,1x,i2.2,a1,i2.2,a)')&
'<selections>',iyear,imon,iday,'-', &
ihour,imin,isec,'Reflectivity',ielv,'.',felv, &
'</selections>'
write(idxunit,'(a)') '</item>'
END IF ! creidx
!
!-----------------------------------------------------------------------
!
! Write tilt netCDF file - Velocity
!
!-----------------------------------------------------------------------
!
write(fname,'(a,i2.2,a,i2.2,a,i4.4,2(i2.2),a,3(i2.2),a)') &
'Velocity_',ielv,'.',felv,'_',iyear,imon,iday,'-', &
ihour,imin,isec,'.netcdf'
varname='RadialVelocity'
CALL wtvrtiltcdf(ngate,maxazim,nazim,fname,outdir,varname, &
radname,radlat,radlon,radelv,vcp,elv, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyqstl,rfrgate, &
azim,beamw,gtspc,vnyq,attradv)
!
!-----------------------------------------------------------------------
!
! Write index record
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 ) THEN
write(idxunit,'(a)') '<item>'
write(idxunit,'(a,f8.6,a,i10,a)') &
'<time fractional="',frtime,'">',itimcdf,'</time>'
write(idxunit,'(a,i4.4,2i2.2,a,3i2.2,1x,a,1x,i2.2,a1,i2.2,a)')&
'<selections>',iyear,imon,iday,'-', &
ihour,imin,isec,'Velocity',ielv,'.',felv, &
'</selections>'
write(idxunit,'(a)') '</item>'
END IF ! creidx
!
write(fname,'(a,i2.2,a,i2.2,a,i4.4,2(i2.2),a,3(i2.2),a)') &
'ReflectUnatten_',ielv,'.',felv,'_',iyear,imon,iday,'-', &
ihour,imin,isec,'.netcdf'
varname='UnattenuatedReflectivity'
CALL wtrftiltcdf(ngate,maxazim,nazim,fname,outdir,varname, &
radname,radlata,radlona,radelv,vcp,elv, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyqstl,rfrgate, &
azim,beamw,gtspc,refl)
!
!-----------------------------------------------------------------------
!
! Write index record
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 ) THEN
write(idxunit,'(a)') '<item>'
write(idxunit,'(a,f8.6,a,i10,a)') &
'<time fractional="',frtime,'">',itimcdf,'</time>'
write(idxunit,'(a,a,a,a)') '<params>netcdf {indexlocation} ', &
TRIM(fname),TRIM(cmprext),'</params>'
write(idxunit,'(a,i4.4,2i2.2,a,3i2.2,1x,a,1x,i2.2,a1,i2.2,a)')&
'<selections>',iyear,imon,iday,'-', &
ihour,imin,isec,'Reflectivity',ielv,'.',felv, &
'</selections>'
write(idxunit,'(a)') '</item>'
END IF ! creidx
!
!-----------------------------------------------------------------------
!
! Write tilt netCDF file - Velocity
!
!-----------------------------------------------------------------------
!
write(fname,'(a,i2.2,a,i2.2,a,i4.4,2(i2.2),a,3(i2.2),a)') &
'VelocityUnatten_',ielv,'.',felv,'_',iyear,imon,iday,'-', &
ihour,imin,isec,'.netcdf'
varname='UnattenuatedVelocity'
CALL wtvrtiltcdf(ngate,maxazim,nazim,fname,outdir,varname, &
radname,radlat,radlon,radelv,vcp,elv, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyqstl,rfrgate, &
azim,beamw,gtspc,vnyq,radv)
!
!-----------------------------------------------------------------------
!
! Write index record
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 ) THEN
write(idxunit,'(a)') '<item>'
write(idxunit,'(a,f8.6,a,i10,a)') &
'<time fractional="',frtime,'">',itimcdf,'</time>'
write(idxunit,'(a,i4.4,2i2.2,a,3i2.2,1x,a,1x,i2.2,a1,i2.2,a)')&
'<selections>',iyear,imon,iday,'-', &
ihour,imin,isec,'Velocity',ielv,'.',felv, &
'</selections>'
write(idxunit,'(a)') '</item>'
END IF ! creidx
!-----------------------------------------------------------------------
!
! Write tilt netCDF file - Vorticity
!
!-----------------------------------------------------------------------
!
write(fname,'(a,i2.2,a,i2.2,a,i4.4,2(i2.2),a,3(i2.2),a)') &
'Vorticity_',ielv,'.',felv,'_',iyear,imon,iday,'-', &
ihour,imin,isec,'.netcdf'
varname='VerticalVorticity'
CALL wtvvtiltcdf(ngate,maxazim,nazim,fname,outdir,varname, &
radname,radlat,radlon,radelv,vcp,elv, &
rmisval,rngfval,itimcdf,frtime,initime, &
vnyqstl,rfrgate, &
azim,beamw,gtspc,vvort)
!
!-----------------------------------------------------------------------
!
! Write index record
!
!-----------------------------------------------------------------------
!
IF( creidx > 0 ) THEN
write(idxunit,'(a)') '<item>'
write(idxunit,'(a,f8.6,a,i10,a)') &
'<time fractional="',frtime,'">',itimcdf,'</time>'
write(idxunit,'(a,i4.4,2i2.2,a,3i2.2,1x,a,1x,i2.2,a1,i2.2,a)')&
'<selections>',iyear,imon,iday,'-', &
ihour,imin,isec,'Vorticity',ielv,'.',felv, &
'</selections>'
write(idxunit,'(a)') '</item>'
END IF ! creidx
END IF ! ifmt=2
IF(tmadvopt == 2) itimtilt=itimtilt+NINT(azmsect/rotrate)
!
END DO ! itilt
IF( ifmt == 1 ) THEN
write(fnamevol,'(4a,i4.4,i2.2,i2.2,a,i2.2,i2.2,i2.2,a)') &
TRIM(outdir),'/',TRIM(radname),'_',year,month,day,'_', &
hour,minute,second,'.vol'
CALL wtradvol(fnamevol,ngate,ngate,maxazim,nelev, &
rngvol,azmvol,elvvol,velvol,uavelvol, &
vorvol,vnyqvol,itimvol, &
rngvol,azmvol,elvvol,refvol,uarefvol, &
wrtuaref,wrtuavel,wrtvort, &
runname,radname,vcp,curtim,beamwid, &
rmisval,rngfval,radelv,radlat,radlon)
END IF
!
!-----------------------------------------------------------------------
!
! Release memory for gridded data
!
!-----------------------------------------------------------------------
WRITE(6,'(a)') ' Deallocating grid memory'
!
deallocate(usc)
deallocate(vsc)
deallocate(wsc)
deallocate(latsc)
deallocate(lonsc)
deallocate(azmsc)
deallocate(sfcr)
deallocate(elvsc)
deallocate(rngsc)
deallocate(vort)
deallocate(vrsc)
deallocate(ref)
deallocate(refz)
deallocate(cmpref)
!
! Release memory for polar coordinate radar data
!
WRITE(6,'(a)') ' Deallocating radar tilt memory'
deallocate(azim)
deallocate(beamw)
deallocate(gtspc)
deallocate(vnyq)
deallocate(radv)
deallocate(attradv)
deallocate(refl)
deallocate(attrefl)
deallocate(stdvr)
deallocate(vvort)
deallocate(sens)
WRITE(6,'(a)') ' Deallocating radar volume memory'
IF( ifmt == 1 ) THEN
deallocate(itimvol)
deallocate(vnyqvol)
deallocate(rngvol)
deallocate(azmvol)
deallocate(elvvol)
deallocate(refvol)
deallocate(uarefvol)
deallocate(velvol)
deallocate(uavelvol)
deallocate(vorvol)
END IF
WRITE(6,'(a)') ' End of file loop'
!
END DO ! ifile loop
!
! Finish index file and close it out
!
IF( creidx > 0 .AND. ifmt == 2 ) THEN
write(idxunit,'(a)') '</codeindex>'
CLOSE(idxunit)
CALL retunit(idxunit)
END IF
STOP
END PROGRAM rademul
SUBROUTINE termvel(refl,hgtmsl,vt) 3
!
!-----------------------------------------------------------------------
!
! Calculates the precipitation terminal velocity from
! the observed radial velocity. From Zeigler, 1978.
!
! An empirical formula based on the reflectivity is
! used to estimate the terminal velocity.
!
! For simplicity, atmospheric density is approxmated
! by a standard atmosphere, rho(z)=rhonot exp(-z/h0)
!
!-----------------------------------------------------------------------
!
!
IMPLICIT NONE
REAL :: refl
REAL :: hgtmsl
REAL :: vt
!
REAL, PARAMETER :: zfrez=3000.
REAL, PARAMETER :: zice=8000.
REAL, PARAMETER :: rho0=1.2250
REAL, PARAMETER :: h0=7000.
REAL, PARAMETER :: denom=(1./(zice-zfrez))
REAL, PARAMETER :: dbzmin=10.
REAL, PARAMETER :: dbzmax=100.
!
REAL :: refz,rhofact,s1,s2
!
refz=10.**(0.1*refl)
rhofact=EXP(0.4*hgtmsl/h0)
IF(hgtmsl < zfrez) THEN
vt=2.6*(refz**0.107)*rhofact
ELSE IF(hgtmsl < zice) THEN
s1=(zice-hgtmsl)*denom
s2=2.*(hgtmsl-zfrez)*denom
vt=s1*2.6*(refz**0.107)*rhofact + s2
ELSE
vt=2.0
END IF
RETURN
END SUBROUTINE termvel
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WTRADVOL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wtradvol(fnamevol,maxgatevel,maxgateref,maxazim,maxelev, & 1,1
rngvvol,azmvvol,elvvvol,velvol,uavelvol, &
vorvol,vnyqvol,itimvol, &
rngrvol,azmrvol,elvrvol,refvol,uarefvol, &
wrtuaref,wrtuavel,wrtvort, &
runname,radname,vcp,curtim,beamwid, &
rmisval,rngfval,radar_alt,radar_lat,radar_lon)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write out radar observations in radar coordinate
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Hu, CAPS
!
! MODIFICATION HISTORY:
! Converted to subroutine WTRADVOL, Keith Brewster
!
! 11/13/2004 Keith Brewster, CAPS
! Updated to write more variables for use by radbin2cdf
!
! 05/24/2005 Keith Brewster, CAPS
! Added unattenuated variables and variable write switches
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! maxgate Maximum gates in a radial
! maxazim Maximum radials per tilt
! maxelev Maximum number of tilts
!
IMPLICIT NONE
!
CHARACTER (LEN=256), INTENT(IN) :: fnamevol
INTEGER, INTENT(IN) :: maxgatevel
INTEGER, INTENT(IN) :: maxgateref
INTEGER, INTENT(IN) :: maxazim
INTEGER, INTENT(IN) :: maxelev
INTEGER, INTENT(IN) :: itimvol(maxelev)
REAL, INTENT(IN) :: rngvvol(maxgatevel,maxelev)
REAL, INTENT(IN) :: azmvvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvvvol(maxazim,maxelev)
REAL, INTENT(IN) :: velvol(maxgatevel,maxazim,maxelev)
REAL, INTENT(IN) :: uavelvol(maxgatevel,maxazim,maxelev)
REAL, INTENT(IN) :: vorvol(maxgatevel,maxazim,maxelev)
REAL, INTENT(IN) :: vnyqvol(maxelev)
REAL, INTENT(IN) :: rngrvol(maxgateref,maxelev)
REAL, INTENT(IN) :: azmrvol(maxazim,maxelev)
REAL, INTENT(IN) :: elvrvol(maxazim,maxelev)
REAL, INTENT(IN) :: refvol(maxgateref,maxazim,maxelev)
REAL, INTENT(IN) :: uarefvol(maxgateref,maxazim,maxelev)
INTEGER, INTENT(IN) :: wrtuaref,wrtuavel,wrtvort
CHARACTER (LEN=80 ), INTENT(IN) :: runname
CHARACTER (LEN=4 ), INTENT(IN) :: radname
REAL, INTENT(IN) :: beamwid,rmisval,rngfval
INTEGER, INTENT(IN) :: vcp
REAL, INTENT(IN) :: curtim
REAL, INTENT(IN) :: radar_alt
REAL, INTENT(IN) :: radar_lat
REAL, INTENT(IN) :: radar_lon
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iyear,imon,iday,ihour,imin,isec,idummy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
idummy=0
CALL abss2ctim(itimvol(1),iyear,imon,iday,ihour,imin,isec)
!
!-----------------------------------------------------------------------
!
! Open file for output
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a,a)') 'Opening: ',TRIM(fnamevol)
OPEN(13,FILE=trim(fnamevol),STATUS='unknown', &
FORM='unformatted')
WRITE(13) runname
WRITE(13) radname
WRITE(13) maxgatevel,maxgateref,maxazim,maxelev
WRITE(13) iyear,imon,iday,ihour,imin,isec
WRITE(13) radar_alt,radar_lat,radar_lon
WRITE(13) vcp,beamwid,rmisval,rngfval,curtim
WRITE(13) wrtuaref,wrtuavel,wrtvort, &
idummy,idummy,idummy,idummy,idummy
WRITE(13) itimvol
WRITE(13) rngrvol
WRITE(13) azmrvol
WRITE(13) elvrvol
WRITE(13) refvol
IF(wrtuaref /= 0) WRITE(13) uarefvol
WRITE(13) vnyqvol
WRITE(13) rngvvol
WRITE(13) azmvvol
WRITE(13) elvvvol
WRITE(13) velvol
IF(wrtuavel /= 0) WRITE(13) uavelvol
IF(wrtvort /= 0) WRITE(13) vorvol
close(13)
WRITE(6,'(a,a)') ' Writing complete for file ',TRIM(fnamevol)
RETURN
END SUBROUTINE wtradvol
!
SUBROUTINE xbatten(ngate,maxazim,nazim, & 2
refl,gatesp,rmisval,refmin,attrefl)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate X-band (3.2 cm) attenuated reflectivity given unattenuated
! reflectivty (from S-band data or from simulation model hydrometeors)
! for a single tilt of data.
!
! Uses attenuation coefficients from Burrows and Atwood (1949) data as
! reported in Doviak and Zrnic', Doppler Radar and Weather Observations,
! 2nd Ed., Eq. 3.16c.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS and Erin Fay, OU/SoM
! November, 2004
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! ngate: Number of gates in each radial
! nazim: Number of radials of data
! refl: Unattenuated reflectivity factor dBZ
! gatesp: Gate spacing (m)
! rmisval: Missing data value, assumed below threshold
! refmin: Threshold minimum reflectivity (dBZ)
!
!
! OUTPUT:
! attrefl: Attenuated reflectivity.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: ngate ! Number of gates
INTEGER, INTENT (IN) :: maxazim ! Maximum number of radials
INTEGER, INTENT (IN) :: nazim ! Number of radials
REAL, INTENT (IN):: refl(ngate,maxazim) ! Reflectivity (dBZ)
REAL, INTENT (IN) :: gatesp ! Gate spacing (m)
REAL, INTENT (IN) :: rmisval
REAL, INTENT (IN) :: refmin
REAL, INTENT (OUT) :: attrefl(ngate,maxazim) ! Attenuated reflectivities
INTEGER :: igate,jazim ! Index.
REAL :: att_tally ! Running tally of attenuation values.
REAL :: att_const1 ! Attenuation constant 1.
REAL :: att_const2 ! Attenuation constant 2.
REAL :: att_exp ! Attenuation exponent
REAL :: atten,atten_last,hlfgatsp
REAL, PARAMETER :: zrconst=400.
REAL, PARAMETER :: zrexp=1.4
REAL, PARAMETER :: krconst=0.01
REAL, PARAMETER :: krexp=1.21
REAL, PARAMETER :: mperkm=0.001
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
att_const1 = 2.0*krconst*mperkm ! leading constant (two-way)
att_const2 = 1.0/zrconst
att_exp = krexp/zrexp
hlfgatsp = 0.5*gatesp
DO jazim=1, nazim
IF(refl(1,jazim) .gt. rmisval) THEN
atten = &
att_const1*((att_const2*(10.0**(0.1*refl(1,jazim))))**att_exp)
att_tally = atten*hlfgatsp
attrefl(1,jazim) = refl(1,jazim)-att_tally
attrefl(1,jazim) = max(attrefl(1,jazim),refmin)
ELSE
atten=0.
att_tally=0.
attrefl(1,jazim)=refl(1,jazim)
END IF
atten_last=atten
DO igate = 2, ngate
IF(refl(igate,jazim) .gt. rmisval) THEN
atten = &
att_const1*((att_const2*(10.0**(0.1*refl(igate,jazim))))**att_exp)
att_tally = att_tally + hlfgatsp*atten_last + hlfgatsp*atten
attrefl(igate,jazim) = refl(igate,jazim)-att_tally
attrefl(igate,jazim) = max(attrefl(igate,jazim),refmin)
atten_last=atten
ELSE
attrefl(igate,jazim) = refl(igate,jazim)
attrefl(igate,jazim) = max(attrefl(igate,jazim),refmin)
atten_last=0.
END IF
END DO
END DO
RETURN
END SUBROUTINE xbatten
SUBROUTINE xbsens(ngate,gatesp,pulselen,beamwid,sens) 2
!
!-----------------------------------------------------------------------
! Subroutine xbsens
! Calculate X-band radar sensitivity (dBZ) as a function of range.
!
! Based on information provided by Francesc Junyent
!
! Keith Brewster and Erin Fay, CAPS
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: ngate
REAL, INTENT(IN) :: gatesp
REAL, INTENT(IN) :: pulselen
REAL, INTENT(IN) :: beamwid
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 :: pt = 12500 ! Peak power (Watts)
REAL, PARAMETER :: G = 38.0 ! Antenna gain (dB)
REAL, PARAMETER :: F = 5.5 ! Noise figure (dB)
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
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,ntilt
REAL :: pulsel,pulsew
REAL :: ln2,pi,pi3,four3,bwrad,rnoise,BHz,Ni,sconst,rlm,rG
REAL :: range
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ln2=log(2.0)
pi=acos(-1.)
pi3=pi*pi*pi
four3=4.*4.*4.
bwrad=beamwid*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*pulselen)) * ((lmbda*lmbda*four3*rlm)/(Pt*rG*rG))
print *, ' pulselen=',pulselen
print *, ' beamwid=',beamwid
print *, ' Sensitivity................'
DO igate = 1, ngate
range=igate*gatesp
sens(igate) = 10.*LOG10(range*range*sconst)
IF(mod(igate,10) == 0) print *, (0.001*range),sens(igate)
END DO
END SUBROUTINE xbsens
FUNCTION EFFBMW(beamwid,rotrate,prt,npulse)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Estimate effective half-power beamwidth given the antenna half-power
! beamwidth, the rotation rate, pulse repetition time and the number of
! pulses used in estimating the moment data.
!
! Based on solutions to Eq. 7.34 in Doviak and Zrnic',
! Doppler Radar and Weather Observations, 2nd Ed. (1993)
!
! AUTHOR: Keith Brewster, CAPS
! November, 2004
!
! INPUT:
! beamwid Half-power antenna beamwidth (degrees)
! rotrate Rotation rate (degrees per second)
! prt Pulse repetition time (seconds)
! npulse Number of pulses averaged to get moment data
!
! OUTPUT:
! effbmw Effective half-power beamwidth (degrees)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL :: effbmw
!
REAL :: beamwid
REAL :: rotrate
REAL :: prt
INTEGER :: npulse
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: ntry=30
REAL,PARAMETER :: eps = 1.0E-06
!
REAL :: sqrtln4,twortln4,bmwinv,amts,erfcst
REAL :: gtheta,theta,test,test0,test2,delth
REAL :: theta1,theta2
INTEGER :: itry
!
!-----------------------------------------------------------------------
!
! External functions
!
!-----------------------------------------------------------------------
!
REAL :: erf
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
sqrtln4=sqrt(alog(4.))
twortln4=2.0*sqrtln4
bmwinv=1.0/beamwid
amts=abs(rotrate)*npulse*prt
erfcst=0.5*erf(sqrtln4*amts*bmwinv)
!
!-----------------------------------------------------------------------
!
! Iteratively solve for zeros. First for positive solution angle.
!
!-----------------------------------------------------------------------
!
gtheta=beamwid+amts
theta=0.
test0=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
delth=gtheta
theta=gtheta
test2=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
DO itry=1,ntry
delth=0.5*sign(delth,(test0*test2))
theta=theta+delth
test2=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
IF(abs(test2) < eps) EXIT
END DO
theta2=theta
WRITE(6,'(a,f10.4,a,g12.4,a,i3)') &
' EFFBMW Positive sol, theta=',theta2,' diffr=',test2, &
' iter=',itry
!
!-----------------------------------------------------------------------
!
! Solve for negative solution angle.
!
!-----------------------------------------------------------------------
!
theta=-gtheta
test0=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
delth=gtheta
theta=0.
test2=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
DO itry=1,ntry
delth=0.5*sign(delth,(test0*test2))
theta=theta+delth
test2=erf(twortln4*theta*bmwinv)-erf(twortln4*(theta-amts)*bmwinv)- &
erfcst
IF(abs(test2) < eps) EXIT
END DO
theta1=theta
WRITE(6,'(a,f10.4,a,g12.4,a,i3)') &
' EFFBMW Negative sol, theta=',theta1,' diffr=',test2, &
' iter=',itry
!
!-----------------------------------------------------------------------
!
! Effective beamwidth is difference between positive and negative
! solution angles.
!
!-----------------------------------------------------------------------
!
effbmw=(theta2-theta1)
WRITE(6,'(a,f12.4,a)') &
' EFFBMW Effective half-power beamwidth =',effbmw,' degrees'
RETURN
END FUNCTION EFFBMW
!
FUNCTION RANDNOR()
!
! Returns a random number with a normal (Gaussian) distribution
! with a mean of zero and standard deviation of one.
! A distribution with a mean of m and standard deviation of s would
! use the result of this function in the following way:
! x = m + s*randnor()
!
! The algorithm used here to convert a uniform random distribution (0,1)
! to a Gaussian is the polar form of the Box-Muller transformation. This
! code based on sample C code and a description of the technique given
! by Dr. Everett Carter of Taygeta Scientific on the following website:
! http://www.taygeta.com/random/gaussian.html
!
! I've added a test to avoid dividing by zero, rsq > epsilon, which doesn't
! seem to affect the statistics of the distribution.
!
! random_number is a fortran-90 intrinsic function. Implementation on
! different machines and compilers may vary.
!
! Keith Brewster, CAPS
! August, 2004
!
! MODIFICATIONS
!
! 02/28/2005 Keith Brewster
! Replaced intrinsic function rand with random_number, that seems to be
! available on more machines than "rand".
!
!
IMPLICIT NONE
REAL, PARAMETER :: epsilon = 1.0E-12
REAL :: randnor
REAL :: xr,yr,rsq,r
REAL :: myranf
DO
xr = 2.*myranf() - 1.
yr = 2.*myranf() - 1.
rsq = xr*xr + yr*yr
IF(rsq < 1.0 .AND. rsq > epsilon) EXIT
END DO
r = sqrt( (-2.0 * log( rsq ) ) / rsq )
randnor = xr * r
RETURN
END FUNCTION RANDNOR
REAL FUNCTION myranf ()
!
! A simple random number generator provided by Vijay Lakamraju UMASS, ECE
! as ranf is not implemented on all compilers.
!
! Generates a random number between 0 and 1.
!
INTEGER L, C, M
PARAMETER ( L = 1029, C = 221591, M = 1048576 )
INTEGER SEED
SAVE SEED
DATA SEED / 0 /
SEED = MOD ( SEED * L + C, M )
myranf = REAL ( SEED ) / M
RETURN
END FUNCTION myranf
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION ERF ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION erf(x)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ERF computes the error function.
!
!
! Definition:
!
! ERF(X) = ( 2 / SQRT ( PI ) ) * Integral ( 0 <= T <= X ) EXP ( -T**2 ) dT
!
!
! Reference:
!
! David Kahaner, Clever Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1988.
!
! INPUT:
!
! real x, the argument of the error function.
!
! OUTPUT:
! real erf, the value of the error function at X.
!
!
! AUTHOR: John Burkardt, Pittsburgh Supercomputing Center
! NMS Fortran Software Package
!
! MODIFICATION HISTORY:
! 23-Oct-2004 Modified for ARPS Standard Fortran-90
! Keith Brewster, CAPS Univ of Oklahoma
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL :: x
REAL :: erf
!
REAL :: csevl
REAL :: erfc
!
REAL :: erfcs(13)
INTEGER inits
REAL :: r1mach
REAL :: sqeps
REAL :: sqrtpi
REAL :: y
!
REAL, SAVE :: xbig = 0.
INTEGER, SAVE :: nterf = 0
!
data erfcs( 1) / -0.049046121234691808 /
data erfcs( 2) / -0.14226120510371364 /
data erfcs( 3) / 0.010035582187599796 /
data erfcs( 4) / -0.000576876469976748 /
data erfcs( 5) / 0.000027419931252196 /
data erfcs( 6) / -0.000001104317550734 /
data erfcs( 7) / 0.000000038488755420 /
data erfcs( 8) / -0.000000001180858253 /
data erfcs( 9) / 0.000000000032334215 /
data erfcs(10) / -0.000000000000799101 /
data erfcs(11) / 0.000000000000017990 /
data erfcs(12) / -0.000000000000000371 /
data erfcs(13) / 0.000000000000000007 /
data sqeps / 0. /
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
sqrtpi=sqrt(acos(-1.))
!
IF ( nterf == 0 ) THEN
nterf = inits ( erfcs, 13, 0.1*r1mach(3) )
xbig = sqrt ( - log ( sqrtpi * r1mach(3) ) )
sqeps = sqrt ( 2. * epsilon ( sqeps ) )
END IF
y = abs ( x )
IF ( y <= sqeps ) THEN
erf = 2.0 * x / sqrtpi
ELSE IF ( y <= 1.0 ) THEN
erf = x * ( 1. + csevl ( (2.*x*x - 1.), erfcs, nterf ) )
ELSE IF ( y <= xbig ) THEN
erf = sign ( 1. - erfc ( y ), x )
ELSE
erf = sign ( 1., x )
END IF
RETURN
END FUNCTION erf
!
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION ERFC ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION erfc ( x )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! erfc(x) calculates the single precision complementary error
! function for single precision argument x.
!
!
! Reference:
!
! David Kahaner, Clever Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1988.
!
!
! series for erf on the interval 0. to 1.00000d+00
! with weighted error 7.10e-18
! log weighted error 17.15
! significant figures required 16.31
! decimal places required 17.71
!
! series for erfc on the interval 0. to 2.50000d-01
! with weighted error 4.81e-17
! log weighted error 16.32
! approx significant figures required 15.0E+00
!
!
! series for erc2 on the interval 2.50000d-01 to 1.00000d+00
! with weighted error 5.22e-17
! log weighted error 16.28
! approx significant figures required 15.0E+00
! decimal places required 16.96
!
!
! AUTHOR: John Burkardt, Pittsburgh Supercomputing Center
! NMS Fortran Software Package
!
! MODIFICATION HISTORY:
! 23-Oct-2004 Modified for ARPS Standard Fortran-90
! Keith Brewster, CAPS Univ of Oklahoma
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL :: erfc
REAL :: x
!
REAL :: csevl
REAL :: r1mach
!
REAL :: erfcs(13)
REAL :: erfccs(24)
REAL :: erc2cs(23)
REAL :: eta
INTEGER :: inits
INTEGER :: nterc2
INTEGER :: nterf
INTEGER :: nterfc
REAL :: sqeps
REAL :: sqrtpi
REAL :: xmax
REAL :: xsml
REAL :: y
!
data erfcs( 1) / -.049046121234691808e0 /
data erfcs( 2) / -.14226120510371364e0 /
data erfcs( 3) / .010035582187599796e0 /
data erfcs( 4) / -.000576876469976748e0 /
data erfcs( 5) / .000027419931252196e0 /
data erfcs( 6) / -.000001104317550734e0 /
data erfcs( 7) / .000000038488755420e0 /
data erfcs( 8) / -.000000001180858253e0 /
data erfcs( 9) / .000000000032334215e0 /
data erfcs(10) / -.000000000000799101e0 /
data erfcs(11) / .000000000000017990e0 /
data erfcs(12) / -.000000000000000371e0 /
data erfcs(13) / .000000000000000007e0 /
data erc2cs( 1) / -.069601346602309501e0 /
data erc2cs( 2) / -.041101339362620893e0 /
data erc2cs( 3) / .003914495866689626e0 /
data erc2cs( 4) / -.000490639565054897e0 /
data erc2cs( 5) / .000071574790013770e0 /
data erc2cs( 6) / -.000011530716341312e0 /
data erc2cs( 7) / .000001994670590201e0 /
data erc2cs( 8) / -.000000364266647159e0 /
data erc2cs( 9) / .000000069443726100e0 /
data erc2cs(10) / -.000000013712209021e0 /
data erc2cs(11) / .000000002788389661e0 /
data erc2cs(12) / -.000000000581416472e0 /
data erc2cs(13) / .000000000123892049e0 /
data erc2cs(14) / -.000000000026906391e0 /
data erc2cs(15) / .000000000005942614e0 /
data erc2cs(16) / -.000000000001332386e0 /
data erc2cs(17) / .000000000000302804e0 /
data erc2cs(18) / -.000000000000069666e0 /
data erc2cs(19) / .000000000000016208e0 /
data erc2cs(20) / -.000000000000003809e0 /
data erc2cs(21) / .000000000000000904e0 /
data erc2cs(22) / -.000000000000000216e0 /
data erc2cs(23) / .000000000000000052e0 /
data erfccs( 1) / 0.0715179310202925e0 /
data erfccs( 2) / -.026532434337606719e0 /
data erfccs( 3) / .001711153977920853e0 /
data erfccs( 4) / -.000163751663458512e0 /
data erfccs( 5) / .000019871293500549e0 /
data erfccs( 6) / -.000002843712412769e0 /
data erfccs( 7) / .000000460616130901e0 /
data erfccs( 8) / -.000000082277530261e0 /
data erfccs( 9) / .000000015921418724e0 /
data erfccs(10) / -.000000003295071356e0 /
data erfccs(11) / .000000000722343973e0 /
data erfccs(12) / -.000000000166485584e0 /
data erfccs(13) / .000000000040103931e0 /
data erfccs(14) / -.000000000010048164e0 /
data erfccs(15) / .000000000002608272e0 /
data erfccs(16) / -.000000000000699105e0 /
data erfccs(17) / .000000000000192946e0 /
data erfccs(18) / -.000000000000054704e0 /
data erfccs(19) / .000000000000015901e0 /
data erfccs(20) / -.000000000000004729e0 /
data erfccs(21) / .000000000000001432e0 /
data erfccs(22) / -.000000000000000439e0 /
data erfccs(23) / .000000000000000138e0 /
data erfccs(24) / -.000000000000000048e0 /
data nterf, nterfc, nterc2, xsml, xmax, sqeps /3*0, 3*0./
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
sqrtpi=sqrt(acos(-1.))
IF ( nterf == 0 ) THEN
eta = 0.1 * r1mach(3)
nterf = inits ( erfcs, 13, eta )
nterfc = inits ( erfccs, 24, eta )
nterc2 = inits ( erc2cs, 23, eta )
xsml = -sqrt ( - log ( sqrtpi * r1mach(3) ) )
xmax = sqrt ( - log ( sqrtpi * r1mach(1) ) )
xmax = xmax - 0.5 * log ( xmax ) / xmax - 0.01
sqeps = sqrt ( 2. * r1mach(3) )
END IF
IF ( x <= xsml ) THEN
erfc = 2.
RETURN
END IF
IF ( x > xmax ) THEN
WRITE(6,'(a)') 'erfc x so big erfc underflows'
erfc = 0.
RETURN
END IF
y = abs(x)
IF (y>1.0) THEN
!
! erfc(x) = 1.0E+00 - erf(x) for 1. < abs(x) <= xmax
!
y = y*y
IF (y<=4.) THEN
erfc = exp(-y)/abs(x) * (0.5 + csevl ((8./y-5.)/3.,erc2cs, nterc2) )
ELSE
erfc = exp(-y)/abs(x) * (0.5 + csevl (8./y-1.,erfccs, nterfc) )
END IF
IF ( x < 0.0 ) erfc = 2.0 - erfc
RETURN
ELSE
!
! erfc(x) = 1. - erf(x) for -1. <= x <= 1.
!
IF ( y < sqeps ) THEN
erfc = 1.0 - 2.0*x/sqrtpi
ELSE IF (y>=sqeps) THEN
erfc = 1.0 - x*(1.0 + csevl((2.0*x*x-1.0), erfcs, nterf) )
END IF
RETURN
END IF
END FUNCTION erfc
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION CSEVL ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION csevl ( x, cs, n )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! CSEVL evaluates an N term Chebyshev series.
!
!
! Reference:
!
! R Broucke,
! algorithm 446, c.a.c.m.,
! volume 16, page 254, 1973.
!
! Fox and Parker,
! chebyshev polynomials in numerical analysis,
! oxford press, page 56.
!
! David Kahaner, Clever Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1988.
!
! INPUT:
!
! x value at which the series is to be evaluated.
! cs array of n terms of a chebyshev series. in eval-
! uating cs, only half the first coefficient is summed.
! n number of terms in array cs.
!
!
! AUTHOR: John Burkardt, Pittsburgh Supercomputing Center
! NMS Fortran Software Package
!
! MODIFICATION HISTORY:
! 23-Oct-2004 Modified for ARPS Standard Fortran-90
! Keith Brewster, CAPS Univ of Oklahoma
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL :: csevl
INTEGER :: n
REAL :: x
REAL :: cs(n)
!
REAL :: b0
REAL :: b1
REAL :: b2
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( n < 1 ) THEN
WRITE ( *, * ) ' '
WRITE ( *, * ) 'CSEVL - Fatal error!'
WRITE ( *, * ) ' Number of terms N is less than 1.'
STOP
END IF
IF ( n > 1000 ) THEN
WRITE ( *, * ) ' '
WRITE ( *, * ) 'CSEVL - Fatal error!'
WRITE ( *, * ) ' The number of terms is more than 1000.'
STOP
END IF
IF ( x < -1.0 .or. x > 1.0 ) THEN
WRITE ( *, * ) ' '
WRITE ( *, * ) 'CSEVL - Fatal error!'
WRITE ( *, * ) ' The input argument X is outside the interval [-1,1].'
STOP
END IF
b1 = 0.
b0 = 0.
DO i = n, 1, -1
b2 = b1
b1 = b0
b0 = 2.*x*b1 - b2 + cs(i)
END DO
csevl = 0.5*(b0-b2)
RETURN
END FUNCTION csevl
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION INITS ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION inits ( os, nos, eta )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!
! INITS estimates the order of an orthogonal series guaranteeing a given accuracy.
!
!
! Reference:
!
! David Kahaner, Clever Moler, Steven Nash,
! Numerical Methods and Software,
! Prentice Hall, 1988.
!
!
! INPUT:
!
! os array of nos coefficients in an orthogonal series.
! nos number of coefficients in os.
! eta requested accuracy of series.
! Ordinarily, eta will be chosen to be one-tenth machine precision.
!
!
! AUTHOR: John Burkardt, Pittsburgh Supercomputing Center
! NMS Fortran Software Package
!
! MODIFICATION HISTORY:
! 23-Oct-2004 Modified for ARPS Standard Fortran-90
! Keith Brewster, CAPS Univ of Oklahoma
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: inits
INTEGER :: nos
REAL :: eta
!
INTEGER :: ii,i
REAL :: err
REAL :: os(nos)
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( nos < 1) THEN
WRITE(6,'(a)') 'inits number of coefficients lt 1'
STOP
END IF
err = 0.
DO ii=1,nos
i = nos + 1 - ii
err = err + abs(os(i))
if (err>eta) EXIT
END DO
if (i==nos) WRITE(6,'(a)') 'inits eta may be too small'
inits = i
RETURN
END FUNCTION inits
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION ERF ######
!###### ######
!##################################################################
!##################################################################
!
!
FUNCTION r1mach ( i )
!
!*******************************************************************************
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! R1MACH returns single precision machine constants.
!
!
! Assume that single precision numbers are stored with a mantissa of T digits
! in base B, with an exponent whose value must lie between EMIN and EMAX. Then
! for values of I between 1 and 5, R1MACH will return the following values:
!
! R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
! R1MACH(2) = B**EMAX*(1-B**(-T)), the largest magnitude.
! R1MACH(3) = B**(-T), the smallest relative spacing.
! R1MACH(4) = B**(1-T), the largest relative spacing.
! R1MACH(5) = log10(B)
!
! To alter this function for a particular environment, the desired set of data
! statements should be activated by removing the C from column 1.
!
! On rare machines a STATIC statement may need to be added. But probably more
! systems prohibit it that require it.
!
! For IEEE-arithmetic machines (binary standard), the first set of constants
! below should be appropriate.
!
! Where possible, octal or hexadecimal constants have been used to specify the
! constants exactly which has in some cases required the use of EQUIVALENCED
! integer arrays. If your compiler uses half-word integers by default
! (sometimes called INTEGER*2), you may need to change INTEGER to INTEGER*4 or
! otherwise instruct your compiler to use full-word integers in the next 5
! declarations.
!
!
! AUTHOR: John Burkardt, Pittsburgh Supercomputing Center
! NMS Fortran Software Package
!
! MODIFICATION HISTORY:
! 23-Oct-2004 Modified for ARPS Standard Fortran-90
! Keith Brewster, CAPS Univ of Oklahoma
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: r1mach
INTEGER :: i
!
INTEGER :: diver(2)
INTEGER :: large(2)
INTEGER :: ilog10(2)
INTEGER :: right(2)
REAL :: rmach(5)
INTEGER ::small(2)
!
equivalence (rmach(1),small(1))
equivalence (rmach(2),large(1))
equivalence (rmach(3),right(1))
equivalence (rmach(4),diver(1))
equivalence (rmach(5),ilog10(1))
!
! IEEE arithmetic machines, such as the ATT 3B series, Motorola 68000 based
! machines such as the SUN 3 and ATT PC 7300, and 8087 based micros such as
! the IBM PC and ATT 6300.
!
! data small(1) / 8388608 /
! data large(1) / 2139095039 /
! data right(1) / 864026624 /
! data diver(1) / 872415232 /
! data ilog10(1) / 1050288283 /
!
! ALLIANT FX/8 UNIX Fortran compiler with the -r8 command line option. This
! option causes all variables declared with 'REAL' to be of type 'REAL*8' or
! DOUBLE PRECISION. This option does not override the 'REAL*4' declarations.
! These R1MACH numbers below and the coresponding I1MACH are simply the DOUBLE
! PRECISION or 'REAL*8' numbers. If you use the -r8 your whole code (and the
! user libraries you link with, the system libraries are taken care of
! automagicly) must be compiled with this option.
!
! data rmach(1) / 2.22507385850721D-308 /
! data rmach(2) / 1.79769313486231D+308 /
! data rmach(3) / 1.1101827117665D-16 /
! data rmach(4) / 2.2203654423533D-16 /
! data rmach(5) / 3.01029995663981E-1 /
!
! AMDAHL machines.
!
! data small(1) / 1048576 /
! data large(1) / 2147483647 /
! data right(1) / 990904320 /
! data diver(1) / 1007681536 /
! data ilog10(1) / 1091781651 /
!
! BURROUGHS 1700 system.
!
! data rmach(1) / Z400800000 /
! data rmach(2) / Z5FFFFFFFF /
! data rmach(3) / Z4E9800000 /
! data rmach(4) / Z4EA800000 /
! data rmach(5) / Z500E730E8 /
!
! BURROUGHS 5700/6700/7700 systems.
!
! data rmach(1) / O1771000000000000 /
! data rmach(2) / O0777777777777777 /
! data rmach(3) / O1311000000000000 /
! data rmach(4) / O1301000000000000 /
! data rmach(5) / O1157163034761675 /
!
! CDC CYBER 170/180 series using NOS
!
! data rmach(1) / O"00014000000000000000" /
! data rmach(2) / O"37767777777777777777" /
! data rmach(3) / O"16404000000000000000" /
! data rmach(4) / O"16414000000000000000" /
! data rmach(5) / O"17164642023241175720" /
!
! CDC CYBER 170/180 series using NOS/VE
!
! data rmach(1) / Z"3001800000000000" /
! data rmach(2) / Z"4FFEFFFFFFFFFFFE" /
! data rmach(3) / Z"3FD2800000000000" /
! data rmach(4) / Z"3FD3800000000000" /
! data rmach(5) / Z"3FFF9A209A84FBCF" /
!
! CDC CYBER 200 series
!
! data rmach(1) / X'9000400000000000' /
! data rmach(2) / X'6FFF7FFFFFFFFFFF' /
! data rmach(3) / X'FFA3400000000000' /
! data rmach(4) / X'FFA4400000000000' /
! data rmach(5) / X'FFD04D104D427DE8' /
!
! CDC 6000/7000 series using FTN4.
!
! data rmach(1) / 00564000000000000000B /
! data rmach(2) / 37767777777777777776B /
! data rmach(3) / 16414000000000000000B /
! data rmach(4) / 16424000000000000000B /
! data rmach(5) / 17164642023241175720B /
!
! CDC 6000/7000 series using FTN5.
!
! data rmach(1) / O"00564000000000000000" /
! data rmach(2) / O"37767777777777777776" /
! data rmach(3) / O"16414000000000000000" /
! data rmach(4) / O"16424000000000000000" /
! data rmach(5) / O"17164642023241175720" /
!
! CONVEX C-1.
!
! data rmach(1) / '00800000'X /
! data rmach(2) / '7FFFFFFF'X /
! data rmach(3) / '34800000'X /
! data rmach(4) / '35000000'X /
! data rmach(5) / '3F9A209B'X /
!
! CONVEX C-120 (native mode) without -R8 option
!
! data rmach(1) / 2.9387360E-39 /
! data rmach(2) / 1.7014117E+38 /
! data rmach(3) / 5.9604645E-08 /
! data rmach(4) / 1.1920929E-07 /
! data rmach(5) / 3.0102999E-01 /
!
! CONVEX C-120 (native mode) with -R8 option
!
! data rmach(1) / 5.562684646268007D-309 /
! data rmach(2) / 8.988465674311577D+307 /
! data rmach(3) / 1.110223024625157D-016 /
! data rmach(4) / 2.220446049250313D-016 /
! data rmach(5) / 3.010299956639812D-001 /
!
! CONVEX C-120 (IEEE mode) without -R8 option
!
! data rmach(1) / 1.1754945E-38 /
! data rmach(2) / 3.4028234E+38 /
! data rmach(3) / 5.9604645E-08 /
! data rmach(4) / 1.1920929E-07 /
! data rmach(5) / 3.0102999E-01 /
!
! CONVEX C-120 (IEEE mode) with -R8 option
!
! data rmach(1) / 2.225073858507202D-308 /
! data rmach(2) / 1.797693134862315D+308 /
! data rmach(3) / 1.110223024625157D-016 /
! data rmach(4) / 2.220446049250313D-016 /
! data rmach(5) / 3.010299956639812D-001 /
!
! CRAY 1, 2, XMP and YMP.
!
! data rmach(1) / 200034000000000000000B /
! data rmach(2) / 577767777777777777776B /
! data rmach(3) / 377224000000000000000B /
! data rmach(4) / 377234000000000000000B /
! data rmach(5) / 377774642023241175720B /
!
! DATA GENERAL ECLIPSE S/200.
! Note - It may be appropriate to include the line: STATIC RMACH(5)
!
! data small /20K,0/
! data large /77777K,177777K/
! data right /35420K,0/
! data diver /36020K,0/
! data ilog10 /40423K,42023K/
!
! ELXSI 6400, assuming REAL*4 is the default real type.
!
! data small(1) / '00800000'X /
! data large(1) / '7F7FFFFF'X /
! data right(1) / '33800000'X /
! data diver(1) / '34000000'X /
! data ilog10(1) / '3E9A209B'X /
!
! HARRIS 220
!
! data small(1),small(2) / '20000000, '00000201 /
! data large(1),large(2) / '37777777, '00000177 /
! data right(1),right(2) / '20000000, '00000352 /
! data diver(1),diver(2) / '20000000, '00000353 /
! data ilog10(1),ilog10(2) / '23210115, '00000377 /
!
! HARRIS SLASH 6 and SLASH 7.
!
! data small(1),small(2) / '20000000, '00000201 /
! data large(1),large(2) / '37777777, '00000177 /
! data right(1),right(2) / '20000000, '00000352 /
! data diver(1),diver(2) / '20000000, '00000353 /
! data ilog10(1),ilog10(2) / '23210115, '00000377 /
!
! HONEYWELL DPS 8/70 and 600/6000 series.
!
! data rmach(1) / O402400000000 /
! data rmach(2) / O376777777777 /
! data rmach(3) / O714400000000 /
! data rmach(4) / O716400000000 /
! data rmach(5) / O776464202324 /
!
! HP 2100, 3 word double precision with FTN4
!
! data small(1), small(2) / 40000B, 1 /
! data large(1), large(2) / 77777B, 177776B /
! data right(1), right(2) / 40000B, 325B /
! data diver(1), diver(2) / 40000B, 327B /
! data ilog10(1), ilog10(2) / 46420B, 46777B /
!
! HP 2100, 4 word double precision with FTN4
!
! data small(1), small(2) / 40000B, 1 /
! data large91), large(2) / 77777B, 177776B /
! data right(1), right(2) / 40000B, 325B /
! data diver(1), diver(2) / 40000B, 327B /
! data ilog10(1), ilog10(2) / 46420B, 46777B /
!
! HP 9000
!
! r1mach(1) = 1.17549435E-38
! r1mach(2) = 1.70141163E+38
! r1mach(3) = 5.960464478E-8
! r1mach(4) = 1.119209290E-7
! r1mach(5) = 3.01030010E-1
!
! data small(1) / 00040000000B /
! data large(1) / 17677777777B /
! data right(1) / 06340000000B /
! data diver(1) / 06400000000B /
! data ilog10(1) / 07646420233B /
!
! IBM 360/370 series, XEROX SIGMA 5/7/9, SEL systems 85/86, PERKIN ELMER 3230,
! and PERKIN ELMER (INTERDATA) 3230.
!
! data rmach(1) / Z00100000 /
! data rmach(2) / Z7FFFFFFF /
! data rmach(3) / Z3B100000 /
! data rmach(4) / Z3C100000 /
! data rmach(5) / Z41134413 /
!
! IBM PC - Microsoft FORTRAN
!
! data small(1) / #00800000 /
! data large(1) / #7F7FFFFF /
! data right(1) / #33800000 /
! data diver(1) / #34000000 /
! data ilog10(1) / #3E9A209A /
!
! IBM PC - Professional FORTRAN and Lahey FORTRAN
!
! data small(1)/ Z'00800000' /
! data large(1)/ Z'7F7FFFFF' /
! data right(1)/ Z'33800000' /
! data diver(1)/ Z'34000000' /
! data ilog10(1)/ Z'3E9A209A' /
!
! INTERDATA 8/32 with the UNIX system FORTRAN 77 compiler.
! For the INTERDATA FORTRAN VII compiler replace the Z'S specifying HEX
! constants with Y'S.
!
! data rmach(1) / Z'00100000' /
! data rmach(2) / Z'7EFFFFFF' /
! data rmach(3) / Z'3B100000' /
! data rmach(4) / Z'3C100000' /
! data rmach(5) / Z'41134413' /
!
! PDP-10 (KA or KI processor).
!
! data rmach(1) / "000400000000 /
! data rmach(2) / "377777777777 /
! data rmach(3) / "146400000000 /
! data rmach(4) / "147400000000 /
! data rmach(5) / "177464202324 /
!
! PDP-11 FORTRANS supporting 32-bit integers (integer version).
!
! data small(1) / 8388608 /
! data large(1) / 2147483647 /
! data right(1) / 880803840 /
! data diver(1) / 889192448 /
! data ilog10(1) / 1067065499 /
!
! PDP-11 FORTRANS supporting 32-bit integers (octal version).
!
! data rmach(1) / O00040000000 /
! data rmach(2) / O17777777777 /
! data rmach(3) / O06440000000 /
! data rmach(4) / O06500000000 /
! data rmach(5) / O07746420233 /
!
! PDP-11 FORTRANS supporting 16-bit integers (integer version).
!
! data small(1),small(2) / 128, 0 /
! data large(1),large(2) / 32767, -1 /
! data right(1),right(2) / 13440, 0 /
! data diver(1),diver(2) / 13568, 0 /
! data ilog10(1),ilog10(2) / 16282, 8347 /
!
! PDP-11 FORTRANS supporting 16-bit integers (octal version).
!
! data small(1),small(2) / O000200, O000000 /
! data large(1),large(2) / O077777, O177777 /
! data right(1),right(2) / O032200, O000000 /
! data diver(1),diver(2) / O032400, O000000 /
! data ilog10(1),ilog10(2) / O037632, O020233 /
!
! SEQUENT BALANCE 8000.
!
! data small(1) / $00800000 /
! data large(1) / $7F7FFFFF /
! data right(1) / $33800000 /
! data diver(1) / $34000000 /
! data ilog10(1) / $3E9A209B /
!
! SUN Microsystems UNIX F77 compiler.
!
data rmach(1) / 1.17549435E-38 /
data rmach(2) / 3.40282347E+38 /
data rmach(3) / 5.96016605E-08 /
data rmach(4) / 1.19203321E-07 /
data rmach(5) / 3.01030010E-01 /
!
! SUN 3 (68881 or FPA)
!
! data small(1) / X'00800000' /
! data large(1) / X'7F7FFFFF' /
! data right(1) / X'33800000' /
! data diver(1) / X'34000000' /
! data ilog10(1) / X'3E9A209B' /
!
! UNIVAC 1100 series.
!
! data rmach(1) / O000400000000 /
! data rmach(2) / O377777777777 /
! data rmach(3) / O146400000000 /
! data rmach(4) / O147400000000 /
! data rmach(5) / O177464202324 /
!
! VAX/ULTRIX F77 compiler.
!
! data small(1) / 128 /
! data large(1) / -32769 /
! data right(1) / 13440 /
! data diver(1) / 13568 /
! data ilog10(1) / 547045274 /
!
! VAX-11 with FORTRAN IV-PLUS compiler.
!
! data rmach(1) / Z00000080 /
! data rmach(2) / ZFFFF7FFF /
! data rmach(3) / Z00003480 /
! data rmach(4) / Z00003500 /
! data rmach(5) / Z209B3F9A /
!
! VAX/VMS version 2.2.
!
! data rmach(1) / '80'X /
! data rmach(2) / 'FFFF7FFF'X /
! data rmach(3) / '3480'X /
! data rmach(4) / '3500'X /
! data rmach(5) / '209B3F9A'X /
!
! VAX/VMS 11/780
!
! data small(1) / Z00000080 /
! data large(1) / ZFFFF7FFF /
! data right(1) / Z00003480 /
! data diver(1) / Z00003500 /
! data ilog10(1) / Z209B3F9A /
!
! Z80 microprocessor.
!
! data small(1), small(2) / 0, 256 /
! data large(1), large(2) / -1, -129 /
! data right(1), right(2) / 0, 26880 /
! data diver(1), diver(2) / 0, 27136 /
! data ilog10(1), ilog10(2) / 8347, 32538 /
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( i > 0 .AND. i < 6 ) THEN
r1mach = rmach(i)
ELSE
WRITE(6,'(/a)') 'R1MACH - Fatal error!'
WRITE(6,'(a,i9)')'I is out of bounds=',i
r1mach=0.
STOP
END IF
RETURN
END FUNCTION r1mach
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GET_GRIDXYZP_HDF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE get_gridxyzp_hdf(nx,ny,nz,filename, & 1,12
x,y,zp, &
itmp,hmax,hmin,ireturn)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid variables from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, after a similar program developed by Ming Xue
! 02/06/2005.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! filename HDF file name of grid/base file.
!
! OUTPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! x x-coordinate data
! y y-coordinate data
! zp zp-coordinate data
!
! SCRATCH:
! itmp Temporary array for hdf compression
! hmin hmin temporary array
! hmax hmax temporary array
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: stat, sd_id
CHARACTER (LEN=*) :: filename
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: x(nx)
REAL :: y(ny)
REAL :: zp(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmin(nz)
REAL :: hmax(nz)
INTEGER :: ireturn ! Return status indicator
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
INTEGER :: istat
REAL :: alatpro(2)
REAL :: sclf,dxscl,dyscl,ctrx,ctry,swx,swy
LOGICAL :: success
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'grid.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "get_gridxyzprho_hdf: ERROR opening ", &
trim(filename)," for reading."
ireturn=-1
RETURN
ELSE
WRITE(6,*) 'File ',filename,' opened.'
END IF
success=.true.
CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
CALL hdfrdr(sd_id,"trulon",trulon,istat)
CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)
CALL hdfrd1d(sd_id,"x",nx,x,istat)
IF (istat /= 0) success=.false.
CALL hdfrd1d(sd_id,"y",ny,y,istat)
IF (istat /= 0) success=.false.
CALL hdfrd3d(sd_id,"zp",nx,ny,nz,zp,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
IF (success) THEN
ireturn = 0
ELSE
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
WRITE(6,'(/a/)') ' Error reading data in GET_GRIDXYZPRHO_HDF.'
ireturn=-2
END IF
CALL hdfclose(sd_id,stat)
RETURN
END SUBROUTINE get_gridxyzp_hdf
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREADUVW ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreaduvw(nx,ny,nz,filename, & 1,28
time,u,v,w, &
itmp,hmax,hmin,ireturn)
!-----------------------------------------------------------------------
! PURPOSE:
! Read in ARPS wind data in the NCSA HDF4 format.
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 2000/04/15
!
! MODIFICATION HISTORY:
! Keith Brewster, CAPS
! 2005/02/06 Created hdfreaduvw, streamlined from hdfreadwind
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! filename Character variable naming the input HDF file
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
CHARACTER (LEN=*) :: filename
REAL :: time
REAL :: u(nx,ny,nz)
REAL :: v(nx,ny,nz)
REAL :: w(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmin(nz) ! Temporary array
REAL :: hmax(nz) ! Temporary array
INTEGER :: ireturn
!-----------------------------------------------------------------------
! Parameters describing routine that wrote the gridded data
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version.
! intver??: an integer to allow faster comparison than fmtver??,
! which are strings.
!
! Verion 5.00: significant change in soil variables since version 4.10.
!
!-----------------------------------------------------------------------
CHARACTER (LEN=40) :: fmtver410,fmtver500
INTEGER :: intver,intver410,intver500
PARAMETER (fmtver410='004.10 HDF4 Coded Data',intver410=410)
PARAMETER (fmtver500='005.00 HDF4 Coded Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
!-----------------------------------------------------------------------
! Misc. local variables
!-----------------------------------------------------------------------
INTEGER :: lchanl
PARAMETER (lchanl=6) ! Channel number for formatted printing.
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin
INTEGER :: istat, sd_id
INTEGER :: varflg, istatus
LOGICAL :: success
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'indtflg.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Beginning of executable code...
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(*,*) 'HDFREADUVW: Reading HDF file: ', trim(filename)
!-----------------------------------------------------------------------
! Open file for reading
!-----------------------------------------------------------------------
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "HDFREADUVW: ERROR opening ", trim(filename)," for reading."
ireturn=-1
RETURN
END IF
fmtverin = fmtver500
WRITE(6,'(/1x,a,a/)') 'Incoming data format, fmtverin=',fmtverin
CALL hdfrdc(sd_id,80,"runname",runname,istat)
CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
IF( nocmnt > 0 ) THEN
CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
END IF
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') trim(runname)
WRITE (6,*) "Comments:"
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE (6,*) " "
CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
CALL hdfrdr(sd_id,"time",time,istat)
!-----------------------------------------------------------------------
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdi(sd_id,"nz",nzin,istat)
IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') ' Dimensions in HDFREADUVW inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz
WRITE(6,'(1x,a)') ' Program aborted in HDFREAD.'
CALL arpsstop('arpsstop called from HDFREADUVW due to nxin...',1)
END IF
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', &
time,' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",grdin,istat)
CALL hdfrdi(sd_id,"basflg",basin,istat)
CALL hdfrdi(sd_id,"varflg",varin,istat)
WRITE(6,'(a)') ' Done reading parameters'
CALL hdfrdi(sd_id,"month",month,istat)
CALL hdfrdi(sd_id,"day",day,istat)
CALL hdfrdi(sd_id,"year",year,istat)
CALL hdfrdi(sd_id,"hour",hour,istat)
CALL hdfrdi(sd_id,"minute",minute,istat)
CALL hdfrdi(sd_id,"second",second,istat)
CALL hdfrdr(sd_id,"umove",umove,istat)
CALL hdfrdr(sd_id,"vmove",vmove,istat)
CALL hdfrdr(sd_id,"xgrdorg",xgrdorg,istat)
CALL hdfrdr(sd_id,"ygrdorg",ygrdorg,istat)
success=.true.
IF( varin == 1 ) then
!-----------------------------------------------------------------------
! Read in total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,"u",nx,ny,nz,u,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
CALL hdfrd3d(sd_id,"v",nx,ny,nz,v,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
CALL hdfrd3d(sd_id,"w",nx,ny,nz,w,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
!-----------------------------------------------------------------------
!
! Friendly exit message
!
!-----------------------------------------------------------------------
IF(success) THEN
WRITE(6,'(/a,F8.1,a/)') &
' Data at time=', time/60,' (min) were successfully read.'
ireturn = 0
ELSE
WRITE(6,'(/a,F8.1,a/)') &
' Error reading u,v,w data at time=', time/60,' (min)'
ireturn = -3
END IF
CALL hdfclose(sd_id,istat)
ELSE
WRITE(6,'(/a/a,a/)') ' Error reading data in HDFREADUVW:', &
' No time-dependent data in file',trim(filename)
ireturn = -2
CALL hdfclose(sd_id,istat)
END IF
RETURN
END SUBROUTINE hdfreaduvw
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREADPPTRHO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreadppt(nx,ny,nz,filename, & 1,23
time,p,pt, &
itmp,hmax,hmin,ireturn)
!-----------------------------------------------------------------------
! PURPOSE:
! Read in ARPS pressure and temperature data in the NCSA HDF4 format.
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster after hdfreaduvw
! 02/06/2005
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! filename Character variable naming the input HDF file
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
CHARACTER (LEN=*) :: filename
REAL :: time
REAL :: p(nx,ny,nz)
REAL :: pt(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmin(nz) ! Temporary array
REAL :: hmax(nz) ! Temporary array
INTEGER :: ireturn
!-----------------------------------------------------------------------
! Parameters describing routine that wrote the gridded data
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version.
! intver??: an integer to allow faster comparison than fmtver??,
! which are strings.
!
! Verion 5.00: significant change in soil variables since version 4.10.
!
!-----------------------------------------------------------------------
CHARACTER (LEN=40) :: fmtver410,fmtver500
INTEGER :: intver,intver410,intver500
PARAMETER (fmtver410='004.10 HDF4 Coded Data',intver410=410)
PARAMETER (fmtver500='005.00 HDF4 Coded Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
!-----------------------------------------------------------------------
! Misc. local variables
!-----------------------------------------------------------------------
INTEGER :: lchanl
PARAMETER (lchanl=6) ! Channel number for formatted printing.
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin
INTEGER :: istat, sd_id
INTEGER :: varflg, istatus
LOGICAL :: success
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'indtflg.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Beginning of executable code...
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(*,*) 'HDFREADPPTRHO: Reading HDF file: ', trim(filename)
!-----------------------------------------------------------------------
! Open file for reading
!-----------------------------------------------------------------------
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,'(3a)') 'HDFREADPPTRHO: ERROR opening ', trim(filename), &
' for reading.'
ireturn=-1
RETURN
END IF
fmtverin = fmtver500
WRITE(6,'(/1x,a,a/)') 'Incoming data format, fmtverin=',fmtverin
CALL hdfrdc(sd_id,80,"runname",runname,istat)
CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
IF( nocmnt > 0 ) THEN
CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
END IF
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') trim(runname)
WRITE (6,*) "Comments:"
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE (6,*) " "
CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
CALL hdfrdr(sd_id,"time",time,istat)
!-----------------------------------------------------------------------
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdi(sd_id,"nz",nzin,istat)
IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') ' Dimensions in HDFREADUVW inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz
WRITE(6,'(1x,a)') ' Program aborted in HDFREAD.'
CALL arpsstop('arpsstop called from HDFREADUVW due to nxin...',1)
END IF
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', &
time,' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",grdin,istat)
CALL hdfrdi(sd_id,"basflg",basin,istat)
CALL hdfrdi(sd_id,"varflg",varin,istat)
WRITE(6,'(a)') ' Done reading parameters'
CALL hdfrdi(sd_id,"month",month,istat)
CALL hdfrdi(sd_id,"day",day,istat)
CALL hdfrdi(sd_id,"year",year,istat)
CALL hdfrdi(sd_id,"hour",hour,istat)
CALL hdfrdi(sd_id,"minute",minute,istat)
CALL hdfrdi(sd_id,"second",second,istat)
success=.true.
IF( varin == 1 ) THEN
!-----------------------------------------------------------------------
! Read in total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,"p",nx,ny,nz,p,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
CALL hdfrd3d(sd_id,"pt",nx,ny,nz,pt,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
!-----------------------------------------------------------------------
!
! Friendly exit message
!
!-----------------------------------------------------------------------
IF(success) THEN
WRITE(6,'(/a,F8.1,a/)') &
' Data at time=', time/60,' (min) were successfully read.'
ireturn = 0
ELSE
WRITE(6,'(/a,F8.1,a/)') &
' Error reading u,v,w data at time=', time/60,' (min)'
ireturn = -3
END IF
CALL hdfclose(sd_id,istat)
ELSE
WRITE(6,'(/a/a,a/)') ' Error reading data in HDFREADPPT:', &
' No time-dependent data in file',trim(filename)
ireturn = -2
CALL hdfclose(sd_id,istat)
END IF
RETURN
END SUBROUTINE hdfreadppt
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREADQRSH ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreadqrsh(nx,ny,nz,filename,time,qr,qs,qh, & 1,25
itmp,hmax,hmin,ireturn)
!-----------------------------------------------------------------------
! PURPOSE:
! Read in ARPS pressure and temperature data in the NCSA HDF4 format.
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster after hdfreaduvw
! 02/06/2005
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! filename Character variable naming the input HDF file
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
CHARACTER (LEN=*) :: filename
REAL :: time
REAL :: qr(nx,ny,nz)
REAL :: qs(nx,ny,nz)
REAL :: qh(nx,ny,nz)
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL :: hmin(nz) ! Temporary array
REAL :: hmax(nz) ! Temporary array
INTEGER :: ireturn
!-----------------------------------------------------------------------
! Parameters describing routine that wrote the gridded data
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version.
! intver??: an integer to allow faster comparison than fmtver??,
! which are strings.
!
! Verion 5.00: significant change in soil variables since version 4.10.
!
!-----------------------------------------------------------------------
CHARACTER (LEN=40) :: fmtver410,fmtver500
INTEGER :: intver,intver410,intver500
PARAMETER (fmtver410='004.10 HDF4 Coded Data',intver410=410)
PARAMETER (fmtver500='005.00 HDF4 Coded Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
!-----------------------------------------------------------------------
! Misc. local variables
!-----------------------------------------------------------------------
INTEGER :: lchanl
PARAMETER (lchanl=6) ! Channel number for formatted printing.
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin
INTEGER :: istat, sd_id
INTEGER :: varflg, istatus
LOGICAL :: success
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'indtflg.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Beginning of executable code...
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(*,*) 'HDFREADQRSH: Reading HDF file: ', trim(filename)
!-----------------------------------------------------------------------
! Open file for reading
!-----------------------------------------------------------------------
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "HDFREADQRSH: ERROR opening ", trim(filename)," for reading."
ireturn=-1
RETURN
END IF
fmtverin = fmtver500
WRITE(6,'(/1x,a,a/)') 'Incoming data format, fmtverin=',fmtverin
CALL hdfrdc(sd_id,80,"runname",runname,istat)
CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
IF( nocmnt > 0 ) THEN
CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
END IF
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') trim(runname)
WRITE (6,*) "Comments:"
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE (6,*) " "
CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
CALL hdfrdr(sd_id,"time",time,istat)
!-----------------------------------------------------------------------
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdi(sd_id,"nz",nzin,istat)
IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') ' Dimensions in HDFREADUVW inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz
WRITE(6,'(1x,a)') ' Program aborted in HDFREAD.'
CALL arpsstop('arpsstop called from HDFREADUVW due to nxin...',1)
END IF
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', &
time,' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",grdin,istat)
CALL hdfrdi(sd_id,"basflg",basin,istat)
CALL hdfrdi(sd_id,"varflg",varin,istat)
CALL hdfrdi(sd_id,"iceflg",icein,istat)
WRITE(6,'(a)') ' Done reading parameters'
CALL hdfrdi(sd_id,"month",month,istat)
CALL hdfrdi(sd_id,"day",day,istat)
CALL hdfrdi(sd_id,"year",year,istat)
CALL hdfrdi(sd_id,"hour",hour,istat)
CALL hdfrdi(sd_id,"minute",minute,istat)
CALL hdfrdi(sd_id,"second",second,istat)
success=.true.
IF( varin == 1 ) THEN
!-----------------------------------------------------------------------
! Read in total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,"qr",nx,ny,nz,qr,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
IF(icein == 1) THEN
CALL hdfrd3d(sd_id,"qs",nx,ny,nz,qs,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
CALL hdfrd3d(sd_id,"qh",nx,ny,nz,qh,istat,itmp,hmax,hmin)
IF (istat /= 0) success=.false.
ELSE
qs=0.
qh=0.
END IF
!-----------------------------------------------------------------------
!
! Friendly exit message
!
!-----------------------------------------------------------------------
IF(success) THEN
WRITE(6,'(/a,F8.1,a/)') &
' Data at time=', time/60,' (min) were successfully read.'
ireturn = 0
ELSE
WRITE(6,'(/a,F8.1,a/)') &
' Error reading qr,qs,qh data at time=', time/60,' (min)'
ireturn = -3
END IF
CALL hdfclose(sd_id,istat)
ELSE
WRITE(6,'(/a/a,a/)') ' Error reading data in HDFREADQRSH:', &
' No time-dependent data in file',trim(filename)
ireturn = -2
CALL hdfclose(sd_id,istat)
END IF
RETURN
END SUBROUTINE hdfreadqrsh
SUBROUTINE radverr(ngate,maxazim,nazim,radv,attrefl,stdvr,gatesp, & 1
vrerropt,sigmavr,samploi,vnyqstl, &
wavelen,pwrxmt,gaindb,lossdb,noisedbm,beamwid, &
pulselen,npulse,rmisval)
!
! Calculate an expected error for reflectivity based in the SNR
! of the attenuated signal and apply to the relectivity arrays.
!
! Calculation of expected radial velocity error based on SNR,
! estimated spectrum width and number of samples is derived
! from a fit to the modelled errors presented in Fig 6.5 of
! Doviak and Zrnic', 1993.
!
! Keith Brewster, CAPS
! 29 Sept 2005
!
IMPLICIT NONE
INTEGER :: ngate
INTEGER :: maxazim
INTEGER :: nazim
REAL :: radv(ngate,maxazim)
REAL :: attrefl(ngate,maxazim)
REAL :: stdvr(ngate,maxazim)
REAL :: gatesp
INTEGER :: vrerropt
REAL :: sigmavr
REAL :: samploi
REAL :: vnyqstl
REAL :: wavelen
REAL :: pwrxmt
REAL :: gaindb
REAL :: lossdb
REAL :: noisedbm
REAL :: beamwid
REAL :: pulselen
INTEGER :: npulse
REAL :: rmisval
!
! Misc. local variables
!
REAL, PARAMETER :: kw2=0.91
REAL :: taums,gatspkm,nplsind,nplsinv,twovnyq,twovninv
REAL :: pi,rkm,dBZ,ze,pwrmw,noisemw,snr,snrinv,pwrcst
REAL :: gain,gain2,loss,wavln2,beamwid2,pi5,ln2,two14,loss2,err
REAL :: stdvrnor,sdnor,sigmav,vr,sigmvcst,sigmvlim
INTEGER :: i,j,mid,ifold
REAL :: randnor
!
! Calculate constants, including power constant needed to
! compute SRM from reflectivity. Change units of radar parameters
! to those used in Eq. 4.35 of Doviak and Zrnic', 1993.
!
mid=(nazim/2)+1
pi=acos(-1.)
twovnyq=2.*vnyqstl
twovninv=1./twovnyq
sigmvcst=twovnyq/sqrt(npulse*samploi)
sigmvlim=0.8*vnyqstl
taums=1.0E06*pulselen
gatspkm=0.001*gatesp
nplsind=samploi*npulse
nplsinv=1./nplsind
noisemw=10.**(0.1*noisedbm)
gain=10.**(0.1*gaindb)
loss=10.**(0.1*lossdb)
gain2=gain*gain
wavln2=wavelen*wavelen
beamwid2=beamwid*beamwid
pi5=pi**5
two14=2.0**14
ln2=log(2.)
pwrcst=(pi5*1.0E-17*pwrxmt*gain2*taums*beamwid2*kw2)/ &
(6.75*two14*ln2*wavln2*loss)
IF(vrerropt == 1) THEN
DO j=1,nazim
DO i=1,ngate
IF(radv(i,j) > rmisval) THEN
vr=radv(i,j)+sigmavr*randnor()
ifold=NINT(vr*twovninv)
radv(i,j)=vr-ifold*twovnyq
END IF
END DO
END DO
ELSE IF(vrerropt > 1) THEN
WRITE(6,'(a,f12.2)') ' vnyqstl=',vnyqstl
WRITE(6,'(a,f12.0,a,f12.3)') ' gatesp=',gatesp,'gatspkm',gatspkm
! OPEN(51,file='snrp20.txt',status='unknown',form='formatted')
! WRITE(51,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
! OPEN(52,file='snrp15.txt',status='unknown',form='formatted')
! WRITE(52,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
! OPEN(53,file='snrp08.txt',status='unknown',form='formatted')
! WRITE(53,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
! OPEN(54,file='snrp00.txt',status='unknown',form='formatted')
! WRITE(54,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
! OPEN(55,file='snrm05.txt',status='unknown',form='formatted')
! WRITE(55,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
! OPEN(56,file='snrm06.txt',status='unknown',form='formatted')
! WRITE(56,'(6x,a)') &
! 'rkm attrefl snr stdvr sigmav vr'
DO j=1,nazim
DO i=1,ngate
IF(radv(i,j) > rmisval) THEN
IF(attrefl(i,j) > rmisval ) THEN
rkm=i*gatspkm
ze=10.**(attrefl(i,j)*0.1)
pwrmw=pwrcst*ze/(rkm*rkm)
snr=pwrmw/noisemw
ELSE
rkm=i*gatspkm
snr=0.01
END IF
snr=max(snr,0.01)
snrinv=1.0/snr
stdvrnor=stdvr(i,j)*twovninv
stdvrnor=min(max(stdvrnor,0.01),0.4)
sdnor=(0.0177*exp(12.0*stdvrnor)) + &
((0.159-0.755*stdvrnor+6.5*stdvrnor*stdvrnor)*snrinv)
sigmav=sigmvcst*sdnor
sigmav=min(max(sigmav,0.1),sigmvlim)
! IF(snr > 100.) THEN ! 20 dB
! WRITE(51,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvrnor,sdnor,sigmav
! ELSE IF(snr > 15.85) THEN ! 12 dB
! WRITE(52,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvrnor,sdnor,sigmav
! ELSE IF(snr > 3.16) THEN ! 5 dB
! WRITE(53,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvrnor,sdnor,sigmav
! ELSE IF(snr > 0.631) THEN ! -2 dB
! WRITE(54,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvrnor,sdnor,sigmav
! ELSE IF(snr > 0.20) THEN ! -7 dB
! WRITE(55,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvrnor,sdnor,sigmav
! ELSE
! WRITE(56,'(f10.1,2x,g10.5,2x,g10.5,2x,g10.5,2x,g10.4,f10.2)') &
! rkm,attrefl(i,j),snr,stdvr(i,j),sigmav,radv(i,j)
! END IF
vr=radv(i,j)+sigmav*randnor()
ifold=NINT(vr*twovninv)
radv(i,j)=vr-ifold*twovnyq
END IF
END DO
END DO
END IF
RETURN
END SUBROUTINE radverr
SUBROUTINE reflerr(ngate,maxazim,nazim,refl,attrefl,gatesp, & 1
rferropt,sigmarf,refmin,samploi, &
wavelen,pwrxmt,gaindb,lossdb,noisedbm,beamwid, &
pulselen,npulse,rmisval)
!
! Calculate an expected error for reflectivity based in the SNR
! of the attenuated signal and apply to the relectivity arrays.
!
! Keith Brewster, CAPS
! 29 Sept 2005
!
IMPLICIT NONE
INTEGER :: ngate
INTEGER :: maxazim
INTEGER :: nazim
REAL :: refl(ngate,maxazim)
REAL :: attrefl(ngate,maxazim)
REAL :: gatesp
INTEGER :: rferropt
REAL :: sigmarf
REAL :: refmin
REAL :: samploi
REAL :: wavelen
REAL :: pwrxmt
REAL :: gaindb
REAL :: lossdb
REAL :: noisedbm
REAL :: beamwid
REAL :: pulselen
INTEGER :: npulse
REAL :: rmisval
!
! Misc. local variables
!
REAL, PARAMETER :: kw2=0.91
REAL :: taums,gatspkm,nplsind,nplsinv,rminze
REAL :: pi,rkm,rkm2,dBZ,ze,pwrmw,noisemw,randerr
REAL :: snr,sigmaz,snrinv,pwrcst,pcstinv
REAL :: gain,gain2,loss,wavln2,beamwid2,pi5,ln2,two14,loss2,err
INTEGER :: i,j,mid
REAL :: randnor
WRITE(6,'(a,i4)') ' Calculating Reflectivity Error, Opt=',rferropt
mid=(nazim/2)+1
pi=acos(-1.)
taums=1.0E06*pulselen
gatspkm=0.001*gatesp
nplsind=samploi*npulse
nplsinv=1./nplsind
noisemw=10.**(0.1*noisedbm)
gain=10.**(0.1*gaindb)
loss=10.**(0.1*lossdb)
gain2=gain*gain
wavln2=wavelen*wavelen
beamwid2=beamwid*beamwid
pi5=pi**5
two14=2.0**14
ln2=log(2.)
pwrcst=(pi5*1.0E-17*pwrxmt*gain2*taums*beamwid2*kw2)/ &
(6.75*two14*ln2*wavln2*loss)
pcstinv=1.0/pwrcst
IF(rferropt == 1) THEN
DO j=1,nazim
DO i=1,ngate
IF(attrefl(i,j) > rmisval) THEN
err=sigmarf*randnor()
refl(i,j)=max(refmin,(refl(i,j)+err))
attrefl(i,j)=max(refmin,(attrefl(i,j)+err))
END IF
END DO
END DO
ELSE IF(rferropt > 1) THEN
! WRITE(6,'(6x,a)') 'rkm attrefl snr sigmaz'
DO j=1,nazim
DO i=1,ngate
IF(attrefl(i,j) > rmisval) THEN
rkm=i*gatspkm
rkm2=rkm*rkm
ze=10.**(attrefl(i,j)*0.1)
pwrmw=pwrcst*ze/rkm2
snr=pwrmw/noisemw
snrinv=1.0/snr
sigmaz= &
1.0+sqrt(nplsinv*((1.+snrinv)*(1.+snrinv)+(snrinv*snrinv)))
! IF(mod(i,100) == 0) &
! WRITE(6,'(/f10.1,2x,g10.4,2x,g10.4,2x,g10.4)') &
! rkm,attrefl(i,j),snr,sigmaz
randerr=randnor()
err=sigmaz*noisemw*randerr
pwrmw=pwrmw+err
IF (pwrmw > 0. ) THEN
ze=rkm2*pcstinv*pwrmw
attrefl(i,j)=10.*log10(ze)
attrefl(i,j)=max(refmin,attrefl(i,j))
ELSE
attrefl(i,j)=refmin
END IF
! IF(mod(i,100) == 0) THEN
! WRITE(6,'(12x,g10.4,2x,g10.4)') &
! attrefl(i,j),pwrmw
! print *, 'randerr, err=',randerr,err
! END IF
ze=10.**(refl(i,j)*0.1)
pwrmw=pwrcst*ze/(rkm*rkm)
pwrmw=pwrmw+err
IF (pwrmw > 0.) THEN
ze=rkm2*pcstinv*pwrmw
refl(i,j)=10.*log10(ze)
refl(i,j)=max(refmin,refl(i,j))
ELSE
refl(i,j)=refmin
END IF
END IF
END DO
END DO
END IF
RETURN
END SUBROUTINE reflerr
SUBROUTINE fake_reflec(nx,ny,nz,ref)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: ref(nx,ny,nz)
INTEGER :: i,j,k
REAL, PARAMETER :: refcst=30.0
WRITE(6,'(//a,/,/a,f10.2//)') '*******WARNING***************', &
' Test code: Setting all reflectivity in domain to:',refcst
DO k=1,nz
DO j=1,ny
DO i=1,nx
ref(i,j,k)=refcst
END DO
END DO
END DO
RETURN
END SUBROUTINE fake_reflec
SUBROUTINE fake_vel(nx,ny,nz,usc,vsc)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: usc(nx,ny,nz)
REAL :: vsc(nx,ny,nz)
INTEGER :: i,j,k
REAL, PARAMETER :: ucst=50.0
REAL, PARAMETER :: vcst=0.0
WRITE(6,'(//a,/,/a,f10.2,/a,f10.2//)') &
'*******WARNING***************', &
' Test code: Setting all u velocity in domain to:',ucst, &
' Test code: Setting all v velocity in domain to:',vcst
DO k=1,nz
DO j=1,ny
DO i=1,nx
usc(i,j,k)=ucst
vsc(i,j,k)=vcst
END DO
END DO
END DO
RETURN
END SUBROUTINE fake_vel