1module radialVelocity_mod
2 ! MODULE radialVelocity_mod (prefix='rdv' category='5. Observation operators')
3 !
4 !:Purpose: Containing commonly used functions for the assimilation of Doppler velocity.
5 !
6 !
7 use codePrecision_mod
8 use earthConstants_mod
9 implicit none
10 save
11 private
12
13 ! public procedures
14 public :: rdv_getlatlonHRfromRange, rdv_getRangefromH
15
16
17contains
18
19
20 subroutine rdv_getlatlonHRfromRange(antennaLat, antennaLon, beamElevation, beamAzimuth, & !in
21 radarAltitude, beamRange, & !in
22 latSlant, lonSlant, beamHeight, beamDistance) !out
23 !
24 !:Purpose: Computation of lat-lon , height of the trajectory
25 ! along the radar beam from range of the radar beam
26 !
27 ! Angles are expressed in radians
28 ! Height and distances in meters
29 !
30 implicit none
31
32 ! Arguments:
33 real(pre_obsReal), intent(in) :: antennaLat
34 real(pre_obsReal), intent(in) :: antennaLon
35 real(pre_obsReal), intent(in) :: beamElevation
36 real(pre_obsReal), intent(in) :: beamAzimuth
37 real(pre_obsReal), intent(in) :: radarAltitude
38 real(pre_obsReal), intent(in) :: beamRange
39 real(pre_obsReal), intent(out) :: LatSlant
40 real(pre_obsReal), intent(out) :: lonSlant
41 real(pre_obsReal), intent(out) :: beamHeight
42 real(pre_obsReal), intent(out) :: beamDistance
43
44 ! Locals:
45 real(pre_obsReal) :: Re
46
47 ! Radius of sphere of equal area from earthconstants_mod.f90
48 ! ec_wgs_R2 = 6371007.1809
49 ! effective radius of the earth
50 Re = ec_wgs_R2 * (4./3.)
51
52 !compute height of radar observation
53 beamHeight = sqrt(beamRange**2.+(Re+radarAltitude)**2.+2.*beamRange*(Re+radarAltitude)*sin(beamElevation))-(Re)
54
55 ! distance following surface of the earth from Doviak and Zrnic (2.28c)
56 beamDistance = atan(beamRange*cos(beamElevation)/(beamRange*sin(beamElevation)+Re+radarAltitude))*Re
57
58 ! lat lon of the path along the radar beam
59 latSlant = asin( sin(antennaLat)*cos(beamDistance/ec_wgs_R2) + cos(antennaLat)*sin(beamDistance/ec_wgs_R2)*cos(beamAzimuth))
60 lonSlant = antennaLon + atan2(sin(beamAzimuth)*sin(beamDistance/ec_wgs_R2)*cos(antennaLat), cos(beamDistance/ec_wgs_R2)-sin(antennaLat)*sin(latSlant))
61
62 end subroutine rdv_getlatlonHRfromRange
63
64 subroutine rdv_getRangefromH(beamHeight, radarAltitude, beamElevation, beamRange)
65 !
66 !:Purpose: Computation of range of the radar beam from height of the radar beam
67 !
68 implicit none
69
70 ! Arguments:
71 real(pre_obsReal) , intent(in) :: beamHeight
72 real(pre_obsReal) , intent(in) :: radarAltitude
73 real(pre_obsReal) , intent(in) :: beamElevation
74 real(pre_obsReal) , intent(out) :: beamRange
75
76 ! Locals:
77 real(pre_obsReal) :: a, b, c, Re
78
79 if ( radarAltitude > beamHeight ) then
80 !beamHeight is below radar antenna which may cause the equation below to return garbage
81 !this happens in a few edge cases where its okay to return zero
82 beamRange = 0.0
83 else
84
85 ! Radius of sphere of equal area from earthconstants_mod.f90
86 ! ec_wgs_R2 = 6371007.1809
87 ! effective radius of the earth
88 Re = ec_wgs_R2*(4./3.)
89
90 a = 1.
91 b = 2.*(Re + radarAltitude)*sin(beamElevation)
92 c = -(Re + beamHeight)**2. + (Re + radarAltitude)**2.
93 ! range of radar beam from height and elevation of the radar beam
94 beamRange = (-b + sqrt( b**2. - 4.*a*c )) / (2.*a)
95
96 end if
97
98 end subroutine rdv_getRangefromH
99
100end module radialVelocity_mod
