UFO
ufo_gnssro_bendmetoffice_utils_mod.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! (C) British Crown Copyright 2020 Met Office
3 !
4 ! This software is licensed under the terms of the Apache Licence Version 2.0
5 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 !-------------------------------------------------------------------------------
7 module ufo_gnssro_ukmo1d_utils_mod
8 
9 !use iso_c_binding
10 use fckit_log_module, only: fckit_log
11 use kinds, only: kind_real
12 
13 ! Generic routines from elsewhere in jedi
14 use missing_values_mod
15 use ufo_constants_mod, only: &
16  rd, & ! Gas constant for dry air
17  cp, & ! Heat capacity at constant pressure for air
18  rd_over_cp, & ! Ratio of gas constant to heat capacity
19  pref, & ! Reference pressure for calculating exner
20  pi, & ! Something to do with circles...
21  grav, & ! Gravitational field strength
22  ecc, & ! eccentricity
23  k_somig, & ! Somigliana's constant
24  g_equat, & ! equatorial gravity (ms-2)
25  a_earth, & ! semi-major axis of earth (m)
26  flatt, & ! flattening
27  m_ratio, & ! gravity ratio
28  mw_ratio, & ! Ratio of molecular weights of water and dry air
29  c_virtual, & ! Related to mw_ratio
30  n_alpha, & ! Refractivity constant a
31  n_beta ! Refractivity constant b
32 
33 implicit none
34 public :: ops_gpsrocalc_alpha
35 public :: ops_gpsrocalc_nr
36 public :: ops_gpsro_geop_geom
37 private
38 
39 contains
40 
41 !-------------------------------------------------------------------------------
42 ! GPSRO 1D bending angle operator.
43 !-------------------------------------------------------------------------------
44 SUBROUTINE ops_gpsrocalc_alpha (nobs, &
45  nlev, &
46  a, &
47  refrac, &
48  nr, &
49  alpha)
50 
51 IMPLICIT NONE
52 
53 ! Subroutine arguments:
54 INTEGER, INTENT(IN) :: nobs ! size of ob. vector
55 INTEGER, INTENT(IN) :: nlev ! no. of refractivity levels
56 REAL(kind_real), INTENT(IN) :: a(nobs) ! observation impact parameters
57 REAL(kind_real), INTENT(IN) :: refrac(nlev) ! refractivity values on model levels
58 REAL(kind_real), INTENT(IN) :: nr(nlev) ! refractive index * radius product
59 REAL(kind_real), INTENT(OUT) :: alpha(nobs) ! bending angle
60 
61 ! Local declarations:
62 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_alpha"
63 INTEGER :: i
64 INTEGER :: n
65 INTEGER :: ibot
66 INTEGER :: jbot
67 INTEGER :: kbot
68 REAL(kind_real) :: kval(nlev - 1) ! exponential decay rate between levels
69 REAL(kind_real) :: root_2pia
70 REAL(kind_real) :: ref_low ! refractivity at lower level
71 REAL(kind_real) :: nr_low ! impact parameter lower level
72 REAL(kind_real) :: tup
73 REAL(kind_real) :: tlow ! upper/lower bounds to error function
74 REAL(kind_real) :: diff_erf ! error function result
75 REAL(kind_real) :: dalpha ! delta bending angle
76 REAL(kind_real) :: erf_tup
77 REAL(kind_real) :: erf_tlow ! error function of tup and tlow
78 REAL(kind_real) :: t ! int. step in approx. to error function
79 REAL(kind_real), PARAMETER :: a1 = 0.3480242 ! consts used in error function approx.
80 REAL(kind_real), PARAMETER :: a2 = -0.0958798
81 REAL(kind_real), PARAMETER :: a3 = 0.7478556
82 REAL(kind_real), PARAMETER :: p = 0.47047
83 
84 jbot = 1
85 
86 DO
87 
88  IF (refrac(jbot) > 0.0 .AND. nr(jbot) > 0.0) EXIT
89 
90  jbot = jbot + 1
91 
92 END DO
93 
94 !-------------------------------------------------------------------------------
95 ! Calculate lowest usable level (because of superrefraction)
96 !-------------------------------------------------------------------------------
97 
98 kbot = nlev
99 
100 DO i = nlev,jbot + 1,-1
101 
102  ! to avoid large gradients
103  IF ((nr(kbot) - nr(kbot - 1)) < 10.0) EXIT
104 
105  kbot = kbot - 1
106 
107 END DO
108 
109 jbot = max(jbot,kbot)
110 
111 !-------------------------------------------------------------------------------
112 ! Calculate the exponential decay rate between levels
113 !-------------------------------------------------------------------------------
114 
115 DO i = jbot,nlev - 1
116 
117  kval(i) = log(refrac(i) / refrac(i + 1)) / &
118  max(1.0,(nr(i + 1) - nr(i)))
119  kval(i) = max(1.0e-6,kval(i))
120 
121 END DO
122 
123 !-------------------------------------------------------------------------------
124 ! Calculate the bending angle values
125 !-------------------------------------------------------------------------------
126 
127 alpha(:) = missing_value(alpha(1))
128 
129 DO n = 1,nobs
130 
131  IF (a(n) < nr(jbot) .OR. a(n) > nr(nlev)) cycle
132 
133  root_2pia = sqrt(2.0 * pi * a(n))
134 
135  ibot = jbot
136 
137  ! Find bottom state vector level
138  !----------------------------------
139 
140  DO
141  ! check more than 1 metre apart to stop large gradients in K code
142  IF (((nr(ibot + 1) - a(n)) > 1.0) .OR. ibot == nlev - 1) EXIT
143 
144  ibot = ibot + 1
145 
146  END DO
147 
148  ! Initialise bending angle value
149  !-------------------------------
150 
151  alpha(n) = 0.0
152 
153  ! Values of refractivity and impact parameter at lower level
154  !-----------------------------------------------------------
155 !The following line is a compiler directive to avoid bugs with intel compilers
156 !DIR$ NOVECTOR
157  DO i = ibot, nlev - 1
158 
159  IF (i == ibot) THEN
160 
161  ref_low = refrac(ibot) * exp(-kval(ibot) * (a(n) - nr(ibot)))
162  nr_low = a(n)
163 
164  ELSE
165 
166  ref_low = refrac(i)
167  nr_low = nr(i)
168 
169  END IF
170 
171  ! Limits used in the error function
172  !----------------------------------
173 
174  IF (i == nlev - 1) THEN
175 
176  ! Simple extrapolation 100km above the uppermost level
177  !-----------------------------------------------------
178  tup = sqrt(kval(i) * (nr(i + 1) + 1.0e5 - a(n)))
179 
180  ELSE
181 
182  tup = sqrt(kval(i) * (nr(i + 1) - a(n)))
183 
184  END IF
185 
186  tlow = 0.0
187 
188  IF (i > ibot) tlow = sqrt(kval(i) * (nr(i) - a(n)))
189 
190  ! Abramowitz and Stegun approx. to error function
191  !------------------------------------------------
192  t = 1.0 / (1.0 + p * tup)
193  erf_tup = 1.0 - exp(-(tup ** 2)) * (a1 + (a2 + a3 * t) * t) * t
194 
195  t = 1.0 / (1.0 + p * tlow)
196  erf_tlow = 1.0 - exp(-(tlow ** 2)) * (a1 + (a2 + a3 * t) * t) * t
197 
198  diff_erf = erf_tup - erf_tlow
199 
200  dalpha = 1.0e-6 * root_2pia * sqrt(kval(i)) * &
201  ref_low * exp(kval(i) * (nr_low - a(n))) * diff_erf
202 
203  ! Bending angle value
204  !--------------------
205  alpha(n) = alpha(n) + dalpha
206 
207  END DO
208 
209 END DO
210 
211 END SUBROUTINE ops_gpsrocalc_alpha
212 
213 
214 !-------------------------------------------------------------------------------
215 ! (C) Crown copyright Met Office. All rights reserved.
216 ! Refer to COPYRIGHT.txt of this distribution for details.
217 !-------------------------------------------------------------------------------
218 SUBROUTINE ops_gpsrocalc_nr (nb, & ! number of levels in zb
219  zb, & ! geopotential heights of model levels
220  refrac, & ! refractivity of model on model levels
221  Rad, & ! radius of curvature of earth at observation
222  lat, & ! latitude at observation
223  und, & ! geoid undulation above WGS-84
224  nr) ! Calculated model impact parameters
225 
226 IMPLICIT NONE
227 
228 ! Subroutine arguments:
229 INTEGER, INTENT(IN) :: nb ! number of levels in zb
230 REAL(kind_real), INTENT(IN) :: zb(nb) ! geopotential heights on zb levels /m
231 REAL(kind_real), INTENT(IN) :: rad ! local radius of curvature of earth /m
232 REAL(kind_real), INTENT(IN) :: lat ! latitude at observation/ degrees
233 REAL(kind_real), INTENT(IN) :: und ! geoid undulation
234 REAL(kind_real), INTENT(IN) :: refrac(nb) ! refractivity on model levels / N
235 REAL(kind_real), INTENT(OUT) :: nr(nb) ! Calculated model impact parameters
236 
237 ! Local declarations:
238 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_nr"
239 REAL(kind_real) :: r(nb) ! radius of model levels /m
240 REAL(kind_real) :: z(nb) ! geopotential heights on zb levels /m, local copy of zb
241 INTEGER :: i
242 
243 nr(:) = missing_value(nr(1))
244 
245 !----------------------------------------------
246 ! 1. Convert zb values to geometric altitudes
247 !---------------------------------------------
248 z = zb + und ! approx. convert to geopotential above WGS-84 ellipsoid
249 CALL ops_gpsro_geop_geom (lat, &
250  z)
251 
252 !--------------------------------------------------
253 ! 2. Calculate x =nr, i.e. model impact parameters
254 !--------------------------------------------------
255 
256 r = rad + z
257 
258 DO i = 1,nb
259  IF (zb(i) > 0.0 .AND. refrac(i) > 0.0) THEN
260  nr(i) = (1.0e0 + 1.0e-6 * refrac(i)) * r(i)
261  END IF
262 END DO
263 
264 END SUBROUTINE ops_gpsrocalc_nr
265 
266 
267 !-------------------------------------------------------------------------------
268 ! Convert geopotential height above ellispoid to geometric height on WGS-84
269 ! (ellipsoid). The method for this calculation is available from
270 ! http://mtp.jpl.nasa.gov/notes/altitude/altitude.html.
271 !-------------------------------------------------------------------------------
272 SUBROUTINE ops_gpsro_geop_geom (lat, &
273  z)
274 
275 IMPLICIT NONE
276 
277 ! Subroutine arguments:
278 REAL(kind_real), INTENT(IN) :: lat ! latitude of observation
279 REAL(kind_real), INTENT(INOUT) :: z(:) ! geopotential height in, geometric height out
280 
281 ! Local declarations:
282 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSRO_geop_geom"
283 REAL(kind_real) :: r_eff ! effective radius of Earth, formulas from MJ Mahoney (2005)
284 REAL(kind_real) :: g_somig ! Somigliana's equation for normal gravity on the surface of an ellipsoid of revolution
285 REAL(kind_real) :: latrad ! latitude in radians
286 
287 !------------------------
288 ! 1. correct units
289 !-----------------------
290 
291 latrad = lat * (pi / 180.0) !convert latitude from degrees to radians
292 
293 !-------------------------
294 ! 2. Calculate r and g
295 !-------------------------
296 
297 r_eff = a_earth / (1.0 + flatt + m_ratio - 2.0 * flatt * (sin(latrad)) ** 2)
298 
299 g_somig = g_equat * (1.0 + k_somig * (sin(latrad)) ** 2) / (sqrt(1.0 - (ecc ** 2) * (sin(latrad)) ** 2))
300 
301 !-------------------------------------------------------------------
302 ! 3. convert z (in geopotential height) to geometric wrt ellipsoid
303 !-------------------------------------------------------------------
304 
305 z(:) = (r_eff * z(:)) / ((g_somig / grav) * r_eff - z(:))
306 
307 END SUBROUTINE ops_gpsro_geop_geom
308 
309 !-------------------------------------------------------------------------
310 end module ufo_gnssro_ukmo1d_utils_mod
ufo_constants_mod
Definition: ufo_constants_mod.F90:2
ufo_constants_mod::rd
real(kind_real), parameter, public rd
Definition: ufo_constants_mod.F90:11
ufo_gnssro_ukmo1d_utils_mod
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:7
ufo_gnssro_ukmo1d_utils_mod::ops_gpsrocalc_nr
subroutine, public ops_gpsrocalc_nr(nb, zb, refrac, Rad, lat, und, nr)
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:225
ufo_gnssro_ukmo1d_utils_mod::ops_gpsrocalc_alpha
subroutine, public ops_gpsrocalc_alpha(nobs, nlev, a, refrac, nr, alpha)
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:50
ufo_gnssro_ukmo1d_utils_mod::ops_gpsro_geop_geom
subroutine, public ops_gpsro_geop_geom(lat, z)
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:274