pe.f90

These routines model tropospheric radiowave propagation over variable terrain and calculates propaga

! ************* THIS FILE CONTAINS TPEM MODEL SUBROUTINES *******************
! Author: Amalia E. Barrios
! NCCOSC RDT&E DIV D883
! 49170 Propagation Path
! San Diego, CA 92152-7385
! e-mail: [e-mail]
! phone: [phone]
! fax: [phone]
! Summary: These routines model tropospheric radiowave propagation over
! variable terrain and calculates propagation loss vs. height and
! range. Propagation loss is displayed in dB contours on a height vs.
! range plot. TPEM is based on the split-step Fourier PE method and
! was originally developed from an early PE model called PEPC, written
! by Fred Tappert. Propagation loss over variable terrain is modeled
! by shifting the field an appropriate number of bin widths correspondc
! ing to the height of the ground. The field is determined using the
! smooth earth PE method.

! *********************************************************************************
! Variables in small letters in parameter lists are variables that are input
! or passed to called subroutines. Variables in CAPS in parameter lists are
! returned from the called subroutines.
! ********************************************************************************

! Main Glossary of Common Blocks used in all subroutines:
! ARRAYS:
!    ENVPR() = Complex array containing refractivity exponential term.
!              i.e. ENVPR() = exp[i * dr * k * le-6 * M(z)
!    FILT() = Cosine-tapered (Tukey) filter array.
!    FRSP() = Complex array containing free-space propagator exponential term.
!             i.e., FRSP() = exp[-i * dr * (k - sqrt(k**2 - p**2)) ].
!    U() = Complex array containing field solution at current PE range.
!    ULST() = Complex array containg field solution at previous PE range.
! HTVAR:
!    YLAST = Height of ground in meters at the last PE range.
!    YCUR = Height of ground in meters at the current PE range.
!    YCURM = Height of ground in meters midway between last and current PE range
!            step. For use when shifting profiles to be relative to the local
!            ground height.
! IMPEDANCE:
!    ALPHAV = Vertical polarization impedance term = i*fko/rng.
!    C1 = Coefficient used in vertical polarization calculations.
!    C2 = Coefficient used in vertical polarization calculations.
!    CIM = Constant for each calculated ALPHAV - used in C1 calculation.
!    C2M = Constant for each calculated ALPHAV - used in C2 calculation.
!    IG = Counter indicating current ground type being modeled.
!    RAV() = Array of ROOT to the i'th power, i.e. RAV(I) = ROOT**I.
!    RK = Coefficient used in Cl and C2 calculations.
!    RNG = Complex refractive index.
!    RNG2 = Complex refractive index squared.
!    ROOT = Complex root of quadratic equation for mixed transform method
!    based on Kuttler's formulation.
! MISCVAR:
!    ANTREF = Transmitting antenna height relative to the reference
!             height HMINTER.

!    CNST = Used in calculating FRSP() in routine PHASE1.
!    CNST = DELP/FKO.
!    DELP = Mesh size in angle- (or p-) space.
!    FNORM = Normalization factor used for DFT.
!    FTER = Logical flag - .TRUE.=terrain case, .FALSE.=smooth surface case
!    HLIM() = Array containing height at each output range at which the
!             last valid loss value exists.
!    HMREF = Height relative to HMINTER. Determined from user-provided
!            minimum height VNP.HMIN. That is, if VNP.HMIN is minimum
!            height input by user with respect to mean sea level,
!            and HMINTER is internally considered the new origin,
!            then HMREF = VNP.HMIN - HMINTER.
!    HTLIM = Maximum desired calculation height with respect to HMINTER,
!            i.e., HTLIM = VNP.HMAX-HMINTER.
!    PLCNST = Constant used in determining propagation loss.
!    PLCNST = 20log(2*FKO).
!    QI = Imaginary i -> complex(O,l).
!    RPE = Range at which valid loss values will begin to be calculated.
!    THETAMAX = Maximum propagation angle in PE calculations.
!    SLP() = Slope of each segment of terrain.


! PARINIT:
!    HT() = PE mesh height array of size N.
!    HTDUM() = Dummy array containing height values for current (horizontally
!              interpolated) profile.
!    IS = Counter for current profile (for range-dependent cases).
!    LVLEP = Number of height/refractivity levels in profile.
!    PROFINT() = M-unit profile interpolated to every DELZ in height.
!    REFDUM() = Dummy array containing M-unit values for current (horizontally
!               interpolated) profile.
!    RV2 = Range of the next refractivity profile (for range-dependent cases).

! PATTERN:
!    AFAC = Constant used in determining antenna pattern factors.
!    AFAC = 1.39157 / sin( bw / 2 ) for SIN(X)/X and height-finder.
!    AFAC = (.5*in(2))/(sin(bw/2))**2 for GAUSSIAN.
!    BW = Antenna pattern beamwidth in radians.
!    ELV = Antenna pattern elevation angle in radians.
!    PELEV = Sine of elevation angle.
!    SBW = Sine of the beamwidth.
!    UMAX = Limiting angle used in 30 dB cut-off point for SIN(X)/X and
!           generic height-finder antenna pattern factors.

! PEVAR:
!    CON = l.e-6 * FKO; Constant used in calculation of ENVPR(.
!    DELZ = Bin width in z-space = WL / (2*sin(THETAMAX)).
!    DZ2 = 2. * DELZ.
!    FKO = Free-space wavenumber = (2*pi) / WL
!    LN = Power of 2 transform size, i.e. N = 2**LN.
!    N = Transform size.
!    N34 = 3/4 * N.
!    NMl = N-1.
!    WL = Wavelength in meters.
!    ZMAX = Maximum height of PE calculation domain = N * DELZ.

! PROFWREF:
!    HREF() = Heights of refractivity profile with respect to local ground
!             height.
!    NLVL = Number of levels in new profile.
!    REFREF() = Corresponding refractivity array for HREF().

! RHSTPS:
!    DR = PE range step in meters.
!    DR2 = 1/2 PE range step in meters.
!    DROUT = Output range step in meters.
!    DZOUT = Output height increment in meters.
!    ZOUT() = Array containing all output height points.

! TRVAR:
!    DMDH() = Gradients of first profile in M-units/meters.
!    JLS = Index of refractivity array at which antenna height is located.
!    RLIM = 90% of maximum range RMAX - used for ray tracing.
!    THETALAUNCH = Angle in radians of launch angle for which, when traced,
! height of the ray at each output range step is stored.
! ZLIM = Height limit for ray trace.


! * SUBROUTINE PEINIT ******************************
! Module Name: PEINIT
! Module Security Classification: UNCLASSIFIED
! Purpose: Initializes all variables used in TPEM subroutines for PE calculactions.
! Version Number: 1.5
! INPUTS:
! Argument list: EF(errorflag) structure, RF(refractivity) structure,
! SV(systemvar) structure, TR(terrain) structure,
! VNP(inputvar) structure
! OUTPUTS:
! Argument list: HMINTER, IERROR, ROUT
! Common: AFAC, ALPHAV, ANTREF, BW, Cl, C2, CIM, C2M, CNST, CON,
! DELP, DELZ, DMDH(), DR, DR2, DROUT, DZ2, DZOUT, ELV, ENVPRO,
! FILT(, FKO, FNORM, FRSP(, FTER, HLIM(, HMREF, HT(, HTDUMO,
! HTLIM, IG, IS, JLS, LN, LVLEP, N, N34, NM1, PELEV, PLCNST,
! PROFINT(), QI, RAVO, REFDUM(), RPE, RK, RLIM, RNG, RNG2, ROOT,
! RV2, SBW, SLP(, THETALAUNCH, THETAMAX, U(), UMAX, WL, YCUR,
! YCURM, YLAST, ZLIM, ZMAX, ZOUT()

! FILES INCLUDED: FFTSIZ.INC, TPEM.INC
! CALLING ROUTINES: MAIN DRIVER PROGRAM or TESS CSCI
! ROUTINES CALLED: DIEINIT, FFT, GETFFTSZ, INTPROF, PHASE1, PHASE2, PROFREF,
! REFINIT, TRACEA, TRACEH, XYINIT
! GLOSSARY:
! EF = Error flag structure for external implementation constants.
! EF.LERR6 = Logical flag that allows for greater flexibility in
! allowing error -6 to be bypassed. If set to .TRUE.
! then trapping for this erroroccurs, otherwise it can
! be totally ignored by main driver program. (Within the
! TPEM program it is handled as a warning). If this
! error is bypassed (EF.LERR6 = .FALSE.) terrain profile
! is extended to RMAX with same elevation height of last
! valid terrain profile point.
! EF.LERR12 = Same as EF.LERR6 - allows for trapping of this error.
! If LERR12 = .FALSE., then (for range-dependent case)
! if range of last refractivity profile entered is less
! than RMAX, the environment is treated as homogeneous
! from the last profile entered to RMAX.
! RF = Refractivity structure for external environmental data elements.
! RF.LVLEP = Number of levels in refractivity profile (for range
! dependent case all profiles must have same number of
! levels).
! RF.REFMSL(,) = 2-dimensional array containing refractivity with
! respect to mean sea level of each profile. Array
! format must be REFMSL(I,J) = M-unit value at Ith
! level of Jth profile. J = 1 for range-independent
! cases.
! RF.HMSL(,) = 2-dimensional array containing heights in meters with
! respect to mean sea level of each profile. Array format
! must be HMSL(I,J) = height of Ith level of Jth profile.
! J = 1 for range-independent cases.
! RF.RNGPROF() = Ranges of each profile in meters, i.e., RNGPROF(I) =

! range of Ith profile. RNGPROF(1) should always be
! equal to 0.
! RF.NPROF = Number of profiles. Equals 1 for range-independent cases.
! SV = System structure for external system data elements.
! SV.FREQ = Frequency in MHz.
! SV.ANTHT Transmitting antenna height above local ground in meters.
! SV.BWIDTH = Half-power (3 dB) antenna pattern beamwidth in degrees
! (.5 to 45.).
! SV.ELEV = Antenna pattern elevation angle in degrees. (-10 to 10).
! SV.POLAR = 1-character string indicating polarization. H-horizontal,
! V-vertical
! SV.IPAT = Integer value indicating type of antenna pattern desired.
! IPAT = 0 ->omni
! IPAT = 1 -> gaussian
! IPAT = 2 ->sinc x
! IPAT = 3 >csc**2 x
! IPAT = 4 -> generic height-finder
! TR = Terrain structure for external terrain data elements.
! TR.TERX() = Range points of terrain profile in meters.
! TR.TERY() = Reight points of terrain profile in meters.
! TR.ITP = Number of height/range pairs in profile.
! TR.IGR = Number of different ground types specified.
! TR.IGRND() = Type of ground composition for given terrain profile -
! can vary with range. Different ground types are:
! 0 = sea water, 1 = fresh water, 2 = wet ground,
! 3 = medium dry ground, 4 = very dry ground,
! 5 = user defined (in which case, values of relative
! permittivity and conductivity must be given).
! TR.RGRND() = Ranges at which the ground types apply.
! TR.DIELEC(,) = 2-dimensional array containing the relative
! permittivity and conductivity; DIELEC(l,i) and
! DIELEC(2,i), respectively. Only needs to be specified
! if using IGRND(i) = 5, otherwise, TPEM will
! calculate based on frequency and ground types 0-4.
! VNP = Inputvar structure for external implementation constants.
! VNP.HMAX = Maximum output height with respect to m.s.l. in meters.
! VNP.HMIN = Minimum output height with respect to m.s.l. in meters.
! VNP.RMAX = Maximum output range in meters.
! VNP.NZOUT = Integer number of output height points desired.
! VNP.NROUT = Integer number of output range points desired.
! VNP.PROPANG = Maximum problem (propagation) angle in degrees
! desired for solution. If set to 0., then TPEM will
! determine it's own.
! HMINTER = Height of the minimum elevation of terrain profile. This
! will be used to adjust entire terrain profile so subsequent
! loss values returned will be referenced to this height.
! IERROR = Integer value that is returned if any errors exist in input data:
! -6 : Last range in terrain profile is less than VNP.RMAX. (Will
! only return this error if error flag EF.LERR6 is set to
! .TRUE.).
! -8 : VNP.HMAX is less than maximum height of terrain profile.
! -12 : Range of last refractivity profile entered (for range depenc
! dent case) is less than RMAX. (This is returned from subrouc
! tine REFINIT). Will only return this error if error flag
! EF.LERR12 is set to .TRUE.).
! -14 : Last gradient in any refractivity profile entered is
! negative.
! (This is returned from REFINIT).
! -17 : Range points of terrain profile is not increasing.
! -18 : First range point is not 0.
! -42 : Minimum height input by user (VNP.HMIN) is greater then
! maximum height (VNP.HMAX).
! ROUT = Output range point (meters) - initialized in this routine
! For output variables in common blocks, look in main glossary for
! associated common blocks for variable definitions.

! Common blocks: variables
! ARRAYS: ENVPR(), FILT(, FRSP(, U()
! PEVAR: CON, DELZ, DZ2, FKO, LN, N, N34, NMl, WL, ZMAX,
! RHSTPS: DR, DR2, DROUT, DZOUT, ZOUT()
! MISCVAR: ANTREF, CNST, DELP, FNORM, FTER, HLIM(, HMREF, HTLIM,
! PLCNST, QI, RPE, SLP(, THETAMAX
! PATTERN: AFAC, BW, ELV, PELEV, SBW, UMAX
! PARINIT: HT(), HTDUM(), IS, LVLEP, PROFINT(), REFDUM(), RV2
! HTVAR: YCUR, YCURM, YLAST
! TRVAR: DMDH(), JLS, RLIM, THETALAUNCH, ZLIM
! IMPEDANCE: ALPHAV, Cl, C2, CIM, C2M, IG, RAVO, RK, RNG, RNG2, ROOT
! Local Variables:
! ACRIT = Critical angle, i.e., angle above which no rays are trapped
! for ducting environment.
! ANGLE = Tangent ray angle from source to each terrain peak along
! profile.
! ANGU = Maximum tangent ray angle from souce to terrain peak along
! terrain profile path
! BUGFIX = Dummy variable used to push/move bits in EOF statement.
! This is a known bug in MS Fortran Powerstation 1.0 and is
! not necessary for other flavors of Fortran.
! CO = Speed of light in m/s
! CIC = Intermediate complex number used to calculate Cl
! C2C = Intermediate complex number used to calculate C2
! CN75 =4 * pi / N
! CX = RAV(I) if I is even or -RAV(I) if I is odd
! DXI = Difference between current and previous "unprocessed" range
! points (i.e.,XI-XMI)
! DX2 = Difference between current and next "unprocessed" range points
! (i.e., XPl-XI)
! GRAD = Gradient of current refractivity layer in profile
! Hl = (I)th height value in first refractivity profile
! H2 = (I+l)th height value in first refractivity profile
! HDEG = 1/2 degree in radians
! HTERMAX = Maximum height in meters along terrain profile.
! ISCR = Unit number for scratch file - used for temporary storage of
! terrain profile points for processing.
! LFLAG = Logical flag indicating if end-of-file has been reached in the
! scratch file for reading back processed terrain points.
! LOPEN = Logical flag returned in INQUIRE statement that checks if a
! file is already attached to unit ISCR.
! N04 = N/4
! PRANG = User-specified propagation angle in radians.
! RADC = PI/180 -> used for converting angles given in degrees to
! radians.
! RHOR = Radar horizon range based on transmitter height ANTREF and
! 0 receiver height.
! RKM = Maximum range in km.
! RLLIM = Various maximum range step limits based on geometry and
! whether terrain or smooth surface case is being performed
! RMl = (I)th M-unit value in first refractivity profile
! RM2 = (I+l)th M-unit value in first refractivity profile
! RMATHT = M-unit value at transmitter height ANTREF
! RMINA = Minimum M-unit value above transmitter height
! RMINB = Minimum M-unit value for all heights below transmitter height
! SCD = Difference in previous and next "unprocessed" slope segments
! (i.e., SCD = second derivative, d^2y/dx^2
! SDEG10 = Sine of 10 degrees.
! SDEG15 = Sine of 15 degrees.
! SLl = Slope of previous "unprocessed" terrain segment
! SL2 = Slope of next "unprocessed" terrain segment
! SLOPE = Slope of current terrain segment (after processing)
! STHETAMAX = Sine of THETAMAX
! THETAFRAC = Fractional number relating THETALAUNCH to THETAMAX
! UI = (I)th value of complex PE field
! UNMI = (N-I)th value of complex PE field

! Xl = Range of Ith terrain point in meters (after processing)
! X2 = Range of (I+l)th terrain point in meters (after processing)
! XI = Range of "unprocessed" (I)th terrain point in meters
! XDIF = Range difference between current and next terrain points after
! processing (i.e., X2-Xl)
! XMI = Range of "unprocessed" (I-l)th terrain point in meters
! XPI = Range of "unprocessed" (I+l)th terrain point in meters
! Yl = Height of Ith terrain point in meters (after processing)
! Y2 = Height of (I+l)th terrain point in meters (after processing)
! YDIF = Height difference between current and next terrain points after
! processing (i.e., Y2-Y1)
! YI = Height of "unprocessed" (I)th terrain point in meters
! YMI = Height of "unprocessed" (I-l)th terrain point in meters
! YPl = Height of "unprocessed" (I+l)th terrain point in meters

subroutine peinit( ef, vnp, rf, sv, tr, HMINTER, ROUT, IERROR)
include 'tpem.inc'
common / arrays / u(0:maxpts), filt(O:maxn4), frsp(0:maxpts), envpr(0:maxpts), ulst(0:maxpts)
common / pevar / wl, fko, delz, n, in, zmax, n34, con, dz2, nml
common / rhstps / dr, drout, dzout, dr2, zout(mxzout)
common / miscvar / fnorm, cnst, delp, thetamax, plcnst, qi, antref, rpe, hlim(mxrout), slp(mxter), fter, hmref, htlim
common / pattern / pelev, afac, bw, elv, umax, sbw
common / parinit / rv2, refdum(mxlvls), htdum(mxlvls), profint(O:maxpts), ht(0:maxpts), is, ivlep
common / htvar / ylast, ycur, ycurm
common / trvar / dmdh(mxlvls), zlim, jls, thetalaunch, rlim
common / impedance / alphav, rav(O:maxpts), rng, rng2, cl, c2, rk, clm, c2m, ig, root
record / errorflag / ef
record / inputvar / vnp
record / refractivity / rf
record / systemvar / sv
record / terrain / tr
complex u, frsp, envpr, ulst, qi, alphav, ray, rng, rng2, clc
complex c2c, c2m, rk, chm, cl, c2, root, ui, unmi, cx
logical fter, lopen, iflag
data radc / 1.74533e-2 / !degree to radian conversion factor
data iscr / 20 / ! Unit number for scratch file
data cO / 299.79245 / !speed of light x le.-6 m/s
data sdeglO / .173648177 / ! Sine of 10 degrees
data sdegl5 / .258819045 / ! Sine of 15 degrees
data hdeg / 8.726646e-3 / ! 1/2 degree
ierror = 0
fter = .false.
thetamax = 0.
hminter = 0.
angu = 0.
antref = sv.antht
! Put lower limit on HMAX and RMAX
vnp.rmax = amaxl( vnp.rmax, 5000. ) !Set max. range to no less than 5 km.
vnp.hmax = amaxl( vnp.hmax, 100. ) !Set max. height to no less than 100 m.
if( vnp.hmin .ge. vnp.hmax ) then
ierror = -42
return
end if

vnp.hmin = aminl( vnp.hmin, vnp.hmax-100. )
dzout = (vnp.hmax-vnp.hmin) / float( vnp.nzout)
drout = vnp.rmax / float( vnp.nrout)
! Setup output height array
do i = 1, vnp.nzout
zout(i) = vnp.hmin + float(i) * dzout
end do
WL = cO / sv.freq
FKo = 2. * pi / WL
con = l.e-6 * fko
qi = cmplx( 0., 1.)
! Calculate constants used to determine antenna pattern factor
! SV.IPAT = 0 -> omni
! SV.IPAT = 1 -> gaussian
! SV.IPAT = 2 -> sinc x
! SV.IPAT = 3 - csc**2 x
! SV.IPAT = 4 -> generic height-finder
sv.bwidth = amaxl( sv.bwidth, .5 ) !Don't let vertical beamwidth go
sv.bwidth = aminl( sv.bwidth, 45. ) !outside .5 to 45 degree limit.
sv.elev = amaxl( sv.elev, -10. ) !Don't let elevation angle go
sv.elev = aminl( sv.elev, 10. ) !outside -10 to 10 degree limit.
bw = sv.bwidth * radc
elv= sv.elev * radc
bw2 = .5 * bw
if( sv.ipat .eq. 1 ) then
afac = .34657359 / (sin( bw2 ))**2
pelev = sin( elv )
elseif( sv.ipat .eq. 3 ) then
sbw = sin( bw )
elseif( sv.ipat .ne. 0 ) then
afac = 1.39157 / sin( bw2)
a = pi / afac
umax = atan( a / sqrt(l. - a**2))
end if
! Discard any unnecessary terrain points. Test on the rate of change of slope,
! i.e., second derivative. If that is below l.e-3 then discard that point.
if( tr.itp .gt. 0 ) fter = .true.
if( fter ) then
! Check that all terrain range points are increasing.
do i = 1, tr.itp-I
ipl = i + 1
if( tr.terx(ipl) .lt. tr.terx(i) ) then ierror = -17
return
end if
end do
! Test to see that first range value is 0.
if( tr.terx(l) .gt. 0. ) then
ierror = -18
return
end if
! Test to see if the last range point in the profile meets or exceeds PMAX. If

! not then return error code.
if( tr.terx(tr.itp) .lt. vnp.rmax ) then