PROGRAM attenuation,38
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Sample program to read history data files produced by ARPS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 12/12/2003: Rewrote based on arpscvt. Input file names will be
! specified in arpsread.input.
!
! MODIFICATION HISTORY: edited arpsread for attenuation program
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nzsoil ! Grid dimensions.
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=256) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
! Arrays to be read in:
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x(:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y(:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z(:) ! The z-coord. of the computational
! grid. Defined at w-point.
REAL, ALLOCATABLE :: zp(:,:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The height of the terrain.
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 :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: 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)
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: qv (:,:,:) ! Water vapor specific humidity (kg/kg)
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type
INTEGER, ALLOCATABLE :: vegtyp(:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3)
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: radswnet (:,:) ! Net solar radiation, SWin - SWout
REAL, ALLOCATABLE :: radlwin (:,:) ! Incoming longwave radiation
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
REAL, ALLOCATABLE :: xsc(:)
REAL, ALLOCATABLE :: ysc(:)
REAL, ALLOCATABLE :: zps(:,:,:)
REAL, ALLOCATABLE :: latsc(:,:)
REAL, ALLOCATABLE :: lonsc(:,:)
REAL, ALLOCATABLE :: azmsc(:,:)
REAL, ALLOCATABLE :: sfcr(:,:)
REAL, ALLOCATABLE :: cmpref(:,:)
REAL, ALLOCATABLE :: elvsc(:,:,:)
REAL, ALLOCATABLE :: rngsc(:,:,:)
REAL, ALLOCATABLE :: reflect(:,:,:) ! reflectivity (dBZ)
REAL, ALLOCATABLE :: refz(:,:,:) ! reflectivity (mm6/m3)
!
!-----------------------------------------------------------------------
! Temporary working arrays:
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem4(:,:,:)
!
!-----------------------------------------------------------------------
! Az-ran weighting arrays:
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: delrng(:)
REAL, ALLOCATABLE :: delazm(:)
REAL, ALLOCATABLE :: delelv(:)
REAL, ALLOCATABLE :: wgtpt(:,:,:)
!
!-----------------------------------------------------------------------
! Radar radial arrays
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: refl(:)
REAL, ALLOCATABLE :: refatt(:)
REAL, ALLOCATABLE :: sens(:)
!
!-----------------------------------------------------------------------
! Hit-Miss bookkeeping arrays
!-----------------------------------------------------------------------
!
INTEGER, ALLOCATABLE :: score(:,:,:)
INTEGER, ALLOCATABLE :: netscore(:,:)
REAL, ALLOCATABLE :: tornlat(:,:)
REAL, ALLOCATABLE :: tornlon(:,:)
!
!-----------------------------------------------------------------------
! Original Observing Radar data variables
!-----------------------------------------------------------------------
!
CHARACTER (len=4) :: radname
REAL :: radlat,radlon,radelv,lowest_elv,bmwidth
INTEGER :: fillopt
NAMELIST /radar_src/ radname,radlat,radlon,radelv, &
lowest_elv,bmwidth,fillopt
!
!-----------------------------------------------------------------------
! Tornado location information
!-----------------------------------------------------------------------
!
REAL :: torn_azim,torn_rngkm,torn_refl
NAMELIST /tornado_loc/ torn_azim,torn_rngkm,torn_refl
!
!-----------------------------------------------------------------------
! Spatial offset variables
!-----------------------------------------------------------------------
!
INTEGER :: adjctr,adjhgt
REAL :: ctrlatnew,ctrlonnew,hgtoffset
REAL :: deltax,deltay
INTEGER :: ioffbgn,ioffend,joffbgn,joffend
NAMELIST /offsetxy/ adjctr,ctrlatnew,ctrlonnew, &
adjhgt,hgtoffset, &
deltax,deltay,ioffbgn,ioffend,joffbgn,joffend
!
!-----------------------------------------------------------------------
! Radar Network variables
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: mxradnet = 40
INTEGER :: nradnet,ngate,pulseopt,dualprf
INTEGER :: npulse1,npulse2
REAL :: elvmin,gatesp
REAL :: pulselen1,pulselen2,pulsewid1,pulsewid2,prt1,prt2
REAL :: rotrate,beamwid
REAL :: radnet_lat(mxradnet),radnet_lon(mxradnet)
REAL :: radnet_elev(mxradnet)
CHARACTER (len=4) :: radnet(mxradnet)
NAMELIST /radar_net/ nradnet,ngate,pulseopt,dualprf, &
elvmin,gatesp,pulselen1,pulselen2, &
pulsewid1,pulsewid2,prt1,prt2, &
npulse1,npulse2,beamwid,rotrate, &
radnet_lat,radnet_lon,radnet_elev,radnet
INTEGER :: nptsrng,nptsazm,nptselv
REAL :: evalwid
NAMELIST /emul_specs/ nptsrng,nptsazm,nptselv,evalwid
!
!-----------------------------------------------------------------------
! Misc. internal variables
!-----------------------------------------------------------------------
!
INTEGER, PARAMETER :: kntrmin=5
REAL, PARAMETER :: rngmin=1000.
REAL, PARAMETER :: rmisval=-99.
REAL, PARAMETER :: c = 2.9979e08 ! speed of light m/s
REAL, PARAMETER :: sumweps=1.0E-04
INTEGER :: nchin,idxdir
INTEGER :: grdbas
INTEGER :: i,j,k,ireturn
INTEGER :: idx,igate,iradius,jradius,irad
INTEGER :: ibgn,iend,jbgn,jend,ipt,jpt,ioff,joff
INTEGER :: iloc1,jloc1,iloc,jloc,ilc,jlc,klc,itgate
INTEGER :: kpt,npulsmn
INTEGER :: kntclr,knthob,kntvob,kntfewr,kntr,kntpt,kntintr
INTEGER :: length, ierr, istatus
INTEGER :: nt,nhits,nattn,nrng,kntrad
REAL :: alatnot(2)
REAL :: rad2deg,deg2rad
REAL :: twdx2inv,twdxinv,twdy2inv,twdyinv
REAL :: prtmn,pulselmn,pulsewmn
REAL :: r11l,r21l,r12l,r22l,r11h,r21h,r12h,r22h
REAL :: time,hgtmsl,bmwrad,bmwm,elv,azimuth,azrad
REAL :: depth,dsdr,srange,hgt,sradius,elradius
REAL :: hlfgtsp,efbmwid,fourln4,wasqinv,wesqinv,da,de,dr
REAL :: xsw,ysw,ctrx,ctry,xrad,yrad,xrada,yrada,radlata,radlona
REAL :: srgpt,azmpt,elvpt,hgtpt,sfcrpt
REAL :: rmax,dhdr,gatspkm,delx,dely
REAL :: gtlat,gtlon,gtx,gty,xrat,yrat,xpt,ypt
REAL :: ipos,jpos,irloc,jrloc
REAL :: w1i,w2i,w1j,w2j,w11,w12,w21,w22
REAL :: wgt0,wgt1,whigh,wlow,a,b,fjp1,fj0,fjm1
REAL :: flow,fhigh,reflz,refdbz
REAL :: reflzkm1,reflzk,reflz0,reflz1,reflz2
REAL :: pi,cinv,hlftau,aconst,b6,hlfctau,sigr2,sr2inv
REAL :: wgtg,wgtr,sumwr,sumr
REAL :: torn_rng,torn0x,torn0y,tornx,torny
REAL :: tornsfcr,tornhgt,tornelv,tornrngi,torngate
REAL :: tornref,tornaref,tornsen
REAL :: hpct,apct,rpct,hsum,asum
CHARACTER(len=24) :: varunits
CHARACTER(len=6) :: varid
CHARACTER(len=24) :: varname
CHARACTER(len=180), PARAMETER :: tornfile='tornscore.txt'
CHARACTER(len=180), PARAMETER :: tornstat='tornstat.txt'
LOGICAL :: echo
LOGICAL :: hit,rng
LOGICAL :: printit
!
!-----------------------------------------------------------------------
! Functions
!-----------------------------------------------------------------------
!
REAL :: effbmw
!
!-----------------------------------------------------------------------
! Variables used in finding height of lowest elv + half beamwidth
!-----------------------------------------------------------------------
!
REAL, PARAMETER :: refmiss = -90.
INTEGER :: klow
REAL :: radarx
REAL :: radary
REAL :: bmtop
REAL :: sfcrng
REAL :: bmhgt,bmhgtmsl
REAL :: reflow,hgtlow
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
pi=acos(-1.)
deg2rad=pi/180.
rad2deg=1./deg2rad
!
!-----------------------------------------------------------------------
! Get the names of the input data files.
!-----------------------------------------------------------------------
!
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = len_trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF (nstyps <= 0) nstyps = 1
nstyp = nstyps
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz ,'nzsoil=',nzsoil
ALLOCATE(x (nx))
ALLOCATE(y (ny))
ALLOCATE(z (nz))
ALLOCATE(zp (nx,ny,nz))
ALLOCATE(zpsoil (nx,ny,nzsoil))
ALLOCATE(uprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "uprt")
ALLOCATE(vprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "vprt")
ALLOCATE(wprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "wprt")
ALLOCATE(ptprt (nx,ny,nz))
ALLOCATE(pprt (nx,ny,nz))
ALLOCATE(qvprt (nx,ny,nz))
ALLOCATE(qc (nx,ny,nz))
ALLOCATE(qr (nx,ny,nz))
ALLOCATE(qi (nx,ny,nz))
ALLOCATE(qs (nx,ny,nz))
ALLOCATE(qh (nx,ny,nz))
ALLOCATE(tke (nx,ny,nz))
ALLOCATE(kmh (nx,ny,nz))
ALLOCATE(kmv (nx,ny,nz))
ALLOCATE(ubar (nx,ny,nz))
ALLOCATE(vbar (nx,ny,nz))
ALLOCATE(wbar (nx,ny,nz))
ALLOCATE(ptbar (nx,ny,nz))
ALLOCATE(pbar (nx,ny,nz))
ALLOCATE(rhobar (nx,ny,nz))
ALLOCATE(qvbar (nx,ny,nz))
ALLOCATE(u (nx,ny,nz))
ALLOCATE(v (nx,ny,nz))
ALLOCATE(w (nx,ny,nz))
ALLOCATE(qv (nx,ny,nz))
ALLOCATE(soiltyp (nx,ny,nstyps))
ALLOCATE(stypfrct(nx,ny,nstyps))
ALLOCATE(vegtyp (nx,ny))
ALLOCATE(lai (nx,ny))
ALLOCATE(roufns (nx,ny))
ALLOCATE(veg (nx,ny))
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps))
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps))
ALLOCATE(wetcanp(nx,ny,0:nstyps))
ALLOCATE(snowdpth(nx,ny))
ALLOCATE(raing(nx,ny))
ALLOCATE(rainc(nx,ny))
ALLOCATE(prcrate(nx,ny,4))
ALLOCATE(radsw (nx,ny))
ALLOCATE(rnflx (nx,ny))
ALLOCATE(radswnet (nx,ny))
ALLOCATE(radlwin (nx,ny))
ALLOCATE(radfrc(nx,ny,nz))
ALLOCATE(usflx (nx,ny))
ALLOCATE(vsflx (nx,ny))
ALLOCATE(ptsflx(nx,ny))
ALLOCATE(qvsflx(nx,ny))
ALLOCATE(tem1(nx,ny,nz))
ALLOCATE(tem2(nx,ny,nz))
ALLOCATE(tem3(nx,ny,nz))
ALLOCATE(tem4(nx,ny,nz))
ALLOCATE(xsc (nx))
ALLOCATE(ysc (ny))
ALLOCATE(latsc(nx,ny))
ALLOCATE(lonsc(nx,ny))
ALLOCATE(azmsc(nx,ny))
ALLOCATE(sfcr(nx,ny))
ALLOCATE(cmpref(nx,ny))
ALLOCATE(zps (nx,ny,nz))
ALLOCATE(reflect(nx,ny,nz))
ALLOCATE(refz(nx,ny,nz))
ALLOCATE(elvsc(nx,ny,nz))
ALLOCATE(rngsc(nx,ny,nz))
x =0.0
y =0.0
z =0.0
zp =0.0
zpsoil =0.0
uprt =0.0
vprt =0.0
wprt =0.0
ptprt =0.0
pprt =0.0
qvprt =0.0
qc =0.0
qr =0.0
qi =0.0
qs =0.0
qh =0.0
tke =0.0
kmh =0.0
kmv =0.0
ubar =0.0
vbar =0.0
wbar =0.0
ptbar =0.0
pbar =0.0
rhobar =0.0
qvbar =0.0
u =0.0
v =0.0
w =0.0
qv =0.0
soiltyp =0.0
stypfrct=0.0
vegtyp =0.0
lai =0.0
roufns =0.0
veg =0.0
tsoil =0.0
qsoil =0.0
wetcanp=0.0
snowdpth=0.0
raing=0.0
rainc=0.0
prcrate=0.0
radfrc=0.0
radsw =0.0
rnflx =0.0
radswnet = 0.0
radlwin = 0.0
usflx =0.0
vsflx =0.0
ptsflx=0.0
qvsflx=0.0
tem1=0.0
tem2=0.0
tem4=0.0
xsc=0.
ysc=0.
latsc=0.
lonsc=0.
azmsc=0.
sfcr=0.
cmpref=0.
zps=0.
reflect=0.
refz=0.
elvsc=0.
rngsc=0.
lengbf=len_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
lenfil=len_trim(hisfile(1))
WRITE(6,'(/a,a,a)') &
' Data set ', trim(hisfile(1)),' to be converted.'
!
!-----------------------------------------------------------------------
! Read all input data arrays
!-----------------------------------------------------------------------
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
hisfile(1)(1:lenfil),lenfil,time, &
x,y,z,zp,zpsoil, uprt ,vprt ,wprt ,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil, qsoil, wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1, tem2, tem3)
CALL gtlfnkey( runname, lfnkey )
idxdir=index(hisfile(1)(1:lenfil),runname(1:lfnkey),.true.)
IF(idxdir<=1) THEN
dirname='.'
ELSE
IF(hisfile(1)(idxdir-1:idxdir-1) == '/') THEN
dirname=hisfile(1)(1:idxdir-2)
ELSE
dirname=hisfile(1)(1:idxdir-1)
END IF
END IF
print *, ' ATTEN dirname= ',dirname
ldirnam=LEN(dirname)
CALL strlnth( dirname , ldirnam)
curtim = time
dx=(x(2)-x(1))
dy=(y(2)-y(1))
dxinv=1./dx
dyinv=1./dy
twdx2inv=1./(2.*dx*dx)
twdxinv=1./(2.*dx)
twdy2inv=1./(2.*dy*dy)
twdyinv=1./(2.*dy)
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)
if( ireturn /= 0 ) then
WRITE(6,'(1x,a,i2,/1x,a)') &
'Data read was unsuccessful. ireturn =', ireturn,' Job stopped.'
STOP
endif
!
radname='KTLX'
radlat=35.33
radlon=-97.28
radelv=350.0
lowest_elv=0.5
WRITE(6,'(a)') ' Reading radar_src namelist...'
READ(5,radar_src)
WRITE(6,'(a/)') ' Namelist radar_src successfully read.'
!
!-----------------------------------------------------------------------
! Get tornado location relative to original observing radar
!-----------------------------------------------------------------------
!
torn_azim=300.
torn_rngkm=65.
torn_refl=40.
WRITE(6,'(a)') ' Reading tornado_loc namelist...'
READ(5,tornado_loc)
WRITE(6,'(a/)') ' Namelist tornado_loc successfully read.'
torn_rng=1000.0*torn_rngkm
!
!-----------------------------------------------------------------------
! Get offset control variables
!-----------------------------------------------------------------------
!
adjctr=0
adjhgt=0
ctrlatnew=35.
ctrlonnew=-98.
hgtoffset=0.
deltax=3000.
deltay=3000.
ioffbgn=-10
ioffend=10
joffbgn=-10
joffend=10
WRITE(6,'(a)') ' Reading offsetxy namelist...'
READ(5,offsetxy)
WRITE(6,'(a/)') ' Namelist offsetxy successfully read.'
!
!-----------------------------------------------------------------------
! Get radar network observing locations
! and operating specs for network radar
!-----------------------------------------------------------------------
!
nradnet=1
ngate=150
elvmin=1.0
gatesp=200.
dualprf=1
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
beamwid = 2.0
rotrate=18.
evalwid=2.0
WRITE(6,'(a)') ' Reading radar_net namelist...'
READ(5,radar_net)
WRITE(6,'(a/)') ' Namelist radar_net successfully read.'
!
gatspkm=0.001*gatesp
hlfgtsp=0.5*gatesp
ALLOCATE(refl(ngate))
refl=0.
ALLOCATE(refatt(ngate))
refatt=0.
ALLOCATE(sens(ngate))
sens=0.
ALLOCATE(score(ioffbgn:ioffend,joffbgn:joffend,nradnet))
score=-99
ALLOCATE(netscore(ioffbgn:ioffend,joffbgn:joffend))
score=-99
ALLOCATE(tornlat(ioffbgn:ioffend,joffbgn:joffend))
tornlat=-999.
ALLOCATE(tornlon(ioffbgn:ioffend,joffbgn:joffend))
tornlon=-999.
!
CALL radcoord(nx,ny,nz, &
x,y,z,zp,xsc,ysc,zps, &
radlat,radlon,radarx,radary)
varid='reflct'
CALL readvar1(nx,ny,nz, reflect, varid, varname, &
varunits, time, runname, dirname, istatus)
elv=elvmin
IF(fillopt > 0 ) THEN
!
! Fill in radar reflectivity below the lowest beam height using
! zero gradient assumption.
!
bmtop=lowest_elv+0.5*bmwidth
DO j=1,ny
DO i=1,nx
!
! find height of lowest elev + half beamwidth
!
delx=xsc(i)-radarx
dely=ysc(j)-radary
sfcrng=sqrt(delx*delx+dely*dely)
CALL bmhgtsfr(bmtop,sfcrng,bmhgt)
bmhgtmsl=bmhgt+radelv
klow=nz
reflow=0.
hgtlow=1.0E06
DO k=nz-1,2,-1
IF(reflect(i,j,k) > refmiss) THEN
klow=k
reflow=reflect(i,j,k)
hgtlow=zps(i,j,k)
END IF
END DO
IF(hgtlow < bmhgtmsl) THEN
DO k=2,klow-1
reflect(i,j,k)=reflow
END DO
END IF
END DO
END DO
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
IF(reflect(i,j,k) < -30.) reflect(i,j,k)=-30.
END DO
END DO
END DO
IF(fillopt > 0 ) THEN
varid='rflctx'
CALL wrtvar1(nx,ny,nz,reflect, varid,varname,varunits, &
time,runname,dirname, istatus)
END IF
DO j=1,ny-1
DO i=1,nx-1
DO k=2,nz-1
IF(reflect(i,j,k) > -30.) THEN
refz(i,j,k)=10.**(0.1*reflect(i,j,k))
ELSE
refz(i,j,k)=0.
END IF
cmpref(i,j)=max(cmpref(i,j),reflect(i,j,k))
END DO
END DO
END DO
!
!
!-----------------------------------------------------------------------
!
! 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
pulsewmn=c*pulselmn
ELSE
pulsewmn=pulsewid1
pulselmn=cinv*pulsewmn
END IF
prtmn=prt1
npulsmn=npulse1
ELSE
IF(pulseopt == 1) THEN
pulselmn=0.5*(pulselen1*pulselen2)
pulsewmn=c*pulselmn
ELSE
pulsewmn=0.5*(pulsewid1*pulsewid2)
pulselmn=cinv*pulsewmn
END IF
prtmn=0.5*(prt1+prt2)
npulsmn=NINT(0.5*(npulse1+npulse2))
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)
WRITE(6,'(a,f10.2,a)') ' Mean pulse width: ',pulsewmn,' m'
WRITE(6,'(a,f10.2,a)') &
' Mean pulse length: ',(pulselmn*1.0E06),' microsec'
WRITE(6,'(a,f10.2,a)') &
' Gaussian filtered pulse half-width: ',sqrt(sigr2),' m'
!
!-----------------------------------------------------------------------
! Calculate constants for antenna gain weighting.
!-----------------------------------------------------------------------
!
IF( abs(rotrate) > 0.) THEN
efbmwid=EFFBMW(beamwid,rotrate,prtmn,npulsmn)
ELSE
efbmwid=beamwid
END IF
fourln4=4.0*alog(4.0)
wasqinv=fourln4/(efbmwid*efbmwid)
wesqinv=fourln4/(beamwid*beamwid)
bmwrad=deg2rad*efbmwid
elradius=evalwid*beamwid
hlfgtsp=0.5*gatesp
nptsrng=max(nptsrng,3)
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
nptsazm=max(nptsazm,3)
ALLOCATE (delazm(nptsazm))
nptselv=max(nptselv,3)
ALLOCATE (delelv(nptselv))
de=(evalwid*beamwid)/float(nptselv-1)
WRITE(6,'(a,f10.2,a,f10.2)') &
' Beam width ',beamwid,' Eval region width:',(evalwid*beamwid)
WRITE(6,'(a)') &
' Cross-beam evaluation pts (elev deg from 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
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)') &
' 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
wasqinv=fourln4/(efbmwid*efbmwid)
bmwrad=deg2rad*efbmwid
ALLOCATE (wgtpt(nptsrng,nptsazm,nptselv))
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 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 ( mapproj == 0 ) THEN
IF ( adjctr > 0 ) THEN
ctrlat = ctrlatnew
ctrlon = ctrlonnew
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 = ctrlatnew
ctrlon = ctrlonnew
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
torn0x=radarx+torn_rng*sin(deg2rad*torn_azim)
torn0y=radary+torn_rng*cos(deg2rad*torn_azim)
!
!-----------------------------------------------------------------------
! Adjust grid heights to be compatible with the netrad domain
!-----------------------------------------------------------------------
!
IF(adjhgt > 0) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
zps(i,j,k)=zps(i,j,k)+hgtoffset
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
! Calculate sensitivity values for all distances from radar.
!-----------------------------------------------------------------------
!
CALL xbsens(ngate,gatesp,sens)
!
!-----------------------------------------------------------------------
! Loop through all possible offsets
!-----------------------------------------------------------------------
!
DO joff=joffbgn,joffend
DO ioff=ioffbgn,ioffend
!
!-----------------------------------------------------------------------
!
! Calculate new tornado location for this storm offset.
!
!-----------------------------------------------------------------------
!
printit=((mod(ioff,10) == 0) .AND. (mod(joff,10) == 0))
tornx=torn0x+(ioff*deltax)
torny=torn0y+(joff*deltay)
CALL xytoll(1,1,tornx,torny,tornlat(ioff,joff),tornlon(ioff,joff))
IF (printit) &
WRITE(6,'(a,f12.4,a,f12.4)') ' Tornado lat: ',tornlat(ioff,joff), &
' Tornado lon: ',tornlon(ioff,joff)
!
DO irad=1,nradnet
IF (printit) &
WRITE(6,'(a,a)') ' Processing network radar: ',radnet(irad)
!
! Set up radar location info.
! For now, consider just one radar.
! ioff, joff refer to an offset of the grid and hence storm location
! This is accounted for by moving the radar locations by an equal
! distance in the opposite direction.
!
CALL lltoxy(1,1,radnet_lat(irad),radnet_lon(irad),xrad,yrad)
xrada=xrad-umove*curtim
yrada=yrad-vmove*curtim
xrada=xrada-(ioff*deltax)
yrada=yrada-(joff*deltay)
CALL xytoll(1,1,xrada,yrada,radlata,radlona)
irloc=((xrada-xsc(1))*dxinv)+1.0
jrloc=((yrada-ysc(1))*dyinv)+1.0
IF (printit) &
WRITE(6,'(a,i3,a,i3,/,a,f10.3,a,f10.3,/a,f10.3,a,f10.3)') &
' Ioffset ',ioff,' Joffset ',joff, &
' x-coord (km) ',(0.001*xrada),' y-coord (km) ',(0.001*yrada),&
' i-location: ',irloc,' j-location:',jrloc
!
!-----------------------------------------------------------------------
!
! Calculate radar parameters at each scalar grid point.
!
!-----------------------------------------------------------------------
!
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 beam heading to the position of the tornado
!
delx=torn0x-xrada
dely=torn0y-yrada
IF(abs(dely) .gt. 1.0E-06) THEN
IF(dely .gt. 0.) THEN
azrad=atan(delx/dely)
ELSE
azrad=pi+atan(delx/dely)
END IF
ELSE IF(delx .gt. 0.) THEN
azrad=0.5*pi
ELSE
azrad=1.5*pi
END IF
azimuth=rad2deg*azrad
IF(azimuth < 0.) azimuth=azimuth+360.
IF(azimuth > 360.) azimuth=azimuth-360.
tornsfcr=sqrt((delx*delx)+(dely*dely))
CALL bmhgtsfr(elv,tornsfcr,tornhgt)
CALL beamelv(tornhgt,tornsfcr,tornelv,tornrngi)
IF (printit) THEN
WRITE(6,'(a,3(f10.2))') ' torn delx,dely,azimuth: ', &
(0.001*delx),(0.001*dely),azimuth
WRITE(6,'(a,3(f10.2))') ' torn elv in, elv out: ', &
elv,tornelv
WRITE(6,'(a,3(f10.2))') ' torn sfcrange, range: ', &
(0.001*tornsfcr),(0.001*tornrngi)
END IF
torngate=(tornrngi/gatesp)
IF(torngate < float(ngate)) THEN
!
! For this radial, fill a 1d-array with the reflectivity data
!
kntpt=0
kntclr=0
knthob=0
kntvob=0
kntintr=0
kntfewr=0
DO igate=1, ngate
kntpt=kntpt+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
!
! 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)
! print *, 'bmwm,pulsewid = ',bmwm,pulsewmn
! print *, 'iradius,jradius = ',iradius,jradius
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) THEN
!
! Calculate beam-weighted sums for an uniform array of points around this
! range-gate using tri-linear interpolation
!
kntr=0
sumr=0.
sumwr=0.
DO kpt=1,nptselv
DO jpt=1,nptsazm
DO ipt=1,nptsrng
srgpt=(igate*gatesp)+delrng(ipt)
azmpt=azimuth+delazm(jpt)
elvpt=elv+delelv(kpt)
CALL beamhgt(elvpt,srgpt,hgtpt,sfcrpt)
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/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
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
hgtmsl=hgtpt+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
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
reflz2=whigh*fhigh+wlow*flow
!
! TEST CODE
!
reflzk=refz(iloc1,jloc1,k)
reflzkm1=refz(iloc1,jloc1,k-1)
reflz0=whigh*reflzk+wlow*reflzkm1
!
! TEST CODE 2
!
reflz1=whigh*(w11*refz(iloc,jloc,k) + &
w21*refz(iloc+1,jloc,k) + &
w12*refz(iloc,jloc+1,k) + &
w22*refz(iloc+1,jloc+1,k)) + &
wlow*(w11*refz(iloc,jloc,k-1) + &
w21*refz(iloc+1,jloc,k-1) + &
w12*refz(iloc,jloc+1,k-1) + &
w22*refz(iloc+1,jloc+1,k-1))
WRITE(36,'(5f12.1)') reflzkm1,reflzk,reflz0,reflz1,reflz2
kntr=kntr+1
sumr=sumr+reflz0*wgtpt(ipt,jpt,kpt)
sumwr=sumwr+wgtpt(ipt,jpt,kpt)
END IF ! hgt in grid
END IF ! iloc in grid
END DO
END DO
END DO
!
! 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)=max(refdbz,0.)
ELSE
kntfewr=kntfewr+1
refl(igate)=0.
END IF
ELSE ! no echo found in area
kntclr=kntclr+1
refl(igate)=0.
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
!-----------------------------------------------------------------------
! Calculate attenuation values for all distances from radar.
!-----------------------------------------------------------------------
!
CALL xbatten(ngate,1,refl,gatesp,rmisval,refatt)
!-----------------------------------------------------------------------
! Output R(km), K'(dB), K(dB), Z(dB), Z'(dB), and Ze_min.
!-----------------------------------------------------------------------
IF(ioff == 0 .AND. joff == 0) THEN
WRITE(6,'(5(1x,a10))') &
'R(km)','Z(dB)','Za(dB)','Za-Sens','Sens'
WRITE(6,'(a)') &
'-----------------------------------------------------------------'
DO idx = 5, ngate, 5
WRITE(6,'(5(1x,f10.3))') &
(idx*gatspkm),refl(idx),refatt(idx), &
(refatt(idx)-sens(idx)),sens(idx)
END DO
END IF
!-----------------------------------------------------------------------
! Decide if the tornado is detectable from this radar.
!-----------------------------------------------------------------------
itgate=INT(torngate)
wgt1=torngate-float(itgate)
wgt0=1.-wgt1
tornref=wgt0*refl(itgate)+wgt1*refl(itgate+1)
tornaref=wgt0*refatt(itgate)+wgt1*refatt(itgate+1)
tornsen=wgt0*sens(itgate)+wgt1*sens(itgate+1)
IF(printit) &
WRITE(6,'(a,f10.2,/a,f10.2,a,f10.2,a,f10.2)') &
' At tornado location: Refl=',tornref, &
' Att Ref=',tornaref,' Sens=',tornsen,' Att Ref-Sens=',&
(tornaref-tornsen)
IF(tornaref > tornsen) THEN
score(ioff,joff,irad)=1 ! hit, detectable
ELSE
score(ioff,joff,irad)=0 ! attenuated, undetectable
END IF
ELSE
score(ioff,joff,irad)=-1 ! out of range.undetectable
END IF
END DO ! i-radar
END DO ! i-offset
END DO ! j-offset
!
! Compute summary scores of hits-n-misses
! First summarize by individual radar
!
print *, ' writing tornstat file'
open(42,file=tornstat,status='unknown',form='formatted')
WRITE(42,'(a)') TRIM(runname)
WRITE(42,'(i4.4,a,i2.2,a,i2.2,1x,i2.2,a,i2.2,a,i2.2,a)') &
year,'-',month,'-',day,hour,':',minute,':',second,' UTC'
WRITE(42,'(a,f8.1,a,f8.1,a,f8.1)') &
' Tornado Azim:',torn_azim,' Range:',torn_rngkm,' Refl:',torn_refl
WRITE(42,'(/a)') ' Individual Radar Statistics'
WRITE(42,'(a)') ' Radar Trials Hits Miss Out-of-Range'
kntrad=0
hsum=0.
asum=0.
DO irad=1,nradnet
nt=0
nhits=0
nattn=0
nrng=0
DO ioff=ioffbgn,ioffend
DO joff=joffbgn,joffend
IF(score(ioff,joff,irad) > -90) THEN
nt=nt+1
IF(score(ioff,joff,irad) == 1) THEN
nhits=nhits+1
ELSE IF (score(ioff,joff,irad) == 0) THEN
nattn=nattn+1
ELSE IF (score(ioff,joff,irad) == -1) THEN
nrng=nrng+1
END IF
END IF
END DO
END DO
IF((nt-nrng) > 0 ) THEN
hpct=100.*float(nhits)/float(nt-nrng)
apct=100.*float(nattn)/float(nt-nrng)
hsum=hsum+hpct
asum=asum+apct
kntrad=kntrad+1
ELSE
hpct=0.
apct=0.
END IF
rpct=100.*float(nrng)/float(nt)
WRITE(42,'(a4,i6,3f12.1)') &
radnet(irad),nt,hpct,apct,rpct
IF(kntrad > 0) THEN
hsum=hsum/float(kntrad)
asum=hsum/float(kntrad)
END IF
WRITE(42,'(a4,i6,2f12.1)') &
'AVG:',kntrad,hsum,asum
END DO
!
! Compute hits and misses for the network as a whole
!
WRITE(6,'(a)') ' Computing hits and misses'
nt=0
nhits=0
nattn=0
nrng=0
DO joff=joffbgn,joffend
DO ioff=ioffbgn,ioffend
nt=nt+1
hit=.false.
rng=.false.
DO irad=1,nradnet
IF(score(ioff,joff,irad) == 1) hit=.true.
IF(score(ioff,joff,irad) /= -1) rng=.true.
END DO
IF(hit) THEN
nhits=nhits+1
netscore(ioff,joff)=1 ! hit, detectable
ELSE IF (.not.rng) THEN
nrng=nrng+1
netscore(ioff,joff)=-1 ! out-of-range
ELSE
nattn=nattn+1
netscore(ioff,joff)=0 ! network miss
END IF
END DO
END DO
IF((nt-nrng) > 0) THEN
hpct=100.*float(nhits)/float(nt-nrng)
apct=100.*float(nattn)/float(nt-nrng)
ELSE
hpct=0.
apct=0.
END IF
rpct=100.*float(nrng)/float(nt)
print *, ' writing tornfile file'
open(41,file=tornfile,status='unknown',form='formatted')
WRITE(41,'(a)') TRIM(runname)
WRITE(41,'(i4.4,a,i2.2,a,i2.2,1x,i2.2,a,i2.2,a,i2.2,a)') &
year,'-',month,'-',day,hour,':',minute,':',second,' UTC'
WRITE(41,'(a,f8.1,a,f8.1,a,f8.1)') &
' Tornado Azim:',torn_azim,' Range:',torn_rngkm,' Refl:',torn_refl
DO joff=joffbgn,joffend
DO ioff=ioffbgn,ioffend
WRITE(41,'(2i6,2f12.4,13i5)') &
ioff,joff,tornlat(ioff,joff),tornlon(ioff,joff), &
netscore(ioff,joff), &
(score(ioff,joff,irad),irad=1,nradnet)
END DO
END DO
CLOSE(41)
WRITE(42,'(//a)') 'Network Statistics'
WRITE(42,'(a)') 'Trials Hits Miss Out-of-Range'
WRITE(42,'(i6,3f12.1)') nt,hpct,apct,rpct
CLOSE(42)
CLOSE(9)
STOP
END PROGRAM attenuation
!-----------------------------------------------------------------------
! Subroutine atten
!-----------------------------------------------------------------------
!
SUBROUTINE atten(all_rngs,all_refs,seg,att_values,att)
IMPLICIT NONE
INTEGER :: idx ! Index.
INTEGER :: seg ! Number of distance segments.
REAL :: att_tally ! Running tally of attenuation values.
REAL :: att_const1 ! Attenuation constant 1.
REAL :: att_const2 ! Attenuation constant 2.
REAL :: att_const3 ! Attenuation constant 3.
REAL, DIMENSION(seg),INTENT (IN):: all_rngs ! Range values for all dist from radar.
REAL, DIMENSION(seg),INTENT (IN):: all_refs ! Reflectivity values for all dist from radar.
REAL, DIMENSION(seg),INTENT (OUT) :: att_values ! Attenuation values for all dist from radar.
REAL :: att(seg) ! Attenuation value for a specific distance.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
att_const1 = 2.0*0.01
att_const2 = 1.0/400.0
att_const3 = 1.21/1.4
att_tally = 0
DO idx = 1, seg
att(idx) = &
att_const1*((att_const2*(10.0**(0.1*all_refs(idx))))**att_const3)
IF(idx == 1) THEN
att_tally = att_tally + att(idx)
att_values(idx) = att_tally
ELSE IF(all_refs(idx) == all_refs(idx-1)) THEN
att_tally = att_tally + att(idx)
att_values(idx) = att_tally
ELSE
att_tally = att_tally + .75*att(idx-1) + &
.25*att_const1*((att_const2*(10.0**(0.1*all_refs(idx))))**att_const3)
att_values(idx) = att_tally
END IF
END DO
END SUBROUTINE atten
!-----------------------------------------------------------------------
! Subroutine sens
!-----------------------------------------------------------------------
!
!
SUBROUTINE xbsens(ngate,gatesp,sens) 2
!
!-----------------------------------------------------------------------
! Subroutine xbsens
! Calculate X-band radar sensitivity (dBZ) as a function of range.
!
! Based on information provided by Francesc Joyvount
!
! Keith Brewster and Erin Fay, CAPS
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: ngate
REAL, INTENT(IN) :: gatesp
REAL, INTENT (OUT):: sens(ngate) ! Sensitivity values for all dist from radar.
!
! Parameters for CASA radar
!
REAL, PARAMETER :: lmbda = 3.0E-02 ! Wavelength (m)
REAL, PARAMETER :: pt = 12500 ! Peak power (Watts)
REAL, PARAMETER :: G = 38.0 ! Antenna gain (dB)
REAL, PARAMETER :: F = 5.5 ! Noise figure (dB)
REAL, PARAMETER :: tau = 0.67E-06 ! Radar pulse length (s)
REAL, PARAMETER :: theta = 2.0 ! Antenna half-power beamwidth (deg)
REAL, PARAMETER :: B = 2.5 ! Bandwidth (MHz)
REAL, PARAMETER :: lm = 1.0 ! Receiver mis-match loss(dB)
REAL, PARAMETER :: Kw2 = 0.91 ! Refractive Index of water squarred
REAL, PARAMETER :: T0 = 300. ! Rx temperature (K)
REAL, PARAMETER :: Ta = 200. ! Antenna temperature (K)
!
! Physical constants
!
REAL, PARAMETER :: K = 1.38E-23 ! Boltsmann's Constant (J/K)
REAL, PARAMETER :: c = 2.99792E08 ! Speed of light (m/s)
REAL, PARAMETER :: rconst =1.0E18 ! m6/m3 to mm6/m3
!
! Misc internal variables
!
INTEGER :: igate
REAL :: ln2,pi,pi3,four3,bwrad,rnoise,BHz,Ni,sconst,rlm,rG
REAL :: range
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ln2=log(2.0)
pi=acos(-1.)
pi3=pi*pi*pi
four3=4.*4.*4.
bwrad=theta*pi/180.
rnoise=10.**(0.1*F)
rlm=10.**(0.1*lm)
rG=10.**(0.1*G)
BHz=1.0E06*b
Ni=K*(Ta+(rnoise-1.0)*T0)*BHz
sconst=rconst *(Ni/(pi3*Kw2)) * ((8.*ln2)/(bwrad*bwrad)) * &
(2./(c*tau)) * ((lmbda*lmbda*four3*rlm)/(Pt*rG*rG))
DO igate = 1, ngate
range=igate*gatesp
sens(igate) = 10.*LOG10(range*range*sconst)
END DO
END SUBROUTINE xbsens
SUBROUTINE xbatten(ngate,nazim,refl,gatesp,rmisval,refatt) 2
!
!-----------------------------------------------------------------------
!
! 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
!
!
! OUTPUT:
! refatt: Attenuated reflectivity.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: ngate ! Number of gates
INTEGER, INTENT (IN) :: nazim ! Number of radials
REAL, INTENT (IN):: refl(ngate,nazim) ! Reflectivity (dBZ)
REAL, INTENT (IN) :: gatesp ! Gate spacing (m)
REAL, INTENT (IN) :: rmisval
REAL, INTENT (OUT) :: refatt(ngate,nazim) ! Attenuated reflectivities
INTEGER :: igate,jazim ! Index.
REAL :: att_tally ! Running tally of attenuation values.
REAL :: att_const ! Attenuation constant
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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
! Use Z-R relationship to develop attenuation constant for rain.
! refl(dBZ)=10 log10(zrconst*R**zrexp)
! So R=(10**(refl*.1))**(1/zrexp)
! Kr(dbZ/km)=krconst*(R**krexp)
! So
! Kr=(krconst/zrconst)*(10**(refl*0.1)**(krexp/zrexp))
!
! Calculating two-way attenuation requires this figure be doubled.
!
att_exp=krexp/zrexp
att_const=2.*krconst/(zrconst**att_exp)
hlfgatsp=0.5*gatesp*mperkm
DO jazim=1, nazim
IF(refl(1,jazim) .gt. rmisval) THEN
atten = att_const*((10.0**(0.1*refl(1,jazim)))**att_exp)
att_tally = atten*hlfgatsp
refatt(1,jazim) = refl(1,jazim)-att_tally
ELSE
atten=0.
att_tally=0.
refatt(1,jazim)=refl(1,jazim)
END IF
atten_last=atten
DO igate = 2, ngate
IF(refl(igate,jazim) .gt. rmisval) THEN
atten = att_const*((10.0**(0.1*refl(igate,jazim)))**att_exp)
att_tally = att_tally + hlfgatsp*atten_last + hlfgatsp*atten
refatt(igate,jazim) = refl(igate,jazim)-att_tally
atten_last=atten
ELSE
refatt(igate,jazim) = refl(igate,jazim)
atten_last=0.
END IF
END DO
END DO
RETURN
END SUBROUTINE xbatten
!
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