UFO
ufo_gnssroonedvarcheck_rootsolv_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 
8 !> Fortran module for gnssro bending angle Met Office forward operator
9 
10 module ufo_gnssroonedvarcheck_rootsolv_mod
11 
12 use kinds, only: kind_real
13 use missing_values_mod, only: missing_value
14 use fckit_log_module, only: fckit_log
15 
16 private
17 public :: ops_gpsro_rootsolv_ba
18 
19 contains
20 
21 !-------------------------------------------------------------------------------
22 ! Solve the 1dvar problem
23 !-------------------------------------------------------------------------------
24 
25 SUBROUTINE ops_gpsro_rootsolv_ba (nstate, & ! size of state vector
26  nlevp, & ! no. of press. levels
27  nb, & ! no. of theta levels
28  nlevq, & ! no. of q levels
29  Nobs, & ! no of obs
30  za, & ! height of rho levels
31  zb, & ! height of theta levels
32  xb, & ! background vector
33  yobs, & ! ob. vector
34  zobs, & ! ob. impact parameters
35  Bsig, & ! standard dev. of b errors
36  Bm1, & ! Inverse Back. cov matrix
37  Om1, & ! Inverse Ob cov matrix
38  it, & ! no of iterations
39  x, & ! solution vector
40  yb, & ! obs first guess
41  ycalc, & ! obs at solution
42  J_pen, & ! penalty value
43  Amat, & ! solution cov. matrix
44  converged, & ! convergence flag
45  BAerr, & ! error flag
46  Do1DVar_error, & ! error flag
47  Iter_max, & ! Max number of iterations
48  Delta, & !
49  Delta_ct2, & !
50  GPSRO_pseudo_ops, & ! Whether to use pseudo levels
51  GPSRO_vert_interp_ops, & ! Whether to vertically interpolate using exner or ln(p)
52  gpsro_min_temp_grad, & ! Minimum vertical temperature gradient allowed
53  capsupersat, &
54  o_bdiff, & ! observed -background bending angle value
55  ro_rad_curv, & ! Radius of curvature of ellipsoid
56  latitude, & ! Latitude of occ
57  ro_geoid_und, & ! geoid undulation
58  tb, &
59  ts, &
60  dfs)
61 
62 
63 USE ufo_gnssro_ukmo1d_utils_mod, only: &
64  ops_gpsrocalc_alpha, &
65  ops_gpsrocalc_nr
66 
67 USE ufo_gnssro_bendmetoffice_tlad_utils_mod, only: &
68  ops_gpsrocalc_alphak, &
69  ops_gpsrocalc_nrk
70 
71 USE ufo_gnssroonedvarcheck_humidcheck_mod, only: &
72  ops_gpsro_humidcheck
73 
74 USE ufo_gnssroonedvarcheck_pen_mod, only: &
75  ops_gpsro_pen
76 
77 USE ufo_gnssroonedvarcheck_eval_derivs_mod, only: &
78  ops_gpsro_eval_derivs_ba
79 
80 USE ufo_utils_mod, only: &
81  invertmatrix, ops_cholesky
82 
83 use ufo_utils_refractivity_calculator, only: &
84  ufo_calculate_refractivity, &
85  ufo_refractivity_kmat
86 
87 IMPLICIT NONE
88 
89 ! Subroutine arguments:
90 INTEGER, INTENT(IN) :: nstate
91 INTEGER, INTENT(IN) :: nlevp
92 INTEGER, INTENT(IN) :: nb
93 INTEGER, INTENT(IN) :: nlevq
94 INTEGER, INTENT(IN) :: nobs
95 INTEGER, INTENT(IN) :: iter_max
96 REAL(kind_real), INTENT(IN) :: delta_ct2
97 REAL(kind_real), INTENT(IN) :: za(:)
98 REAL(kind_real), INTENT(IN) :: zb(:)
99 REAL(kind_real), INTENT(IN) :: xb(:)
100 REAL(kind_real), INTENT(IN) :: yobs(:)
101 REAL(kind_real), INTENT(IN) :: zobs(:)
102 REAL(kind_real), INTENT(IN) :: bsig(:)
103 REAL(kind_real), INTENT(IN) :: bm1(:,:)
104 REAL(kind_real), INTENT(IN) :: om1(:,:)
105 REAL(kind_real), INTENT(IN) :: delta
106 LOGICAL, INTENT(IN) :: gpsro_pseudo_ops
107 LOGICAL, INTENT(IN) :: gpsro_vert_interp_ops
108 REAL(kind_real), INTENT(IN) :: gpsro_min_temp_grad
109 LOGICAL, INTENT(IN) :: capsupersat
110 INTEGER, INTENT(OUT) :: it
111 REAL(kind_real), INTENT(OUT) :: x(:)
112 REAL(kind_real), INTENT(OUT) :: yb(:)
113 REAL(kind_real), INTENT(OUT) :: ycalc(:)
114 REAL(kind_real), INTENT(OUT) :: j_pen
115 REAL(kind_real), INTENT(OUT) :: amat(:,:)
116 LOGICAL, INTENT(OUT) :: converged
117 LOGICAL, INTENT(OUT) :: baerr
118 LOGICAL, INTENT(OUT) :: do1dvar_error
119 REAL(kind_real), INTENT(OUT) :: o_bdiff ! Measure of O-B for whole profile
120 REAL(kind_real), INTENT(IN) :: ro_rad_curv
121 REAL(kind_real), INTENT(IN) :: latitude
122 REAL(kind_real), INTENT(IN) :: ro_geoid_und
123 REAL(kind_real), INTENT(INOUT) :: tb(nlevq)
124 REAL(kind_real), INTENT(INOUT) :: ts(nlevq)
125 REAL(kind_real), INTENT(INOUT) :: dfs ! Measure of degrees of freesom of signal for whole profile
126 
127 ! Local declarations:
128 CHARACTER(len=*), PARAMETER :: routinename = "Ops_GPSRO_rootsolv_BA"
129 INTEGER, PARAMETER :: max_string = 800
130 
131 LOGICAL :: marq
132 INTEGER :: i
133 INTEGER :: it_marq
134 INTEGER :: errorcode
135 REAL(kind_real) :: j_old
136 REAL(kind_real) :: j_min
137 REAL(kind_real) :: pen_ob
138 REAL(kind_real) :: pen_back
139 REAL(kind_real) :: xold(nstate)
140 REAL(kind_real) :: x_temp(nstate)
141 REAL(kind_real) :: xmin(nstate)
142 REAL(kind_real) :: ymin(nobs)
143 REAL(kind_real) :: lambda
144 REAL(kind_real) :: lamp1
145 REAL(kind_real) :: d2j_dx2(nstate,nstate)
146 REAL(kind_real) :: dj_dx(nstate)
147 REAL(kind_real) :: diag_d2j(nstate)
148 REAL(kind_real) :: kmat(nobs,nstate)
149 REAL(kind_real) :: dx(nstate)
150 REAL(kind_real) :: conv_test
151 REAL(kind_real) :: ct2
152 REAL(kind_real) :: ct3
153 REAL(kind_real) :: d2 ! measure of step taken
154 REAL(kind_real) :: sdx(nstate)
155 REAL(kind_real) :: ok(nobs,nstate) ! O^-1 * Kmat
156 REAL(kind_real) :: kok(nstate,nstate) ! KT *O^-1 *K
157 REAL(kind_real) :: akok(nstate,nstate) ! Amat * above
158 REAL(kind_real), ALLOCATABLE :: nr(:)
159 REAL(kind_real), ALLOCATABLE :: dref_dp(:,:)
160 REAL(kind_real), ALLOCATABLE :: dref_dq(:,:)
161 REAL(kind_real), ALLOCATABLE :: dnr_dref(:,:)
162 REAL(kind_real), ALLOCATABLE :: dalpha_dref(:,:)
163 REAL(kind_real), ALLOCATABLE :: dalpha_dnr(:,:)
164 REAL(kind_real), ALLOCATABLE :: m1(:,:)
165 REAL(kind_real) :: t(nlevq)
166 REAL(kind_real) :: pressure(1:nlevp)
167 REAL(kind_real) :: humidity(1:nlevq)
168 
169 REAL(kind_real), ALLOCATABLE :: model_heights(:)
170 REAL(kind_real), ALLOCATABLE :: refractivity(:)
171 INTEGER :: nreflevels
172 CHARACTER(LEN=max_string) :: message
173 
174 !--------------
175 ! 1. Initialise
176 !--------------
177 
178 IF (gpsro_pseudo_ops) THEN
179  ALLOCATE(nr(2 * nlevq - 1))
180  ALLOCATE(dref_dp(2 * nlevq - 1,nlevp))
181  ALLOCATE(dref_dq(2 * nlevq - 1,nlevq))
182  ALLOCATE(dnr_dref(2 * nlevq - 1,2 * nlevq - 1))
183  ALLOCATE(dalpha_dref(nobs,2 * nlevq - 1))
184  ALLOCATE(dalpha_dnr(nobs,2 * nlevq - 1))
185  ALLOCATE(m1(nobs,2 * nlevq - 1))
186 ELSE
187  ALLOCATE(nr(nlevq))
188  ALLOCATE(dref_dp(nlevq,nlevp))
189  ALLOCATE(dref_dq(nlevq,nlevq))
190  ALLOCATE(dnr_dref(nlevq,nlevq))
191  ALLOCATE(dalpha_dref(nobs,nlevq))
192  ALLOCATE(dalpha_dnr(nobs,nlevq))
193  ALLOCATE(m1(nobs,nlevq))
194 END IF
195 
196 x(:) = xb(:) ! first guess = background
197 xold(:) = xb(:)
198 xmin(:) = xb(:)
199 
200 t(:) = missing_value(t(1))
201 tb(:) = missing_value(tb(1))
202 ts(:) = missing_value(ts(1))
203 yb(:) = missing_value(yb(1))
204 ycalc(:) = missing_value(ycalc(1))
205 
206 ! Logicals
207 
208 converged = .false.
209 do1dvar_error = .false.
210 marq = .false.
211 
212 ! Set convergence + penalty values to a large no.
213 
214 conv_test = 1.0e30
215 ct2 = 1.0e30
216 ct3 = 1.0e30
217 j_old = 1.0e30
218 j_pen = 1.0e30
219 j_min = 1.0e30
220 
221 ! Set lambda used in the Marqu. Lev. soln -initially a small value
222 
223 lambda = 1.0e-4_kind_real
224 
225 it = 0
226 it_marq = 0 ! no. of Marquardt iterations
227 
228 o_bdiff = missing_value(o_bdiff)
229 dfs = missing_value(dfs)
230 ycalc(:) = missing_value(ycalc(1))
231 ymin(:) = missing_value(ymin(1))
232 errorcode = missing_value(errorcode)
233 
234 !-----------------------
235 ! 2. Main iteration loop
236 !-----------------------
237 
238 ! Data to stdout on convergence of iteration loop
239 CALL fckit_log % info('J_pen|Conv_test|ct2|lambda|d2|(dJ/dx)^2|')
240 
241 iteration_loop: DO
242 
243  IF (it > iter_max .OR. &
244  (conv_test < delta .AND. &
245  ct2 < delta_ct2 * (nobs / 200.0_kind_real))) EXIT ! exit the iteration
246 
247  ! Count no. of iterations
248 
249  it = it + 1
250 
251  ! Calculate the bending angle profile for current guess x
252 
253  ! Call the 1D bending angle forward model
254 
255  ! 1. First calculate model refractivity on theta levels
256 
257  ! Unpack the solution to p and q, changing units
258  pressure = 100 * x(1:nlevp)
259  humidity = 0.001 * x(nlevp+1:nlevp+nlevq)
260 
261  CALL ufo_calculate_refractivity (nlevp, &
262  nlevq, &
263  za, &
264  zb, &
265  pressure, &
266  humidity, &
267  gpsro_pseudo_ops, &
268  gpsro_vert_interp_ops, &
269  gpsro_min_temp_grad, &
270  baerr, &
271  nreflevels, &
272  refractivity, &
273  model_heights, &
274  t)
275 
276  ! no point proceeding further if ...
277  IF (baerr) EXIT
278 
279  ! 2. Calculate the refractive index * radius on theta model levels (or model impact parameter)
280  CALL ops_gpsrocalc_nr (nreflevels, & ! number of refractivity levels
281  model_heights, & ! geopotential heights of refractivity levels
282  refractivity, & ! calculated refractivity
283  ro_rad_curv, & ! radius of curvature of earth at observation
284  latitude, & ! latitude at observation
285  ro_geoid_und, & ! geoid undulation above WGS-84
286  nr) ! Calculated model impact parameters
287 
288  ! 3. Calculate model bending angle on observation impact parameters
289  CALL ops_gpsrocalc_alpha (nobs, & ! size of ob. vector
290  nreflevels, & ! no. of refractivity levels
291  zobs, & ! obs impact parameters
292  refractivity, & ! refractivity values on model+pseudo levels
293  nr, & ! index * radius product
294  ycalc) ! forward modelled bending angle
295 
296  ! Store the bending angle values calculated with `x(:)=xb(:)'
297 
298  IF (it == 1) yb(:) = ycalc(:)
299 
300  ! For writing the background temperature to a file
301  IF (it == 1) tb(:) = t(:)
302 
303  ! Calculate the penalty (cost) function
304 
305  CALL ops_gpsro_pen (nstate, &
306  nobs, &
307  x, &
308  xb, &
309  yobs, &
310  ycalc, &
311  bm1, &
312  om1, &
313  pen_ob, & !out obs penalty
314  pen_back, & !out back penalty
315  j_pen)
316 
317  !store O-B value i.e. 2J/m on first iteration when dx=0
318  IF (it == 1) o_bdiff = 2.0_kind_real * j_pen / real(nobs)
319 
320  ! store the lowest cost value calculated
321 
322  IF (j_min > j_pen) THEN
323 
324  j_min = j_pen
325  xmin(:) = x(:)
326  ymin(:) = ycalc(:)
327 
328  END IF
329 
330  ! Convergence test
331 
332  conv_test = abs((j_old - j_pen) / min( j_old, real(nobs, kind=kind_real)))
333 
334  IF (it == 1 .OR. j_pen - j_old < 0.1_kind_real .OR. lambda > 1.0e6_kind_real) THEN
335 
336  ! Normal Newtonian iteration
337 
338  j_old = j_pen
339 
340  IF (.NOT. marq) lambda = 0.1_kind_real * lambda
341 
342  marq = .false.
343 
344  ! Evaluate the K matrix for current x
345  ! 1. Calculate the gradient of ref wrt p (on rho levels) and q (on theta levels)
346  ! Note: pressure and humidity were unpacked from x earlier in the loop
347 
348  CALL ufo_refractivity_kmat(nlevp, &
349  nlevq, &
350  nreflevels, &
351  za, &
352  zb, &
353  pressure, &
354  humidity, &
355  gpsro_pseudo_ops, &
356  gpsro_vert_interp_ops, &
357  gpsro_min_temp_grad, &
358  dref_dp, &
359  dref_dq)
360 
361  ! Change the units for the K-matrices
362  dref_dp(:,:) = 1.0e2 * dref_dp(:,:) ! hPa
363  dref_dq(:,:) = 1.0e-3 * dref_dq(:,:) ! g/kg
364 
365  ! 2. Calculate the gradient of nr wrt ref
366  CALL ops_gpsrocalc_nrk (model_heights, & ! geopotential heights of pseudo levels
367  nreflevels, & ! number of pseudo levels
368  ro_rad_curv, & ! radius of curvature of earth at observation
369  latitude, & ! latitude at observation
370  ro_geoid_und, & ! geoid undulation above WGS-84
371  refractivity, & ! refractivity of model on pseudo levels
372  dnr_dref) ! out
373 
374  ! 3. Calculate the gradient of bending angle wrt ref and nr
375  CALL ops_gpsrocalc_alphak (nobs, & ! size of ob. vector
376  nreflevels, & ! no. of refractivity pseudo levels
377  zobs, & ! obs impact parameters
378  refractivity, & ! refractivity values on pseudo levels
379  nr, & ! index * radius product
380  dalpha_dref, & ! out
381  dalpha_dnr) ! out
382 
383  ! Calculate overall gradient of bending angle wrt p and q
384 
385  m1 = matmul(dalpha_dnr,dnr_dref)
386  kmat(1:nobs,1:nlevp) = matmul(dalpha_dref,dref_dp) + matmul(m1,dref_dp) !P part
387  kmat(1:nobs,nlevp + 1:nstate) = matmul(dalpha_dref,dref_dq) + matmul(m1,dref_dq) !q part
388 
389  ! Store the state vector in xold
390 
391  xold(:) = x(:)
392 
393  ! Evaluate the -dJ_dx(vector NOTE SIGN!!) and d2J_dx2(matrix) at x
394 
395  CALL ops_gpsro_eval_derivs_ba (nstate, &
396  nobs, &
397  x, &
398  xb, &
399  yobs, &
400  ycalc, &
401  bm1, &
402  om1, &
403  kmat, &
404  dj_dx, &
405  d2j_dx2, &
406  diag_d2j)
407 
408  ! Store inverse of soln. cov matrix
409  amat(:,:) = d2j_dx2(:,:)
410 
411  ELSE
412 
413  ! M.Lev iteration.The previous increment increased the value
414  ! of the penalty function Use previous values of -dJ_dx, d2J_dx2 and xold
415 
416  marq = .true.
417  it_marq = it_marq + 1
418  lambda = 10.0_kind_real * lambda
419 
420  END IF
421 
422  ! Marq.Lev adjustment to diagonal terms
423 
424  lamp1 = lambda + 1.0_kind_real
425 
426  DO i = 1,nstate
427 
428  d2j_dx2(i,i) = diag_d2j(i)
429 
430  ! Marq.Lev. modification to diagonals of Hessian
431 
432  d2j_dx2(i,i) = d2j_dx2(i,i) * lamp1
433 
434  END DO
435 
436  ! Solve the matrix equation
437  !
438  ! d2J_dx2(:,:) . dx (:) = - dJ_dx (:)
439  !
440  ! to find dx(:) using the Cholesky decomposition routine
441  !
442 
443  CALL ops_cholesky (d2j_dx2, &
444  dj_dx, &
445  nstate, &
446  dx, & ! The answer, i.e. the "increment"
447  errorcode)
448 
449  IF (errorcode /= 0) EXIT
450 
451  ! Update estimate, but limit magnitude of increment with expected background
452  ! error
453 
454  x(:) = xold(:) + sign(min(abs(dx(:)), 2.0_kind_real * bsig(:), xold(:) / 2.0_kind_real), dx(:))
455 
456  ! Check values of humidity -limit to supersat and set <0.0 to = 0.0
457 
458  CALL ops_gpsro_humidcheck (nstate, &
459  nlevp, &
460  nlevq, &
461  za, &
462  zb, &
463  capsupersat, &
464  x)
465 
466  ! Second convergence criteria in terms of gradient
467 
468  ct2 = dot_product(abs(x(:) - xold(:)), abs(dj_dx(:)))
469  ct3 = dot_product(dj_dx(:), dj_dx(:))
470 
471  !-----------------------------------------
472  ! Save iteration info to standard output
473  !-----------------------------------------
474  sdx(:) = matmul(amat(:,:) , (x(:) - xold(:))) !S^-1.dx
475  d2 = dot_product((x(:) - xold(:)) , sdx(:)) !d^2=dx(S^-1)dx, size of step normalized by error size
476 
477  WRITE (message,'(6E14.6)') j_pen, conv_test, ct2, lambda, d2, ct3
478  CALL fckit_log % info(message)
479 
480 END DO iteration_loop
481 
482 ts(:) = t(:) !1DVAR solution temperature
483 
484 WRITE (message, '(A,I0)') 'Number of iterations ', it !write out number of iterations done
485 CALL fckit_log % info(message)
486 WRITE (message, '(A,F16.4)') 'O-B size ', o_bdiff
487 CALL fckit_log % info(message)
488 
489 ! Output the x(:) that gave the lowest cost function
490 
491 IF (j_min < j_pen) THEN
492 
493  j_pen = j_min
494  x(:) = xmin(:)
495  ycalc(:) = ymin(:)
496 
497 END IF
498 
499 IF (it <= iter_max) converged = .true.
500 
501 ! Calculate the solution error cov. matrix
502 
503 IF (errorcode == 0 .AND. it <= iter_max) THEN
504 
505  CALL invertmatrix (nstate, &
506  nstate, &
507  amat(:,:), &
508  errorcode)
509 
510  ! Check there isn't an error during the matrix inversion
511 
512  IF (errorcode /= 0) do1dvar_error = .true.
513 
514  !calculate degrees of freedom for signal
515  ok = matmul(om1(:,:) , kmat(:,:)) ! O^-1K
516  kok = matmul(transpose(kmat(:,:)) , ok(:,:)) ! KTO^-1K
517  akok = matmul(amat(:,:) , kok(:,:)) ! AKTO^-1K
518  dfs = 0.0 ! initialise DFS
519  !DFS is equal to trace of AKOK
520  DO i = 1, nstate
521  dfs = dfs + akok(i,i)
522  END DO
523 
524  WRITE (message,'(A,F16.4)') 'DFS', dfs
525  CALL fckit_log % info(message)
526 ELSE
527 
528  do1dvar_error = .true.
529 
530 END IF
531 
532 IF (ALLOCATED (nr)) DEALLOCATE (nr)
533 IF (ALLOCATED (dref_dp)) DEALLOCATE (dref_dp)
534 IF (ALLOCATED (dref_dq)) DEALLOCATE (dref_dq)
535 IF (ALLOCATED (dnr_dref)) DEALLOCATE (dnr_dref)
536 IF (ALLOCATED (dalpha_dref)) DEALLOCATE (dalpha_dref)
537 IF (ALLOCATED (dalpha_dnr)) DEALLOCATE (dalpha_dnr)
538 IF (ALLOCATED (m1)) DEALLOCATE (m1)
539 
540 END SUBROUTINE ops_gpsro_rootsolv_ba
541 
542 end module ufo_gnssroonedvarcheck_rootsolv_mod
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_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_gnssroonedvarcheck_eval_derivs_mod
Evaluate the 1st and 2nd deriv. of the cost function.
Definition: ufo_gnssroonedvarcheck_eval_derivs_mod.f90:10
ufo_gnssroonedvarcheck_eval_derivs_mod::ops_gpsro_eval_derivs_ba
subroutine, public ops_gpsro_eval_derivs_ba(Nstate, Nobs, x, xb, yobs, ycalc, BM1, OM1, Kmat, dJ_dx, d2J_dx2, diag_d2J)
Definition: ufo_gnssroonedvarcheck_eval_derivs_mod.f90:32
ufo_gnssroonedvarcheck_humidcheck_mod
Definition: ufo_gnssroonedvarcheck_humidcheck_mod.f90:8
ufo_gnssroonedvarcheck_humidcheck_mod::ops_gpsro_humidcheck
subroutine, public ops_gpsro_humidcheck(nstate, nlevp, nlevq, za, zb, capsupersat, x)
Definition: ufo_gnssroonedvarcheck_humidcheck_mod.f90:24
ufo_gnssroonedvarcheck_pen_mod
Calculates GPSRO penalty.
Definition: ufo_gnssroonedvarcheck_pen_mod.f90:10
ufo_gnssroonedvarcheck_pen_mod::ops_gpsro_pen
subroutine, public ops_gpsro_pen(Nstate, Nobs, x, xb, yobs, ycalc, BM1, OM1, pen_ob, pen_back, pen_func)
Definition: ufo_gnssroonedvarcheck_pen_mod.f90:31
ufo_gnssroonedvarcheck_rootsolv_mod
Fortran module for gnssro bending angle Met Office forward operator.
Definition: ufo_gnssroonedvarcheck_rootsolv_mod.f90:10
ufo_gnssroonedvarcheck_rootsolv_mod::ops_gpsro_rootsolv_ba
subroutine, public ops_gpsro_rootsolv_ba(nstate, nlevp, nb, nlevq, Nobs, za, zb, xb, yobs, zobs, Bsig, Bm1, Om1, it, x, yb, ycalc, J_pen, Amat, converged, BAerr, Do1DVar_error, Iter_max, Delta, Delta_ct2, GPSRO_pseudo_ops, GPSRO_vert_interp_ops, GPSRO_min_temp_grad, capsupersat, O_Bdiff, RO_Rad_Curv, Latitude, RO_geoid_und, Tb, Ts, DFS)
Definition: ufo_gnssroonedvarcheck_rootsolv_mod.f90:61
ufo_utils_mod
Fortran module with various useful routines.
Definition: ufo_utils_mod.F90:10
ufo_utils_mod::ops_cholesky
subroutine, public ops_cholesky(U, V, N, Q, ErrorCode)
Do cholesky decomposition.
Definition: ufo_utils_mod.F90:1110
ufo_utils_mod::invertmatrix
subroutine, public invertmatrix(n, m, a, status, matrix)
Definition: ufo_utils_mod.F90:1238
ufo_utils_refractivity_calculator
Module for containing a general refractivity forward operator and its K-matrix.
Definition: RefractivityCalculator.F90:17
ufo_utils_refractivity_calculator::ufo_calculate_refractivity
subroutine, public ufo_calculate_refractivity(nlevP, nlevq, za, zb, P, q, vert_interp_ops, pseudo_ops, min_temp_grad, refracerr, nRefLevels, refractivity, model_heights, temperature, interp_pressure)
Calculation of the refractivity from the pressure and specific humidity.
Definition: RefractivityCalculator.F90:81
ufo_utils_refractivity_calculator::ufo_refractivity_kmat
subroutine, public ufo_refractivity_kmat(nlevP, nlevq, nRefLevels, za, zb, P, q, pseudo_ops, vert_interp_ops, min_temp_grad, dref_dP, dref_dq, dPb_dP, refractivity)
Calculate general refractivity K matrix.
Definition: RefractivityCalculator.F90:314