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 public :: ops_gpsro_refrack
37 private
38 
39 contains
40 
41 !-------------------------------------------------------------------------------
42 ! (C) Crown copyright Met Office. All rights reserved.
43 ! Refer to COPYRIGHT.txt of this distribution for details.
44 !-------------------------------------------------------------------------------
45 ! Calculate GPSRO refractivity K matrix.
46 !-------------------------------------------------------------------------------
47 
48 SUBROUTINE ops_gpsro_refrack (nstate, &
49  nlevP, &
50  nb, &
51  nlevq, &
52  za, &
53  zb, &
54  x, &
55  GPSRO_pseudo_ops, &
56  GPSRO_vert_interp_ops, &
57  dref_dP, &
58  dref_dq)
59 
60 IMPLICIT NONE
61 
62 ! Subroutine arguments:
63 INTEGER, INTENT(IN) :: nstate ! Number of variables in full model state
64 INTEGER, INTENT(IN) :: nlevp ! Number of p levels in state vector
65 INTEGER, INTENT(IN) :: nb ! Number of theta levels
66 INTEGER, INTENT(IN) :: nlevq ! Number of specific humidity levels (=nb)
67 REAL(kind_real), INTENT(IN) :: za(:) ! Heights of the pressure levels
68 REAL(kind_real), INTENT(IN) :: zb(:) ! Heights of the theta (specific humidity) levels
69 REAL(kind_real), INTENT(IN) :: x(:) ! Input model state
70 LOGICAL, INTENT(IN) :: gpsro_pseudo_ops ! Whether to use pseudo-levels in the calculation
71 LOGICAL, INTENT(IN) :: gpsro_vert_interp_ops ! Whether to use log(p) for the vertical interpolation
72 REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: dref_dp(:,:) ! K-matrix for p
73 REAL(kind_real), ALLOCATABLE, INTENT(OUT) :: dref_dq(:,:) ! K-matrix for q
74 
75 ! Local declarations:
76 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSRO_refracK"
77 INTEGER :: i
78 INTEGER :: counter
79 REAL(kind_real) :: p(nlevp)
80 REAL(kind_real) :: exner(nlevp)
81 REAL(kind_real) :: q(nlevq)
82 REAL(kind_real) :: pb(nlevq)
83 REAL(kind_real) :: tv(nlevq)
84 REAL(kind_real) :: t(nlevq)
85 REAL(kind_real) :: extheta
86 REAL(kind_real) :: pwt1
87 REAL(kind_real) :: pwt2
88 REAL(kind_real) :: ndry
89 REAL(kind_real) :: nwet
90 REAL(kind_real) :: refrac(nlevq)
91 REAL(kind_real) :: dex_dp(nlevp,nlevp)
92 REAL(kind_real) :: dpb_dp(nlevq,nlevp)
93 REAL(kind_real) :: dextheta_dpb(nlevq,nlevq)
94 REAL(kind_real) :: dtv_dextheta(nlevq,nlevq)
95 REAL(kind_real) :: dtv_dex(nlevq,nlevp)
96 REAL(kind_real) :: dt_dtv(nlevq,nlevq)
97 REAL(kind_real) :: dt_dq(nlevq,nlevq)
98 REAL(kind_real) :: dref_dpb(nlevq,nlevq)
99 REAL(kind_real) :: dref_dt(nlevq,nlevq)
100 REAL(kind_real) :: m1(nb,nlevq)
101 REAL(kind_real) :: m2(nb,nlevp)
102 REAL(kind_real) :: m3(nb,nlevq)
103 REAL(kind_real) :: m4(nb,nlevq)
104 REAL(kind_real), ALLOCATABLE :: p_pseudo(:)
105 REAL(kind_real), ALLOCATABLE :: q_pseudo(:)
106 REAL(kind_real), ALLOCATABLE :: t_pseudo(:)
107 REAL(kind_real), ALLOCATABLE :: z_pseudo(:)
108 REAL(kind_real), ALLOCATABLE :: n_pseudo(:)
109 REAL(kind_real) :: gamma
110 REAL(kind_real) :: beta
111 REAL(kind_real) :: c !continuity constant for hydrostatic pressure
112 REAL(kind_real) :: g_rb ! Frequently used term
113 REAL(kind_real) :: c_zz ! Frequently used term
114 REAL(kind_real) :: dpp_dt1 ! dP_pseudo / dT_below
115 REAL(kind_real) :: dpp_dtp ! dP_pseudo / dT_pseudo
116 REAL(kind_real) :: dtp_dt1 ! dT_pseudo / dT_below
117 REAL(kind_real) :: dpp_dbeta ! dP_pseudo / dbeta
118 REAL(kind_real) :: dbeta_dt1 ! dbeta / dT_below
119 REAL(kind_real) :: dtp_dbeta ! dT_pseudo / dbeta
120 REAL(kind_real) :: dbeta_dt2 ! dbeta / dT_above
121 REAL(kind_real) :: dpp_dc ! dP_pseudo / dc
122 REAL(kind_real) :: dc_dt1 ! dc / dT_below
123 REAL(kind_real) :: dc_dbeta ! dc / dbeta
124 REAL(kind_real) :: dc_dt2 ! dc / dT_above
125 REAL(kind_real) :: dpp_dp1 ! dP_pseudo / dP_below
126 REAL(kind_real) :: dc_dp1 ! dc / dP_below
127 REAL(kind_real) :: dc_dp2 ! dc / dP_above
128 REAL(kind_real), ALLOCATABLE :: dref_dppseudo(:,:)
129 REAL(kind_real), ALLOCATABLE :: dref_dtpseudo(:,:)
130 REAL(kind_real), ALLOCATABLE :: dref_dqpseudo(:,:)
131 REAL(kind_real), ALLOCATABLE :: dppseudo_dpb(:,:)
132 REAL(kind_real), ALLOCATABLE :: dtpseudo_dtb(:,:)
133 REAL(kind_real), ALLOCATABLE :: dqpseudo_dqb(:,:)
134 REAL(kind_real), ALLOCATABLE :: dppseudo_dt(:,:)
135 REAL(kind_real), ALLOCATABLE :: dtb_dp(:,:)
136 REAL(kind_real), ALLOCATABLE :: dppseudo_dp(:,:)
137 REAL(kind_real), ALLOCATABLE :: dtpseudo_dp(:,:)
138 
139 IF (gpsro_pseudo_ops) THEN
140  ALLOCATE (dref_dp(nb,nlevp))
141  ALLOCATE (dref_dq(nb,nlevq))
142 
143  ALLOCATE (p_pseudo(nb))
144  ALLOCATE (q_pseudo(nb))
145  ALLOCATE (t_pseudo(nb))
146  ALLOCATE (z_pseudo(nb))
147  ALLOCATE (n_pseudo(nb))
148 
149  ALLOCATE (dref_dppseudo(nb,nb))
150  ALLOCATE (dref_dtpseudo(nb,nb))
151  ALLOCATE (dref_dqpseudo(nb,nb))
152  ALLOCATE (dppseudo_dpb(nb,nlevq))
153  ALLOCATE (dppseudo_dt(nb,nlevq))
154  ALLOCATE (dtpseudo_dtb(nb,nlevq))
155  ALLOCATE (dqpseudo_dqb(nb,nlevq))
156 
157  ALLOCATE (dtb_dp(nlevq,nlevp))
158  ALLOCATE (dppseudo_dp(nb,nlevp))
159  ALLOCATE (dtpseudo_dp(nb,nlevp))
160 
161  dref_dppseudo(:,:) = 0.0
162  dref_dtpseudo(:,:) = 0.0
163  dref_dqpseudo(:,:) = 0.0
164  dppseudo_dpb(:,:) = 0.0
165  dppseudo_dt(:,:) = 0.0
166  dtpseudo_dtb(:,:) = 0.0
167  dqpseudo_dqb(:,:) = 0.0
168  dppseudo_dp(:,:) = 0.0
169  dtpseudo_dp(:,:) = 0.0
170 ELSE
171  ALLOCATE (dref_dp(nlevq,nlevp))
172  ALLOCATE (dref_dq(nlevq,nlevq))
173 END IF
174 
175 !-----------------------
176 ! 1. Initialise matrices
177 !-----------------------
178 
179 dpb_dp(:,:) = 0.0
180 dextheta_dpb(:,:) = 0.0
181 dex_dp(:,:) = 0.0
182 dtv_dextheta(:,:) = 0.0
183 dtv_dex(:,:) = 0.0
184 dt_dtv(:,:) = 0.0
185 dt_dq(:,:) = 0.0
186 dref_dpb(:,:) = 0.0
187 dref_dt(:,:) = 0.0
188 dref_dq(:,:) = 0.0
189 dref_dp(:,:) = 0.0
190 
191 ! Set up the P and q vectors from x
192 
193 p(:) = 1.0e2 * x(1:nlevp)
194 q(:) = 1.0e-3 * x(nlevp + 1:nstate)
195 
196 ! Calculate exner on rho levels.
197 
198 exner(:) = (p(:) / pref) ** rd_over_cp
199 
200 DO i = 1,nlevp
201 
202  dex_dp(i,i) = rd_over_cp / pref * (p(i) / pref) ** (rd_over_cp - 1.0)
203 
204 END DO
205 
206 !----------------------------------------------
207 ! 2. Calculate the refractivity on the b levels
208 !----------------------------------------------
209 
210 DO i = 1, nlevq
211 
212  ! Calc. pressure on b levels
213 
214  pwt1 = (za(i + 1) - zb(i)) / (za(i + 1) - za(i))
215 
216  pwt2 = 1.0 - pwt1
217 
218  ! calculate the pressure on the theta level.
219  IF (gpsro_vert_interp_ops) THEN
220  pb(i) = exp(pwt1 * log(p(i)) + pwt2 * log(p(i + 1)))
221 
222  dpb_dp(i,i) = pb(i) * pwt1 / p(i)
223  dpb_dp(i,i + 1) = pb(i) * pwt2 / p(i + 1)
224  ELSE
225  ! Assume Exner varies linearly with height
226  pb(i) = pref * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp)
227 
228  dpb_dp(i,i) = pwt1 * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * &
229  (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp - 1.0) * (p(i) / pref) ** (rd_over_cp - 1.0)
230  dpb_dp(i,i + 1) = pwt2 * (pwt1 * (p(i) / pref) ** rd_over_cp + pwt2 * &
231  (p(i + 1) / pref) ** rd_over_cp) ** (1.0 / rd_over_cp - 1.0) * (p(i + 1) / pref) ** (rd_over_cp-1.0)
232  END IF
233 
234  ! calculate Exner on the theta level.
235 
236  extheta = (pb(i) / pref) ** rd_over_cp
237 
238  dextheta_dpb(i,i) = rd_over_cp * (pb(i) ** (rd_over_cp - 1.0)) / (pref ** rd_over_cp)
239 
240  ! Calculate mean layer Tv using ND definition
241 
242  tv(i) = grav * (za(i + 1) - za(i)) * extheta / (cp * (exner(i) - exner(i + 1)))
243 
244  dtv_dextheta(i,i) = tv(i) / extheta
245 
246  dtv_dex(i,i) = -tv(i) / (exner(i) - exner(i + 1))
247 
248  dtv_dex(i,i + 1) = tv(i) / (exner(i) - exner(i + 1))
249 
250  IF (i > nlevq) THEN
251 
252  t(i) = tv(i)
253 
254  dt_dtv(i,i) = 1.0
255 
256  ! no wet component
257 
258  nwet = 0.0
259 
260  ELSE
261 
262  t(i) = tv(i) / (1.0 + c_virtual * q(i))
263 
264  dt_dtv(i,i) = 1.0 / (1.0 + c_virtual * q(i))
265 
266  dt_dq(i,i) = -c_virtual * t(i) / (1.0 + c_virtual * q(i))
267 
268  ! wet compontent
269 
270  nwet = n_beta * pb(i) * q(i) / (t(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q(i)))
271 
272  dref_dq(i,i) = n_beta * pb(i) * mw_ratio / (t(i) * (mw_ratio + (1.0 - mw_ratio) * q(i))) ** 2
273 
274  END IF
275 
276  ndry = n_alpha * pb(i) / t(i)
277 
278  refrac(i) = ndry + nwet
279 
280  dref_dpb(i,i) = refrac(i) / pb(i)
281 
282  dref_dt(i,i) = -(ndry + 2.0 * nwet) / t(i)
283 
284 END DO
285 
286 IF (gpsro_pseudo_ops) THEN
287  !----------------------------------!
288  !- Add intermediate pseudo-levels -!
289  !----------------------------------!
290  dtb_dp = matmul(dt_dtv, matmul(matmul(dtv_dextheta, dextheta_dpb), dpb_dp) + matmul(dtv_dex, dex_dp))
291  counter = 1
292  DO i = 1, nb
293  ! Odd 'i' (i.e. copies of actual model level values)
294  IF (mod(i, 2) > 0) THEN
295  z_pseudo(i) = zb(counter)
296  p_pseudo(i) = pb(counter)
297  q_pseudo(i) = q(counter)
298  t_pseudo(i) = t(counter)
299  n_pseudo(i) = n_alpha * p_pseudo(i) / t_pseudo(i) + n_beta * p_pseudo(i) * q_pseudo(i) / &
300  (t_pseudo(i) ** 2 * (mw_ratio + (1.0 - mw_ratio) * q_pseudo(i)))
301 
302  dref_dppseudo(i,i) = dref_dpb(counter,counter)
303  dref_dtpseudo(i,i) = dref_dt(counter,counter)
304  dref_dqpseudo(i,i) = dref_dq(counter,counter)
305 
306  dppseudo_dpb(i,counter) = 1.0
307  dtpseudo_dtb(i,counter) = 1.0
308  dqpseudo_dqb(i,counter) = 1.0
309 
310  counter = counter + 1
311 
312  ! Even 'i' (i.e. intermediate pseudo-levels)
313  ELSE
314  z_pseudo(i) = (zb(counter - 1) + zb(counter)) / 2.0
315  IF (min(q(counter - 1), q(counter)) > 0.0) THEN
316  gamma = log(q(counter - 1) / q(counter)) / (zb(counter) - zb(counter - 1))
317  q_pseudo(i) = q(counter - 1) * exp(-gamma * (z_pseudo(i) - z_pseudo(i - 1)))
318  ELSE
319  q_pseudo(i) = q(counter - 1) + (q(counter) - q(counter - 1)) / &
320  (zb(counter) - zb(counter - 1)) * (z_pseudo(i) - zb(counter - 1))
321  END IF
322 
323  beta = (t(counter) - t(counter - 1)) / (zb(counter) - zb(counter - 1))
324  t_pseudo(i) = t(counter - 1) + beta * (z_pseudo(i) - zb(counter - 1))
325  IF (abs(t(counter) - t(counter - 1)) > 1.0e-10) THEN
326  c = ((pb(counter) / pb(counter - 1)) * (t(counter) / t(counter - 1)) ** (grav / (rd * beta)) - 1.0) / &
327  (zb(counter) - zb(counter - 1))
328  p_pseudo(i) = (pb(counter - 1) * (t_pseudo(i) / t(counter - 1)) ** &
329  (-grav / (rd * beta))) * (1.0 + c * (z_pseudo(i) - zb(counter - 1)))
330  ELSE
331  p_pseudo(i) = pb(counter - 1) * exp(log(pb(counter) / pb(counter - 1)) * &
332  ((z_pseudo(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))))
333  END IF
334 
335  ndry = n_alpha * p_pseudo(i) / t_pseudo(i)
336  nwet = n_beta * p_pseudo(i) * q_pseudo(i) / (t_pseudo(i) ** 2 &
337  *(mw_ratio + (1.0 - mw_ratio) * q_pseudo(i)))
338 
339  n_pseudo(i) = ndry + nwet
340 
341  dref_dppseudo(i,i) = n_pseudo(i) / p_pseudo(i)
342  dref_dtpseudo(i,i) = -(ndry + 2.0 * nwet) / t_pseudo(i)
343  dref_dqpseudo(i,i) = n_beta * p_pseudo(i) * mw_ratio / (t_pseudo(i) * (mw_ratio + &
344  (1.0 - mw_ratio) * q_pseudo(i))) ** 2
345 
346  ! Transform P, q and T to pseudo-levels
347  IF (min(q(counter - 1), q(counter)) > 0.0) THEN
348  dqpseudo_dqb(i,counter - 1) = (q_pseudo(i) / q(counter - 1)) * (1.0 - &
349  (z_pseudo(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1)))
350 
351  dqpseudo_dqb(i,counter) = (q_pseudo(i) / q(counter)) * ((z_pseudo(i) - zb(counter - 1)) / &
352  (zb(counter) - zb(counter - 1)))
353  ELSE
354  dqpseudo_dqb(i,counter - 1) = (zb(counter) - z_pseudo(i)) / (zb(counter) - zb(counter - 1))
355 
356  dqpseudo_dqb(i,counter) = (z_pseudo(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))
357  END IF
358 
359  dtpseudo_dtb(i,counter - 1) = (zb(counter) - z_pseudo(i)) / (zb(counter) - zb(counter - 1))
360 
361  dtpseudo_dtb(i,counter) = (z_pseudo(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))
362 
363  ! Items that need changing for dPpseudo/dTb:
364 
365  g_rb = grav / (rd * beta)
366  c_zz = c + 1.0 / (zb(counter) - zb(counter - 1))
367 
368  dpp_dt1 = (g_rb / t(counter - 1)) * p_pseudo(i)
369  dpp_dtp = -(g_rb / t_pseudo(i)) * p_pseudo(i)
370  dtp_dt1 = dtpseudo_dtb(i,counter - 1)
371  dpp_dbeta = (g_rb / beta) * log(t_pseudo(i) / t(counter - 1)) * p_pseudo(i)
372  dbeta_dt1 = -1.0 / (zb(counter) - zb(counter - 1))
373  dtp_dbeta = z_pseudo(i) - zb(counter - 1)
374  dbeta_dt2 = 1.0 / (zb(counter) - zb(counter - 1))
375 
376  IF (abs(t(counter) - t(counter - 1)) > 1.0e-10) THEN
377  ! Incomplete computation of dPpseudo_dPb. Temperature derivatives done below
378  dpp_dc = (z_pseudo(i) - zb(counter - 1)) * pb(counter - 1) * (t_pseudo(i) / t(counter - 1)) ** (-g_rb)
379  dc_dt1 = -(g_rb / (t(counter - 1))) * c_zz
380  dc_dbeta = -(g_rb / beta) * log(t(counter) / t(counter - 1)) * c_zz
381  dc_dt2 = (g_rb / t(counter)) * c_zz
382  dpp_dp1 = p_pseudo(i) / pb(counter - 1)
383  dppseudo_dt(i,counter - 1) = dpp_dt1 + dpp_dtp * dtp_dt1 + dpp_dbeta * dbeta_dt1 + dpp_dc * &
384  (dc_dt1 + dc_dbeta * dbeta_dt1)
385  dppseudo_dt(i,counter) = (dpp_dbeta + dpp_dtp * dtp_dbeta + dpp_dc * dc_dbeta) * &
386  dbeta_dt2 + dpp_dc * dc_dt2
387 
388  dc_dp1 = -(1.0 / pb(counter - 1)) * c_zz
389  dc_dp2 = (1.0 / pb(counter)) * c_zz
390 
391  dppseudo_dpb(i,counter - 1) = dpp_dp1 + dpp_dc * dc_dp1
392  dppseudo_dpb(i,counter) = dpp_dc * dc_dp2
393  ELSE
394  dppseudo_dpb(i,counter - 1) = exp(log(pb(counter) / pb(counter - 1)) * ((z_pseudo(i) - zb(counter - 1)) / &
395  (zb(counter) - zb(counter - 1)))) * (1.0 - ((z_pseudo(i) - zb(counter - 1)) / (zb(counter) - zb(counter - 1))))
396  dppseudo_dpb(i,counter) = (pb(counter - 1) / pb(counter)) * ((z_pseudo(i) - zb(counter - 1)) / &
397  (zb(counter) - zb(counter - 1))) * exp(log(pb(counter) / pb(counter - 1)) * ((z_pseudo(i) - zb(counter - 1)) / &
398  (zb(counter) - zb(counter - 1))))
399  END IF
400 
401  END IF
402 
403  END DO
404 
405  ! temperature derivatives:
406  dppseudo_dp = matmul(dppseudo_dpb, dpb_dp) + matmul(dppseudo_dt, dtb_dp)
407  dtpseudo_dp = matmul(dtpseudo_dtb, dtb_dp)
408 
409  !-------------------------------------------------
410  ! 3. Evaluate the Kmatrix by matrix multiplication
411  !-------------------------------------------------
412 
413  ! calc K matrix for P on rho levels
414  dref_dp(:,:) = matmul(dref_dppseudo, dppseudo_dp) + matmul(dref_dtpseudo, dtpseudo_dp)
415 
416  ! calc Kmatrix for q on theta levels
417  dref_dq(:,:) = matmul(dref_dqpseudo, dqpseudo_dqb) + matmul(matmul(dref_dtpseudo, dtpseudo_dtb), dt_dq) + &
418  matmul(matmul(dref_dppseudo, dppseudo_dt), dt_dq)
419 
420  ! Put the K matrices in correct units
421 
422  dref_dp(:,:) = 1.0e2 * dref_dp(:,:) ! hPa
423 
424  dref_dq(:,:) = 1.0e-3 * dref_dq(:,:) ! g/kg
425 
426 ! Normal model levels
427 ELSE
428 
429  !-------------------------------------------------
430  ! 3. Evaluate the Kmatrix by matrix multiplication
431  !-------------------------------------------------
432 
433  ! calc K matrix for P on rho levels
434  ! dNmod/dP = (dNmod/dPb * dPb/dP) + ....
435  dref_dp(:,:) = matmul(dref_dpb, dpb_dp)
436 
437  ! .... (dNmod/dT * dT/dTv * dTv/dEx *dEx/dP) + .....
438  m1(:,:) = matmul(dref_dt, dt_dtv)
439  m2(:,:) = matmul(m1, dtv_dex)
440  dref_dp(:,:) = dref_dp(:,:) + matmul(m2, dex_dp)
441 
442  ! .... (dNmod/dT * dT/dTv * dTv/dExtheta *dExtheta/dPb*dPb/dP)
443  m3(:,:) = matmul(m1, dtv_dextheta)
444  m4(:,:) = matmul(m3, dextheta_dpb)
445  dref_dp(:,:) = dref_dp(:,:) + matmul(m4, dpb_dp)
446 
447  ! calc Kmatrix for q on theta levels
448  ! dNmod/dq = (dNmod/dq) + (dNmod/dT*dT/dq)
449  dref_dq(:,:) = dref_dq(:,:) + matmul(dref_dt, dt_dq)
450 
451  ! Put the K matrices in correct units
452 
453  dref_dp(:,:) = 1.0e2 * dref_dp(:,:) ! hPa
454 
455  dref_dq(:,:) = 1.0e-3 * dref_dq(:,:) ! g/kg
456 
457 END IF
458 
459 IF (ALLOCATED (p_pseudo)) DEALLOCATE (p_pseudo)
460 IF (ALLOCATED (q_pseudo)) DEALLOCATE (q_pseudo)
461 IF (ALLOCATED (t_pseudo)) DEALLOCATE (t_pseudo)
462 IF (ALLOCATED (z_pseudo)) DEALLOCATE (z_pseudo)
463 IF (ALLOCATED (n_pseudo)) DEALLOCATE (n_pseudo)
464 IF (ALLOCATED (dref_dppseudo)) DEALLOCATE (dref_dppseudo)
465 IF (ALLOCATED (dref_dtpseudo)) DEALLOCATE (dref_dtpseudo)
466 IF (ALLOCATED (dref_dqpseudo)) DEALLOCATE (dref_dqpseudo)
467 IF (ALLOCATED (dppseudo_dpb)) DEALLOCATE (dppseudo_dpb)
468 IF (ALLOCATED (dtpseudo_dtb)) DEALLOCATE (dtpseudo_dtb)
469 IF (ALLOCATED (dqpseudo_dqb)) DEALLOCATE (dqpseudo_dqb)
470 IF (ALLOCATED (dppseudo_dt)) DEALLOCATE (dppseudo_dt)
471 IF (ALLOCATED (dtb_dp)) DEALLOCATE (dtb_dp)
472 IF (ALLOCATED (dppseudo_dp)) DEALLOCATE (dppseudo_dp)
473 IF (ALLOCATED (dtpseudo_dp)) DEALLOCATE (dtpseudo_dp)
474 
475 END SUBROUTINE ops_gpsro_refrack
476 
477 
478 SUBROUTINE ops_gpsrocalc_nrk (zb, & ! geopotential heights of model levels
479  nb, & ! number of levels in zb
480  Rad, & ! radius of curvature of earth at observation
481  lat, & ! latitude at observation
482  und, & ! geoid undulation above WGS-84
483  refrac, & ! refractivity of model on model levels
484  dnr_dref) ! Calculated gradient of nr
485 
486 
487 USE ufo_gnssro_ukmo1d_utils_mod, only: ops_gpsro_geop_geom
488 
489 IMPLICIT NONE
490 
491 ! Subroutine arguments:
492 INTEGER, INTENT(IN) :: nb ! number of levels in zb
493 REAL(kind_real), INTENT(IN) :: zb(nb) ! geopotential heights on zb levels /m
494 REAL(kind_real), INTENT(IN) :: rad ! local radius of curvature of earth /m
495 REAL(kind_real), INTENT(IN) :: lat ! latitude at observation/ degrees
496 REAL(kind_real), INTENT(IN) :: und ! geoid undulation
497 REAL(kind_real), INTENT(IN) :: refrac(nb) ! refractivity on model levels / N
498 REAL(kind_real), INTENT(OUT) :: dnr_dref(nb,nb) ! Calculated gradient of nr wrt ref
499 
500 ! Local declarations:
501 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_nrK"
502 REAL(kind_real) :: r(nb) ! radius of model levels /m
503 REAL(kind_real) :: z(nb) ! geopotential heights on zb levels /m, local copy of zb
504 INTEGER :: i
505 
506 !Initialise the matrix
507 dnr_dref = 0.0
508 
509 !----------------------------------------------
510 ! 1. Convert zb values to geometric altitudes
511 !---------------------------------------------
512 z= zb + und ! approx. convert to geopotential above WGS-84 ellipsoid
513 CALL ops_gpsro_geop_geom (lat, &
514  z)
515 
516 !--------------------------------------------------
517 ! 2. Calculate dnr/dref
518 !--------------------------------------------------
519 
520 r = rad + z
521 
522 DO i = 1,nb
523  IF (zb(i) > 0.0 .AND. refrac(i) > 0.0) THEN
524  dnr_dref(i,i) = 1.0e-6 * r(i)
525  END IF
526 END DO
527 
528 END SUBROUTINE ops_gpsrocalc_nrk
529 
530 !-------------------------------------------------------------------------------
531 ! GPSRO 1D bending angle operator K code.
532 !-------------------------------------------------------------------------------
533 SUBROUTINE ops_gpsrocalc_alphak (nobs, &
534  nlev, &
535  a, &
536  refrac, &
537  nr, &
538  Kmat_ref, &
539  Kmat_nr)
540 
541 IMPLICIT NONE
542 
543 ! Subroutine arguments:
544 INTEGER, INTENT(IN) :: nobs ! size of ob. vector
545 INTEGER, INTENT(IN) :: nlev ! no. of refractivity levels
546 REAL(kind_real), INTENT(IN) :: a(nobs) ! observation impact parameters
547 REAL(kind_real), INTENT(IN) :: refrac(nlev) ! refractivity values on model levels
548 REAL(kind_real), INTENT(IN) :: nr(nlev) ! refractive index * radius product
549 REAL(kind_real), INTENT(OUT) :: kmat_ref(nobs,nlev) ! BA gradient wrt refractivity
550 REAL(kind_real), INTENT(OUT) :: kmat_nr(nobs,nlev) ! BA gradient wrt index * radius product
551 
552 ! Local declarations:
553 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSROcalc_alphaK"
554 INTEGER :: i
555 INTEGER :: n
556 INTEGER :: ibot
557 INTEGER :: jbot
558 INTEGER :: kbot
559 REAL(kind_real) :: kval(nlev - 1) ! exponential decay rate between levels
560 REAL(kind_real) :: root_2pia
561 REAL(kind_real) :: ref_low ! refractivity at lower level
562 REAL(kind_real) :: nr_low ! impact parameter lower level
563 REAL(kind_real) :: tup
564 REAL(kind_real) :: tlow ! upper/lower bounds to error function
565 REAL(kind_real) :: dalpha ! delta bending angle
566 REAL(kind_real) :: diff_erf ! error function result
567 REAL(kind_real) :: erf_tup
568 REAL(kind_real) :: erf_tlow ! error function of tup and tlow
569 REAL(kind_real) :: t ! int. step in approx. to error function
570 REAL(kind_real), PARAMETER :: a1 = 0.3480242 ! consts used in error function approx.
571 REAL(kind_real), PARAMETER :: a2 = -0.0958798
572 REAL(kind_real), PARAMETER :: a3 = 0.7478556
573 REAL(kind_real), PARAMETER :: p = 0.47047
574 REAL(kind_real) :: dkval_dref(nlev - 1,2)
575 REAL(kind_real) :: dkval_dnr(nlev - 1,2)
576 REAL(kind_real) :: dalpha_dref(2)
577 REAL(kind_real) :: dalpha_dnr(2)
578 REAL(kind_real) :: dalpha_dk
579 REAL(kind_real) :: dalpha_drlow
580 REAL(kind_real) :: dalpha_derf
581 REAL(kind_real) :: dalpha_dnrlow
582 REAL(kind_real) :: dnrlow_dref(2)
583 REAL(kind_real) :: dnrlow_dnr(2)
584 REAL(kind_real) :: drlow_dref(2)
585 REAL(kind_real) :: drlow_dnr(2)
586 REAL(kind_real) :: drlow_dk
587 REAL(kind_real) :: derf_dref(2)
588 REAL(kind_real) :: derf_dnr(2)
589 REAL(kind_real) :: derf_dtup
590 REAL(kind_real) :: derf_dtlow
591 REAL(kind_real) :: dtup_dnr(2)
592 REAL(kind_real) :: dtlow_dnr(2)
593 REAL(kind_real) :: dtup_dref(2)
594 REAL(kind_real) :: dtlow_dref(2)
595 REAL(kind_real) :: dtup_dk
596 REAL(kind_real) :: dtlow_dk
597 
598 !-------------------------------------------------------------------------------
599 ! Initialise the K matrices
600 !-------------------------------------------------------------------------------
601 
602 kmat_ref(:,:) = 0.0
603 kmat_nr(:,:) = 0.0
604 dkval_dref(:,:) = 0.0
605 dkval_dnr(:,:) = 0.0
606 
607 jbot = 1
608 
609 DO
610 
611  IF (refrac(jbot) > 0.0 .AND. nr(jbot) > 0.0) EXIT
612 
613  jbot = jbot + 1
614 
615 END DO
616 
617 !-------------------------------------------------------------------------------
618 ! Calculate lowest usable level (because of superrefraction)
619 !-------------------------------------------------------------------------------
620 
621 kbot = nlev
622 
623 DO i = nlev,jbot + 1,-1
624 
625  ! to avoid large gradients
626  IF ((nr(kbot) - nr(kbot-1)) < 10.0) EXIT
627 
628  kbot = kbot - 1
629 
630 END DO
631 
632 jbot = max(jbot,kbot)
633 
634 !-------------------------------------------------------------------------------
635 ! Calculate the exponential decay rate between levels
636 !-------------------------------------------------------------------------------
637 
638 DO i = jbot,nlev - 1
639 
640  kval(i) = log(refrac(i) / refrac(i + 1)) / &
641  max(1.0,(nr(i + 1) - nr(i)))
642 
643  kval(i) = max(1.0e-6,kval(i))
644 
645  IF (kval(i) > 1.0e-6) THEN
646 
647  dkval_dref(i,1) = 1.0 / (refrac(i) * max(1.0,(nr(i + 1) - nr(i))))
648  dkval_dref(i,2) = -1.0 / (refrac(i + 1) * max(1.0,(nr(i + 1) - nr(i))))
649 
650  dkval_dnr(i,1) = kval(i) / max(1.0,(nr(i + 1) - nr(i)))
651  dkval_dnr(i,2) = -kval(i) / max(1.0,(nr(i + 1) - nr(i)))
652 
653  END IF
654 
655 END DO
656 
657 !-------------------------------------------------------------------------------
658 ! Calculate the bending angle gradients
659 !-------------------------------------------------------------------------------
660 
661 DO n = 1,nobs
662 
663  IF (a(n) < nr(jbot) .OR. a(n) > nr(nlev)) cycle
664 
665  root_2pia = sqrt(2.0 * pi * a(n))
666 
667  ibot = jbot
668 
669  ! Find bottom state vector level
670  !----------------------------------
671  DO
672 
673  ! check more than 1 metre apart to stop large gradients in K code
674  ! ---------------------------------------------------------------
675  IF (((nr(ibot + 1) - a(n)) > 1.0) .OR. ibot == nlev - 1) EXIT
676 
677  ibot = ibot + 1
678 
679  END DO
680 
681  tlow = 0.0
682 
683  DO i = ibot, nlev - 1
684 
685  ! initialise matrices
686  !---------------------
687 
688  dalpha_dref(:) = 0.0
689  dalpha_dnr(:) = 0.0
690  drlow_dref(:) = 0.0
691  drlow_dnr(:) = 0.0
692  dnrlow_dref(:) = 0.0
693  dnrlow_dnr(:) = 0.0
694  dtup_dnr(:) = 0.0
695  dtup_dref(:) = 0.0
696  dtlow_dnr(:) = 0.0
697  dtlow_dref(:) = 0.0
698  derf_dtup = 0.0
699  derf_dtlow = 0.0
700  dalpha_drlow = 0.0
701  dalpha_dnrlow = 0.0
702  dalpha_dk = 0.0
703  dalpha_derf = 0.0
704  drlow_dk = 0.0
705  dtup_dk = 0.0
706  dtlow_dk = 0.0
707 
708  ! Values of refractivity and impact parameter at lower level
709  !-----------------------------------------------------------
710  IF (i == ibot) THEN
711 
712  ref_low = refrac(i) * exp(-kval(i) * (a(n) - nr(i)))
713 
714  drlow_dref(1)= ref_low / refrac(i)
715  drlow_dk = -ref_low * (a(n) - nr(i))
716  drlow_dnr(1) = ref_low * kval(i)
717 
718  nr_low = a(n)
719 
720  dnrlow_dnr(1) = 0.0
721 
722  ELSE
723 
724  ref_low = refrac(i)
725 
726  drlow_dref(1) = 1.0
727 
728  nr_low = nr(i)
729 
730  dnrlow_dnr(1) = 1.0
731 
732  END IF
733 
734  drlow_dref(:) = drlow_dref(:) + drlow_dk * dkval_dref(i,:)
735  drlow_dnr(:) = drlow_dnr(:) + drlow_dk * dkval_dnr(i,:)
736 
737 
738  ! Limits used in the error function
739  !----------------------------------
740  IF (i == nlev - 1) THEN
741 
742  ! simple extrapolation 100km above the uppermost level.
743  !-----------------------------------------------------
744  tup = sqrt(kval(i) * (nr(i + 1) + 1.0e5 - a(n)))
745 
746  dtup_dk = 0.5 * (nr(i + 1) + 1.0e5 - a(n)) / tup
747  dtup_dnr(2) = 0.5 * kval(i) / tup
748 
749  ELSE
750 
751  tup = sqrt(kval(i) * (nr(i + 1) - a(n)))
752 
753  dtup_dk = 0.5 * (nr(i + 1) - a(n)) / tup
754  dtup_dnr(2) = 0.5 * kval(i) / tup
755 
756  END IF
757 
758  dtup_dref(:) = dtup_dref(:) + dtup_dk * dkval_dref(i,:)
759  dtup_dnr(:) = dtup_dnr(:) + dtup_dk * dkval_dnr(i,:)
760 
761  tlow = 0.0
762 
763  IF (i > ibot) THEN
764  tlow = sqrt(kval(i) * (nr(i) - a(n)))
765  dtlow_dk = 0.5 * (nr(i) - a(n)) / tlow
766  dtlow_dnr(1) = 0.5 * kval(i) / tlow
767  END IF
768 
769  dtlow_dref(:) = dtlow_dref(:) + dtlow_dk * dkval_dref(i,:)
770  dtlow_dnr(:) = dtlow_dnr(:) + dtlow_dk * dkval_dnr(i,:)
771 
772  ! Abramowitz and Stegun approx. to error function
773  !------------------------------------------------
774  t = 1.0 / (1.0 + p * tup)
775  erf_tup = 1.0 - exp(-(tup ** 2)) * (a1 + (a2 + a3 * t) * t) * t
776  derf_dtup = exp(-(tup ** 2)) * ((2.0 * tup * (a1 + (a2 + a3 * t) * t) * t) + &
777  ((a1 + (2.0 * a2 + 3.0 * a3 * t) * t) * p * t ** 2))
778 
779  t = 1.0 / (1.0 + p * tlow)
780  erf_tlow = 1.0 - exp(-(tlow ** 2)) * (a1 + (a2 + a3 * t) * t) * t
781 
782  ! Multiplied by -1.0 to account for the usage in derf_dref and derf_dnr
783  !---------------------------------------------------------------------
784  derf_dtlow = -1.0 * exp(-(tlow ** 2)) * ((2.0 * tlow * (a1 + (a2 + a3 * t) * t) * t) + &
785  ((a1 + (2.0 * a2 + 3.0 * a3 * t) * t) * p * t ** 2))
786 
787  diff_erf = erf_tup - erf_tlow
788 
789  derf_dref(:) = derf_dtup * dtup_dref(:) + &
790  derf_dtlow * dtlow_dref(:)
791 
792  derf_dnr(:) = derf_dtup * dtup_dnr(:) + &
793  derf_dtlow * dtlow_dnr(:)
794 
795  dalpha = 1.0e-6 * root_2pia * sqrt(kval(i)) * &
796  ref_low * exp(kval(i) * (nr_low - a(n))) * diff_erf
797 
798  dalpha_drlow = dalpha / max(1.0e-10,ref_low)
799 
800  dalpha_derf = dalpha / max(1.0e-10,diff_erf)
801 
802  dalpha_dnrlow = dalpha * kval(i)
803 
804  dalpha_dk = dalpha * (nr_low - a(n) + 0.5 / kval(i))
805 
806  ! Now apply chain rule
807  !----------------------
808 
809  dalpha_dref(:) = dalpha_dref(:) + &
810  dalpha_drlow * drlow_dref(:) + &
811  dalpha_derf * derf_dref(:) + &
812  dalpha_dnrlow * dnrlow_dref(:) + &
813  dalpha_dk * dkval_dref(i,:)
814 
815  dalpha_dnr(:) = dalpha_dnr(:) + &
816  dalpha_drlow * drlow_dnr(:) + &
817  dalpha_derf * derf_dnr(:) + &
818  dalpha_dnrlow * dnrlow_dnr(:) + &
819  dalpha_dk * dkval_dnr(i,:)
820 
821  ! Now update matrices
822  !---------------------
823 
824  kmat_ref(n,i) = kmat_ref(n,i) + dalpha_dref(1)
825  kmat_nr(n,i) = kmat_nr(n,i) + dalpha_dnr(1)
826 
827  kmat_ref(n,i + 1) = kmat_ref(n,i + 1) + dalpha_dref(2)
828  kmat_nr(n,i + 1) = kmat_nr(n,i + 1) + dalpha_dnr(2)
829 
830  END DO
831 
832 END DO
833 
834 END SUBROUTINE ops_gpsrocalc_alphak
835 
836 end module ufo_gnssro_bendmetoffice_tlad_utils_mod
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:540
ufo_gnssro_ukmo1d_utils_mod
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:7
ufo_gnssro_bendmetoffice_tlad_utils_mod::ops_gpsro_refrack
subroutine, public ops_gpsro_refrack(nstate, nlevP, nb, nlevq, za, zb, x, GPSRO_pseudo_ops, GPSRO_vert_interp_ops, dref_dP, dref_dq)
Definition: ufo_gnssro_bendmetoffice_tlad_utils_mod.F90:59
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:485
ufo_constants_mod
Definition: ufo_constants_mod.F90:2
ufo_gnssro_ukmo1d_utils_mod::ops_gpsro_geop_geom
subroutine, public ops_gpsro_geop_geom(lat, z)
Definition: ufo_gnssro_bendmetoffice_utils_mod.F90:275
ufo_gnssro_bendmetoffice_tlad_utils_mod
Definition: ufo_gnssro_bendmetoffice_tlad_utils_mod.F90:7
ufo_constants_mod::rd
real(kind_real), parameter, public rd
Definition: ufo_constants_mod.F90:12