UFO
ufo_gnssro_bendmetoffice_tlad_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_bendmetoffice_tlad_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_alphak
35 public :: ops_gpsrocalc_nrk
36 private
37 
38 contains
39 
40 SUBROUTINE ops_gpsrocalc_nrk (zb, & ! geopotential heights of model levels
41  nb, & ! number of levels in zb
42  Rad, & ! radius of curvature of earth at observation
43  lat, & ! latitude at observation
44  und, & ! geoid undulation above WGS-84
45  refrac, & ! refractivity of model on model levels
46  dnr_dref) ! Calculated gradient of nr
47 
48 
49 USE ufo_gnssro_ukmo1d_utils_mod, only: ops_gpsro_geop_geom
50 
51 IMPLICIT NONE
52 
53 ! Subroutine arguments:
54 INTEGER, INTENT(IN) :: nb ! number of levels in zb
55 REAL(kind_real), INTENT(IN) :: zb(nb) ! geopotential heights on zb levels /m
56 REAL(kind_real), INTENT(IN) :: rad ! local radius of curvature of earth /m
57 REAL(kind_real), INTENT(IN) :: lat ! latitude at observation/ degrees
58 REAL(kind_real), INTENT(IN) :: und ! geoid undulation
59 REAL(kind_real), INTENT(IN) :: refrac(nb) ! refractivity on model levels / N
60 REAL(kind_real), INTENT(OUT) :: dnr_dref(nb,nb) ! Calculated gradient of nr wrt ref
61 
62 ! Local declarations:
63 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_nrK"
64 REAL(kind_real) :: r(nb) ! radius of model levels /m
65 REAL(kind_real) :: z(nb) ! geopotential heights on zb levels /m, local copy of zb
66 INTEGER :: i
67 
68 !Initialise the matrix
69 dnr_dref = 0.0
70 
71 !----------------------------------------------
72 ! 1. Convert zb values to geometric altitudes
73 !---------------------------------------------
74 z= zb + und ! approx. convert to geopotential above WGS-84 ellipsoid
75 CALL ops_gpsro_geop_geom (lat, &
76  z)
77 
78 !--------------------------------------------------
79 ! 2. Calculate dnr/dref
80 !--------------------------------------------------
81 
82 r = rad + z
83 
84 DO i = 1,nb
85  IF (zb(i) > 0.0 .AND. refrac(i) > 0.0) THEN
86  dnr_dref(i,i) = 1.0e-6 * r(i)
87  END IF
88 END DO
89 
90 END SUBROUTINE ops_gpsrocalc_nrk
91 
92 !-------------------------------------------------------------------------------
93 ! GPSRO 1D bending angle operator K code.
94 !-------------------------------------------------------------------------------
95 SUBROUTINE ops_gpsrocalc_alphak (nobs, &
96  nlev, &
97  a, &
98  refrac, &
99  nr, &
100  Kmat_ref, &
101  Kmat_nr)
102 
103 IMPLICIT NONE
104 
105 ! Subroutine arguments:
106 INTEGER, INTENT(IN) :: nobs ! size of ob. vector
107 INTEGER, INTENT(IN) :: nlev ! no. of refractivity levels
108 REAL(kind_real), INTENT(IN) :: a(nobs) ! observation impact parameters
109 REAL(kind_real), INTENT(IN) :: refrac(nlev) ! refractivity values on model levels
110 REAL(kind_real), INTENT(IN) :: nr(nlev) ! refractive index * radius product
111 REAL(kind_real), INTENT(OUT) :: kmat_ref(nobs,nlev) ! BA gradient wrt refractivity
112 REAL(kind_real), INTENT(OUT) :: kmat_nr(nobs,nlev) ! BA gradient wrt index * radius product
113 
114 ! Local declarations:
115 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_alphaK"
116 INTEGER :: i
117 INTEGER :: n
118 INTEGER :: ibot
119 INTEGER :: jbot
120 INTEGER :: kbot
121 REAL(kind_real) :: kval(nlev - 1) ! exponential decay rate between levels
122 REAL(kind_real) :: root_2pia
123 REAL(kind_real) :: ref_low ! refractivity at lower level
124 REAL(kind_real) :: nr_low ! impact parameter lower level
125 REAL(kind_real) :: tup
126 REAL(kind_real) :: tlow ! upper/lower bounds to error function
127 REAL(kind_real) :: dalpha ! delta bending angle
128 REAL(kind_real) :: diff_erf ! error function result
129 REAL(kind_real) :: erf_tup
130 REAL(kind_real) :: erf_tlow ! error function of tup and tlow
131 REAL(kind_real) :: t ! int. step in approx. to error function
132 REAL(kind_real), PARAMETER :: a1 = 0.3480242 ! consts used in error function approx.
133 REAL(kind_real), PARAMETER :: a2 = -0.0958798
134 REAL(kind_real), PARAMETER :: a3 = 0.7478556
135 REAL(kind_real), PARAMETER :: p = 0.47047
136 REAL(kind_real) :: dkval_dref(nlev - 1,2)
137 REAL(kind_real) :: dkval_dnr(nlev - 1,2)
138 REAL(kind_real) :: dalpha_dref(2)
139 REAL(kind_real) :: dalpha_dnr(2)
140 REAL(kind_real) :: dalpha_dk
141 REAL(kind_real) :: dalpha_drlow
142 REAL(kind_real) :: dalpha_derf
143 REAL(kind_real) :: dalpha_dnrlow
144 REAL(kind_real) :: dnrlow_dref(2)
145 REAL(kind_real) :: dnrlow_dnr(2)
146 REAL(kind_real) :: drlow_dref(2)
147 REAL(kind_real) :: drlow_dnr(2)
148 REAL(kind_real) :: drlow_dk
149 REAL(kind_real) :: derf_dref(2)
150 REAL(kind_real) :: derf_dnr(2)
151 REAL(kind_real) :: derf_dtup
152 REAL(kind_real) :: derf_dtlow
153 REAL(kind_real) :: dtup_dnr(2)
154 REAL(kind_real) :: dtlow_dnr(2)
155 REAL(kind_real) :: dtup_dref(2)
156 REAL(kind_real) :: dtlow_dref(2)
157 REAL(kind_real) :: dtup_dk
158 REAL(kind_real) :: dtlow_dk
159 
160 !-------------------------------------------------------------------------------
161 ! Initialise the K matrices
162 !-------------------------------------------------------------------------------
163 
164 kmat_ref(:,:) = 0.0
165 kmat_nr(:,:) = 0.0
166 dkval_dref(:,:) = 0.0
167 dkval_dnr(:,:) = 0.0
168 
169 jbot = 1
170 
171 DO
172 
173  IF (refrac(jbot) > 0.0 .AND. nr(jbot) > 0.0) EXIT
174 
175  jbot = jbot + 1
176 
177 END DO
178 
179 !-------------------------------------------------------------------------------
180 ! Calculate lowest usable level (because of superrefraction)
181 !-------------------------------------------------------------------------------
182 
183 kbot = nlev
184 
185 DO i = nlev,jbot + 1,-1
186 
187  ! to avoid large gradients
188  IF ((nr(kbot) - nr(kbot-1)) < 10.0) EXIT
189 
190  kbot = kbot - 1
191 
192 END DO
193 
194 jbot = max(jbot,kbot)
195 
196 !-------------------------------------------------------------------------------
197 ! Calculate the exponential decay rate between levels
198 !-------------------------------------------------------------------------------
199 
200 DO i = jbot,nlev - 1
201 
202  kval(i) = log(refrac(i) / refrac(i + 1)) / &
203  max(1.0,(nr(i + 1) - nr(i)))
204 
205  kval(i) = max(1.0e-6,kval(i))
206 
207  IF (kval(i) > 1.0e-6) THEN
208 
209  dkval_dref(i,1) = 1.0 / (refrac(i) * max(1.0,(nr(i + 1) - nr(i))))
210  dkval_dref(i,2) = -1.0 / (refrac(i + 1) * max(1.0,(nr(i + 1) - nr(i))))
211 
212  dkval_dnr(i,1) = kval(i) / max(1.0,(nr(i + 1) - nr(i)))
213  dkval_dnr(i,2) = -kval(i) / max(1.0,(nr(i + 1) - nr(i)))
214 
215  END IF
216 
217 END DO
218 
219 !-------------------------------------------------------------------------------
220 ! Calculate the bending angle gradients
221 !-------------------------------------------------------------------------------
222 
223 DO n = 1,nobs
224 
225  IF (a(n) < nr(jbot) .OR. a(n) > nr(nlev)) cycle
226 
227  root_2pia = sqrt(2.0 * pi * a(n))
228 
229  ibot = jbot
230 
231  ! Find bottom state vector level
232  !----------------------------------
233  DO
234 
235  ! check more than 1 metre apart to stop large gradients in K code
236  ! ---------------------------------------------------------------
237  IF (((nr(ibot + 1) - a(n)) > 1.0) .OR. ibot == nlev - 1) EXIT
238 
239  ibot = ibot + 1
240 
241  END DO
242 
243  tlow = 0.0
244 
245  DO i = ibot, nlev - 1
246 
247  ! initialise matrices
248  !---------------------
249 
250  dalpha_dref(:) = 0.0
251  dalpha_dnr(:) = 0.0
252  drlow_dref(:) = 0.0
253  drlow_dnr(:) = 0.0
254  dnrlow_dref(:) = 0.0
255  dnrlow_dnr(:) = 0.0
256  dtup_dnr(:) = 0.0
257  dtup_dref(:) = 0.0
258  dtlow_dnr(:) = 0.0
259  dtlow_dref(:) = 0.0
260  derf_dtup = 0.0
261  derf_dtlow = 0.0
262  dalpha_drlow = 0.0
263  dalpha_dnrlow = 0.0
264  dalpha_dk = 0.0
265  dalpha_derf = 0.0
266  drlow_dk = 0.0
267  dtup_dk = 0.0
268  dtlow_dk = 0.0
269 
270  ! Values of refractivity and impact parameter at lower level
271  !-----------------------------------------------------------
272  IF (i == ibot) THEN
273 
274  ref_low = refrac(i) * exp(-kval(i) * (a(n) - nr(i)))
275 
276  drlow_dref(1)= ref_low / refrac(i)
277  drlow_dk = -ref_low * (a(n) - nr(i))
278  drlow_dnr(1) = ref_low * kval(i)
279 
280  nr_low = a(n)
281 
282  dnrlow_dnr(1) = 0.0
283 
284  ELSE
285 
286  ref_low = refrac(i)
287 
288  drlow_dref(1) = 1.0
289 
290  nr_low = nr(i)
291 
292  dnrlow_dnr(1) = 1.0
293 
294  END IF
295 
296  drlow_dref(:) = drlow_dref(:) + drlow_dk * dkval_dref(i,:)
297  drlow_dnr(:) = drlow_dnr(:) + drlow_dk * dkval_dnr(i,:)
298 
299 
300  ! Limits used in the error function
301  !----------------------------------
302  IF (i == nlev - 1) THEN
303 
304  ! simple extrapolation 100km above the uppermost level.
305  !-----------------------------------------------------
306  tup = sqrt(kval(i) * (nr(i + 1) + 1.0e5 - a(n)))
307 
308  dtup_dk = 0.5 * (nr(i + 1) + 1.0e5 - a(n)) / tup
309  dtup_dnr(2) = 0.5 * kval(i) / tup
310 
311  ELSE
312 
313  tup = sqrt(kval(i) * (nr(i + 1) - a(n)))
314 
315  dtup_dk = 0.5 * (nr(i + 1) - a(n)) / tup
316  dtup_dnr(2) = 0.5 * kval(i) / tup
317 
318  END IF
319 
320  dtup_dref(:) = dtup_dref(:) + dtup_dk * dkval_dref(i,:)
321  dtup_dnr(:) = dtup_dnr(:) + dtup_dk * dkval_dnr(i,:)
322 
323  tlow = 0.0
324 
325  IF (i > ibot) THEN
326  tlow = sqrt(kval(i) * (nr(i) - a(n)))
327  dtlow_dk = 0.5 * (nr(i) - a(n)) / tlow
328  dtlow_dnr(1) = 0.5 * kval(i) / tlow
329  END IF
330 
331  dtlow_dref(:) = dtlow_dref(:) + dtlow_dk * dkval_dref(i,:)
332  dtlow_dnr(:) = dtlow_dnr(:) + dtlow_dk * dkval_dnr(i,:)
333 
334  ! Abramowitz and Stegun approx. to error function
335  !------------------------------------------------
336  t = 1.0 / (1.0 + p * tup)
337  erf_tup = 1.0 - exp(-(tup ** 2)) * (a1 + (a2 + a3 * t) * t) * t
338  derf_dtup = exp(-(tup ** 2)) * ((2.0 * tup * (a1 + (a2 + a3 * t) * t) * t) + &
339  ((a1 + (2.0 * a2 + 3.0 * a3 * t) * t) * p * t ** 2))
340 
341  t = 1.0 / (1.0 + p * tlow)
342  erf_tlow = 1.0 - exp(-(tlow ** 2)) * (a1 + (a2 + a3 * t) * t) * t
343 
344  ! Multiplied by -1.0 to account for the usage in derf_dref and derf_dnr
345  !---------------------------------------------------------------------
346  derf_dtlow = -1.0 * exp(-(tlow ** 2)) * ((2.0 * tlow * (a1 + (a2 + a3 * t) * t) * t) + &
347  ((a1 + (2.0 * a2 + 3.0 * a3 * t) * t) * p * t ** 2))
348 
349  diff_erf = erf_tup - erf_tlow
350 
351  derf_dref(:) = derf_dtup * dtup_dref(:) + &
352  derf_dtlow * dtlow_dref(:)
353 
354  derf_dnr(:) = derf_dtup * dtup_dnr(:) + &
355  derf_dtlow * dtlow_dnr(:)
356 
357  dalpha = 1.0e-6 * root_2pia * sqrt(kval(i)) * &
358  ref_low * exp(kval(i) * (nr_low - a(n))) * diff_erf
359 
360  dalpha_drlow = dalpha / max(1.0e-10,ref_low)
361 
362  dalpha_derf = dalpha / max(1.0e-10,diff_erf)
363 
364  dalpha_dnrlow = dalpha * kval(i)
365 
366  dalpha_dk = dalpha * (nr_low - a(n) + 0.5 / kval(i))
367 
368  ! Now apply chain rule
369  !----------------------
370 
371  dalpha_dref(:) = dalpha_dref(:) + &
372  dalpha_drlow * drlow_dref(:) + &
373  dalpha_derf * derf_dref(:) + &
374  dalpha_dnrlow * dnrlow_dref(:) + &
375  dalpha_dk * dkval_dref(i,:)
376 
377  dalpha_dnr(:) = dalpha_dnr(:) + &
378  dalpha_drlow * drlow_dnr(:) + &
379  dalpha_derf * derf_dnr(:) + &
380  dalpha_dnrlow * dnrlow_dnr(:) + &
381  dalpha_dk * dkval_dnr(i,:)
382 
383  ! Now update matrices
384  !---------------------
385 
386  kmat_ref(n,i) = kmat_ref(n,i) + dalpha_dref(1)
387  kmat_nr(n,i) = kmat_nr(n,i) + dalpha_dnr(1)
388 
389  kmat_ref(n,i + 1) = kmat_ref(n,i + 1) + dalpha_dref(2)
390  kmat_nr(n,i + 1) = kmat_nr(n,i + 1) + dalpha_dnr(2)
391 
392  END DO
393 
394 END DO
395 
396 END SUBROUTINE ops_gpsrocalc_alphak
397 
398 end module ufo_gnssro_bendmetoffice_tlad_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_bendmetoffice_tlad_utils_mod
Definition: ufo_gnssro_bendmetoffice_tlad_utils_mod.F90:7
ufo_gnssro_bendmetoffice_tlad_utils_mod::ops_gpsrocalc_nrk
subroutine, public ops_gpsrocalc_nrk(zb, nb, Rad, lat, und, refrac, dnr_dref)
Definition: ufo_gnssro_bendmetoffice_tlad_utils_mod.F90:47
ufo_gnssro_bendmetoffice_tlad_utils_mod::ops_gpsrocalc_alphak
subroutine, public ops_gpsrocalc_alphak(nobs, nlev, a, refrac, nr, Kmat_ref, Kmat_nr)
Definition: ufo_gnssro_bendmetoffice_tlad_utils_mod.F90:102
ufo_gnssro_ukmo1d_utils_mod
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:7
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