UFO
ufo_gnssro_bendmetoffice_tlad_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 !> Fortran module for gnssro bending angle Met Office's tangent linear and adjoint
8 
9 module ufo_gnssro_bendmetoffice_tlad_mod
10 
11 use iso_c_binding
12 use kinds
13 use ufo_vars_mod
14 use ufo_geovals_mod
15 use ufo_geovals_mod_c, only: ufo_geovals_registry
16 use vert_interp_mod
17 use ufo_basis_tlad_mod, only: ufo_basis_tlad
18 use obsspace_mod
19 use gnssro_mod_conf
20 use missing_values_mod
21 use fckit_log_module, only : fckit_log
22 use ufo_gnssro_bendmetoffice_tlad_utils_mod, only: &
23  ops_gpsrocalc_alphak, ops_gpsrocalc_nrk
24 use ufo_gnssro_ukmo1d_utils_mod, only: ops_gpsrocalc_nr
25 use ufo_utils_refractivity_calculator, only: &
26  ufo_calculate_refractivity, ufo_refractivity_kmat
27 
28 
29 integer, parameter :: max_string=800
30 
31 !> Fortran derived type for gnssro trajectory
32 type, extends(ufo_basis_tlad) :: ufo_gnssro_bendmetoffice_tlad
33  private
34  logical :: vert_interp_ops
35  logical :: pseudo_ops
36  real(kind_real) :: min_temp_grad
37  integer :: nlevp, nlevq, nlocs, iflip
38  real(kind_real), allocatable :: k(:,:)
39  contains
40  procedure :: setup => ufo_gnssro_bendmetoffice_setup
41  procedure :: delete => ufo_gnssro_bendmetoffice_tlad_delete
42  procedure :: settraj => ufo_gnssro_bendmetoffice_tlad_settraj
43  procedure :: simobs_tl => ufo_gnssro_bendmetoffice_simobs_tl
44  procedure :: simobs_ad => ufo_gnssro_bendmetoffice_simobs_ad
45 end type ufo_gnssro_bendmetoffice_tlad
46 
47 contains
48 
49 ! ------------------------------------------------------------------------------
50 ! Get the optional settings for the forward model, and save them in the object
51 ! so that they can be used in the code.
52 ! ------------------------------------------------------------------------------
53 subroutine ufo_gnssro_bendmetoffice_setup(self, vert_interp_ops, pseudo_ops, min_temp_grad)
54 
55 implicit none
56 
57 class(ufo_gnssro_bendmetoffice_tlad), intent(inout) :: self
58 logical(c_bool), intent(in) :: vert_interp_ops
59 logical(c_bool), intent(in) :: pseudo_ops
60 real(c_float), intent(in) :: min_temp_grad
61 
62 self % vert_interp_ops = vert_interp_ops
63 self % pseudo_ops = pseudo_ops
64 self % min_temp_grad = min_temp_grad
65 
66 end subroutine ufo_gnssro_bendmetoffice_setup
67 
68 
69 ! ------------------------------------------------------------------------------
70 ! Calculate the K-matrix (Jacobian) for the observation. It is necessary to run
71 ! this routine before calling the TL or AD routines.
72 ! ------------------------------------------------------------------------------
73 subroutine ufo_gnssro_bendmetoffice_tlad_settraj(self, geovals, obss)
74 
75  implicit none
76 ! Subroutine arguments
77  class(ufo_gnssro_bendmetoffice_tlad), intent(inout) :: self ! The object that we use to save data in
78  type(ufo_geovals), intent(in) :: geovals ! The input geovals
79  type(c_ptr), value, intent(in) :: obss ! The input observations
80 
81 ! Local parameters
82  character(len=*), parameter :: myname_="ufo_gnssro_bendmetoffice_tlad_settraj"
83 
84 ! Local variables
85  character(max_string) :: err_msg ! Messages to be output to the user
86  type(ufo_geoval), pointer :: q ! The model geovals - specific humidity
87  type(ufo_geoval), pointer :: prs ! The model geovals - atmospheric pressure
88  type(ufo_geoval), pointer :: rho_heights ! The model geovals - heights of the pressure-levels
89  type(ufo_geoval), pointer :: theta_heights ! The model geovals - heights of the theta-levels (stores q)
90  integer :: iobs ! Loop variable, observation number
91 
92  real(kind_real), allocatable :: obsLat(:) ! Latitude of the observation
93  real(kind_real), allocatable :: impact_param(:) ! Impact parameter of the observation
94  real(kind_real), allocatable :: obsLocR(:) ! Earth's radius of curvature at the observation tangent point
95  real(kind_real), allocatable :: obsGeoid(:) ! Undulation - height of the geoid above the ellipsoid
96 
97  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_tlad_settraj: begin"
98  call fckit_log%info(err_msg)
99 
100 ! Make sure that any previous values of geovals don't get carried over
101  call self%delete()
102 
103 ! get model state variables from geovals
104  call ufo_geovals_get_var(geovals, var_q, q) ! specific humidity
105  call ufo_geovals_get_var(geovals, var_prsi, prs) ! pressure
106  call ufo_geovals_get_var(geovals, var_z, theta_heights) ! Geopotential height of the normal model levels
107  call ufo_geovals_get_var(geovals, var_zi, rho_heights) ! Geopotential height of the pressure levels
108 
109 ! Keep copy of dimensions
110  self % nlevp = prs % nval
111  self % nlevq = q % nval
112  self % nlocs = obsspace_get_nlocs(obss)
113 
114 ! Check that the pressure values are decreasing (highest pressure at level 1)
115  self%iflip = 0
116  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
117  self%iflip = 1
118  write(err_msg,'(a)') ' ufo_gnssro_bendmetoffice_tlad_settraj:'//new_line('a')// &
119  ' Model vertical height profile is in descending order,'//new_line('a')// &
120  ' The data will be flipped for processing'
121  call fckit_log%info(err_msg)
122  end if
123 
124 ! Get the meta-data from the observations
125  allocate(obslat(self%nlocs))
126  allocate(impact_param(self%nlocs))
127  allocate(obslocr(self%nlocs))
128  allocate(obsgeoid(self%nlocs))
129  call obsspace_get_db(obss, "MetaData", "latitude", obslat)
130  call obsspace_get_db(obss, "MetaData", "impact_parameter", impact_param)
131  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", obslocr)
132  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", obsgeoid)
133  ALLOCATE(self % K(1:self%nlocs, 1:prs%nval + q%nval))
134 
135 ! For each observation, calculate the K-matrix
136  obs_loop: do iobs = 1, self % nlocs
137  if (self%iflip == 1) then
138  CALL jacobian_interface(prs % nval, & ! Number of pressure levels
139  q % nval, & ! Number of specific humidity levels
140  rho_heights % vals(rho_heights%nval:1:-1, iobs), & ! Heights of the pressure levels
141  theta_heights % vals(theta_heights%nval:1:-1, iobs), & ! Heights of the specific humidity levels
142  q % vals(q%nval:1:-1, iobs), & ! Values of the specific humidity
143  prs % vals(prs%nval:1:-1, iobs), & ! Values of the pressure
144  self % pseudo_ops, & ! Whether to use pseudo-levels in the calculation
145  self % vert_interp_ops, & ! Whether to interpolate using log(pressure)
146  self % min_temp_grad, & ! Minimum allowed vertical temperature gradient
147  obslocr(iobs), & ! Local radius of curvature of the earth
148  obslat(iobs), & ! Latitude of the observation
149  obsgeoid(iobs), & ! Geoid undulation at the tangent point
150  1, & ! Number of observations in the profile
151  impact_param(iobs:iobs), & ! Impact parameter for this observation
152  self % K(iobs:iobs,1:prs%nval+q%nval)) ! K-matrix (Jacobian of the observation with respect to the inputs)
153  else
154  CALL jacobian_interface(prs % nval, & ! Number of pressure levels
155  q % nval, & ! Number of specific humidity levels
156  rho_heights % vals(:,iobs), & ! Heights of the pressure levels
157  theta_heights % vals(:,iobs), & ! Heights of the specific humidity levels
158  q % vals(:,iobs), & ! Values of the specific humidity
159  prs % vals(:,iobs), & ! Values of the pressure
160  self % pseudo_ops, & ! Whether to use pseudo-levels in the calculation
161  self % vert_interp_ops, & ! Whether to interpolate using log(pressure)
162  self % min_temp_grad, & ! Minimum allowed vertical temperature gradient
163  obslocr(iobs), & ! Local radius of curvature of the earth
164  obslat(iobs), & ! Latitude of the observation
165  obsgeoid(iobs), & ! Geoid undulation at the tangent point
166  1, & ! Number of observations in the profile
167  impact_param(iobs:iobs), & ! Impact parameter for this observation
168  self % K(iobs:iobs,1:prs%nval+q%nval)) ! K-matrix (Jacobian of the observation with respect to the inputs)
169  end if
170  end do obs_loop
171 
172 ! Note that this routine has been run.
173  self%ltraj = .true.
174 
175  deallocate(obslat)
176  deallocate(impact_param)
177  deallocate(obslocr)
178  deallocate(obsgeoid)
179 
180 end subroutine ufo_gnssro_bendmetoffice_tlad_settraj
181 
182 ! ------------------------------------------------------------------------------
183 ! Given and increment to the model state, calculate an increment to the
184 ! observation
185 ! ------------------------------------------------------------------------------
186 subroutine ufo_gnssro_bendmetoffice_simobs_tl(self, geovals, hofx, obss)
187 
188  implicit none
189 
190 ! Subroutine arguments
191  class(ufo_gnssro_bendmetoffice_tlad), intent(in) :: self ! Object which is being used to transfer information
192  type(ufo_geovals), intent(in) :: geovals ! Model perturbations
193  real(kind_real), intent(inout) :: hofx(:) ! Increment to the observations
194  type(c_ptr), value, intent(in) :: obss ! Input - the observations
195 
196 ! Local parameters
197  character(len=*), parameter :: myname_="ufo_gnssro_bendmetoffice_simobs_tl"
198 
199 ! Local variables
200  integer :: iobs ! Loop variable, observation number
201  integer :: nlocs ! Number of observations
202  character(max_string) :: err_msg ! Message to be output
203  type(ufo_geoval), pointer :: q_d ! Increment to the specific humidity
204  type(ufo_geoval), pointer :: prs_d ! Increment to the air pressure
205  real(kind_real), allocatable :: x_d(:) ! Increment to the complete state
206 
207  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs_tl: begin"
208  call fckit_log%info(err_msg)
209 
210 ! Check if trajectory was set
211  if (.not. self%ltraj) then
212  write(err_msg,*) myname_, ' trajectory wasnt set!'
213  call abor1_ftn(err_msg)
214  endif
215 
216 ! Check if nlocs is consistent in geovals & hofx
217  if (geovals%nlocs /= size(hofx)) then
218  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
219  call abor1_ftn(err_msg)
220  endif
221 
222 ! Get variables from geovals
223  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
224  call ufo_geovals_get_var(geovals, var_prsi, prs_d) ! pressure on rho levels
225 
226  nlocs = self % nlocs ! number of observations
227 
228  allocate(x_d(1:prs_d%nval+q_d%nval))
229 ! Loop through the obs, calculating the increment to the observation
230  obs_loop: do iobs = 1, nlocs ! order of loop doesn't matter
231 
232  x_d(1:prs_d%nval) = prs_d % vals(:,iobs)
233  x_d(prs_d%nval+1:prs_d%nval+q_d%nval) = q_d % vals(:,iobs)
234  hofx(iobs) = sum(self % K(iobs,:) * x_d)
235 
236  end do obs_loop
237 
238  deallocate(x_d)
239 
240  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs_tl: complete"
241  call fckit_log%info(err_msg)
242 
243  return
244 
245 end subroutine ufo_gnssro_bendmetoffice_simobs_tl
246 
247 ! ------------------------------------------------------------------------------
248 ! Given an increment to the observation, find the equivalent increment to the
249 ! model state
250 ! ------------------------------------------------------------------------------
251 subroutine ufo_gnssro_bendmetoffice_simobs_ad(self, geovals, hofx, obss)
252 
253  use typesizes, only: wp => eightbytereal
254 
255  implicit none
256 
257 ! Subroutine arguments
258  class(ufo_gnssro_bendmetoffice_tlad), intent(in) :: self ! Object which is being used to transfer information
259  type(ufo_geovals), intent(inout) :: geovals ! Calculated perturbations to model state
260  real(kind_real), intent(in) :: hofx(:) ! Increment to the observations
261  type(c_ptr), value, intent(in) :: obss ! Input - the observations
262 
263 ! Local parameters
264  character(len=*), parameter :: myname_="ufo_gnssro_bendmetoffice_simobs_ad"
265 
266 ! Local variables
267  real(c_double) :: missing ! Missing data values
268  type(ufo_geoval), pointer :: q_d ! Pointer to the specific humidity perturbations
269  type(ufo_geoval), pointer :: prs_d ! Pointer to the pressure perturbations
270  integer :: iobs ! Loop variable, observation number
271  real(kind_real), allocatable :: x_d(:) ! Perturbation to the full model state
272  character(max_string) :: err_msg ! Message to be output
273 
274  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs_ad: begin"
275  call fckit_log%info(err_msg)
276 
277 ! Check if trajectory was set
278  if (.not. self%ltraj) then
279  write(err_msg,*) myname_, ' trajectory wasnt set!'
280  call abor1_ftn(err_msg)
281  endif
282 
283 ! Check if nlocs is consistent in geovals & hofx
284  if (geovals%nlocs /= size(hofx)) then
285  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
286  call abor1_ftn(err_msg)
287  endif
288 
289 ! Get variables from geovals
290  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
291  call ufo_geovals_get_var(geovals, var_prsi, prs_d) ! pressure
292 
293  missing = missing_value(missing)
294  allocate(x_d(1:prs_d%nval + q_d%nval))
295 
296 ! Loop through the obs, calculating the increment to the model state
297  obs_loop: do iobs = 1, self % nlocs
298 
299  if (hofx(iobs) /= missing) then
300  x_d = self % K(iobs,:) * hofx(iobs)
301  prs_d % vals(:,iobs) = x_d(1:prs_d%nval)
302  q_d % vals(:,iobs) = x_d(prs_d%nval+1:prs_d%nval+q_d%nval)
303  end if
304 
305  end do obs_loop
306 
307  deallocate(x_d)
308 
309  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs_ad: complete"
310  call fckit_log%info(err_msg)
311 
312  return
313 
314 end subroutine ufo_gnssro_bendmetoffice_simobs_ad
315 
316 !-------------------------------------------------------------------------
317 ! Tidy up the variables that are used for passing information
318 !-------------------------------------------------------------------------
319 subroutine ufo_gnssro_bendmetoffice_tlad_delete(self)
320 
321  implicit none
322  class(ufo_gnssro_bendmetoffice_tlad), intent(inout) :: self
323  character(len=*), parameter :: myname_="ufo_gnssro_bendmetoffice_tlad_delete"
324 
325  self%nlocs = 0
326  self%nlevp = 0
327  self%nlevq = 0
328  if (allocated(self%K)) deallocate(self%K)
329  self%ltraj = .false.
330 
331 end subroutine ufo_gnssro_bendmetoffice_tlad_delete
332 
333 !-------------------------------------------------------------------------
334 ! Interface for calculating the K-matrix for calculating TL/AD
335 !-------------------------------------------------------------------------
336 SUBROUTINE jacobian_interface(nlevp, &
337  nlevq, &
338  za, &
339  zb, &
340  q, &
341  prs, &
342  pseudo_ops, &
343  vert_interp_ops, &
344  min_temp_grad, &
345  ro_rad_curv, &
346  latitude, &
347  ro_geoid_und, &
348  nobs, &
349  zobs, &
350  K)
351 
352 IMPLICIT NONE
353 
354 INTEGER, INTENT(IN) :: nlevp ! The number of model pressure levels
355 INTEGER, INTENT(IN) :: nlevq ! The number of model theta levels
356 REAL(kind_real), INTENT(IN) :: za(:) ! The geometric height of the model pressure levels
357 REAL(kind_real), INTENT(IN) :: zb(:) ! The geometric height of the model theta levels
358 REAL(kind_real), INTENT(IN) :: q(1:nlevq) ! The model values that are being perturbed
359 REAL(kind_real), INTENT(IN) :: prs(1:nlevp) ! The model values that are being perturbed
360 LOGICAL, INTENT(IN) :: pseudo_ops ! Whether to use pseudo levels in the calculation
361 LOGICAL, INTENT(IN) :: vert_interp_ops ! Whether to use exner for the vertical interpolation
362 REAL(kind_real), INTENT(IN) :: min_temp_grad ! The minimum allowed vertical temperature gradient
363 REAL(kind_real), INTENT(IN) :: ro_rad_curv ! The earth's radius of curvature at the ob location
364 REAL(kind_real), INTENT(IN) :: latitude ! The latitude of the ob location
365 REAL(kind_real), INTENT(IN) :: ro_geoid_und ! The geoid undulation at the ob location
366 INTEGER, INTENT(IN) :: nobs ! The number of observations in this column
367 REAL(kind_real), INTENT(IN) :: zobs(:) ! The impact parameters of the column of observations
368 REAL(kind_real), INTENT(INOUT) :: K(:,:) ! The calculated K matrix
369 !
370 ! Things that may need to be output, as they are used by the TL/AD calculation
371 !
372 REAL(kind_real), ALLOCATABLE :: model_heights(:) ! Heights of the pseudo levels
373 REAL(kind_real), ALLOCATABLE :: refractivity(:) ! Refractivity on the pseudo levels
374 INTEGER :: nRefLevels ! Number of pseudo levels
375 REAL(kind_real) :: T(1:nlevq) ! Temperature on model levels
376 REAL(kind_real), ALLOCATABLE :: nr(:) ! Model calculation of impact parameters
377 !
378 ! Local variables
379 !
380 INTEGER :: num_pseudo ! Number of levels, including pseudo levels
381 REAL(kind_real) :: x(1:nlevp+nlevq) ! state vector
382 LOGICAL :: BAErr ! Whether we encountered an error in calculating the refractivity
383 CHARACTER(LEN=200) :: err_msg ! Output message
384 
385 ! Set up the size of the state
386 x(1:nlevp) = prs
387 x(nlevp+1:nlevp+nlevq) = q
388 
389 baerr = .false.
390 
391 CALL ufo_calculate_refractivity (nlevp, &
392  nlevq, &
393  za, &
394  zb, &
395  prs, &
396  q, &
397  pseudo_ops, &
398  vert_interp_ops, &
399  min_temp_grad, &
400  baerr, &
401  nreflevels, &
402  refractivity, &
403  model_heights)
404 
405 ALLOCATE(nr(1:nreflevels))
406 
407 IF (.NOT. baerr) THEN
408  ! 2. Calculate the impact parameter (refractive index * radius) on refractivity levels
409  CALL ops_gpsrocalc_nr (nreflevels, & ! number of model+pseudo-levels
410  model_heights, & ! geopotential heights of pseudo levels
411  refractivity, & ! refractivity of model on model+pseudo levels
412  ro_rad_curv, & ! radius of curvature of earth at observation
413  latitude, & ! latitude at observation
414  ro_geoid_und, & ! geoid undulation above WGS-84
415  nr) ! Calculated model impact parameters
416 
417  ! Calculate the K-matrix (Jacobian)
418  CALL ops_gpsro_getk(nlevp, &
419  nreflevels, &
420  nlevq, &
421  za, &
422  zb, &
423  model_heights, &
424  pseudo_ops, &
425  vert_interp_ops, &
426  min_temp_grad, &
427  prs, &
428  q, &
429  ro_rad_curv, &
430  latitude, &
431  ro_geoid_und, &
432  refractivity, &
433  nobs, &
434  zobs, &
435  nr, &
436  k)
437 ELSE
438  k = 0
439  write(err_msg,*) "Error in refractivity calculation"
440  CALL fckit_log % warning(err_msg)
441 END IF
442 
443 DEALLOCATE(nr)
444 DEALLOCATE(refractivity)
445 DEALLOCATE(model_heights)
446 
447 END SUBROUTINE jacobian_interface
448 
449 
450 !-------------------------------------------------------------------------
451 ! Calculate the K-matrix (Jacobian)
452 !-------------------------------------------------------------------------
453 SUBROUTINE ops_gpsro_getk(nlevp, &
454  nRefLevels, &
455  nlevq, &
456  za, &
457  zb, &
458  model_heights, &
459  pseudo_ops, &
460  vert_interp_ops, &
461  min_temp_grad, &
462  pressure, &
463  humidity, &
464  ro_rad_curv, &
465  latitude, &
466  ro_geoid_und, &
467  ref_model, &
468  nobs, &
469  zobs, &
470  nr, &
471  K)
472 !
473 ! Return the K-matrix for calculating TL/AD
474 !
475  IMPLICIT NONE
476 
477  INTEGER, INTENT(IN) :: nlevp ! The number of model pressure levels
478  INTEGER, INTENT(IN) :: nRefLevels ! Number of refractivity levels
479  INTEGER, INTENT(IN) :: nlevq ! The number of model theta levels
480  REAL(kind_real), INTENT(IN) :: za(:) ! The geometric height of the model pressure levels
481  REAL(kind_real), INTENT(IN) :: zb(:) ! The geometric height of the model theta levels
482  REAL(kind_real), INTENT(IN) :: model_heights(:) ! The geometric height of the refractivity levels
483  LOGICAL, INTENT(IN) :: pseudo_ops ! Whether to use pseudo levels in the calculation
484  LOGICAL, INTENT(IN) :: vert_interp_ops ! Whether to use exner for the vertical interpolation
485  REAL(kind_real), INTENT(IN) :: min_temp_grad ! Minimum allowed vertical temperature gradient
486  REAL(kind_real), INTENT(IN) :: pressure(nlevp) ! Model pressure
487  REAL(kind_real), INTENT(IN) :: humidity(nlevq) ! Model specific humidity
488  REAL(kind_real), INTENT(IN) :: ro_rad_curv ! The earth's radius of curvature at the ob location
489  REAL(kind_real), INTENT(IN) :: latitude ! The latitude of the ob location
490  REAL(kind_real), INTENT(IN) :: ro_geoid_und ! The geoid undulation at the ob location
491  REAL(kind_real), INTENT(IN) :: ref_model(nRefLevels) ! Model refractivity on theta levels - returned from forward model
492  INTEGER, INTENT(IN) :: nobs ! The number of observations in this column
493  REAL(kind_real), INTENT(IN) :: zobs(:) ! The impact parameters of the column of observations
494  REAL(kind_real), INTENT(IN) :: nr(nRefLevels) ! The impact parameters of the model data
495  REAL(kind_real), INTENT(OUT) :: K(nobs,nlevp+nlevq) ! The calculated K matrix
496 
497  REAL(kind_real) :: m1(nobs, nRefLevels) ! Intermediate term in the K-matrix calculation
498  REAL(kind_real), ALLOCATABLE :: dref_dp(:, :) ! Partial derivative of refractivity wrt. pressure
499  REAL(kind_real), ALLOCATABLE :: dref_dq(:, :) ! Partial derivative of refractivity wrt. specific humidity
500  REAL(kind_real) :: dnr_dref(nRefLevels, nRefLevels) ! Partial derivative of impact parameter wrt. refractivity
501  REAL(kind_real) :: dalpha_dref(nobs, nRefLevels) ! Partial derivative of bending angle wrt. refractivity
502  REAL(kind_real) :: dalpha_dnr(nobs, nRefLevels) ! Partial derivative of bending angle wrt. impact parameter
503 
504  ! 1. Calculate the gradient of ref wrt p (on rho levels) and q (on theta levels)
505  CALL ufo_refractivity_kmat(nlevp, &
506  nlevq, &
507  nreflevels, &
508  za, &
509  zb, &
510  pressure, &
511  humidity, &
512  pseudo_ops, &
513  vert_interp_ops, &
514  min_temp_grad, &
515  dref_dp, & !out
516  dref_dq) !out
517 
518  ! 2. Calculate the gradient of nr wrt ref
519  CALL ops_gpsrocalc_nrk (model_heights, & ! geopotential heights of pseudo levels
520  nreflevels, & ! number of refractivity levels
521  ro_rad_curv, & ! radius of curvature of earth at observation
522  latitude, & ! latitude at observation
523  ro_geoid_und, & ! geoid undulation above WGS-84
524  ref_model, & ! refractivity of model on model levels
525  dnr_dref) ! out
526 
527  ! 3. Calculate the gradient of bending angle wrt ref and nr
528  CALL ops_gpsrocalc_alphak (nobs, & ! size of ob. vector
529  nreflevels, & ! no. of refractivity levels
530  zobs, & ! obs impact parameters
531  ref_model, & ! refractivity values on model levels
532  nr, & ! index * radius product
533  dalpha_dref, & ! out
534  dalpha_dnr) ! out
535 
536  ! Calculate overall gradient of bending angle wrt p and q
537  m1 = matmul(dalpha_dnr,dnr_dref)
538  k(1:nobs, 1:nlevp) = matmul(dalpha_dref,dref_dp) + matmul(m1,dref_dp) !P part
539  k(1:nobs, nlevp+1:nlevp+nlevq) = matmul(dalpha_dref,dref_dq) + matmul(m1,dref_dq) !q part
540 
541  IF (ALLOCATED(dref_dp)) DEALLOCATE(dref_dp)
542  IF (ALLOCATED(dref_dq)) DEALLOCATE(dref_dq)
543 
544 END SUBROUTINE ops_gpsro_getk
545 
546 end module ufo_gnssro_bendmetoffice_tlad_mod
gnssro_mod_conf
Definition: gnssro_mod_conf.F90:2
ufo_basis_tlad_mod
Definition: ufo_basis_tlad_mod.F90:6
ufo_geovals_mod_c
Definition: GeoVaLs.interface.F90:7
ufo_geovals_mod_c::ufo_geovals_registry
type(registry_t), public ufo_geovals_registry
Linked list interface - defines registry_t type.
Definition: GeoVaLs.interface.F90:27
ufo_geovals_mod
Definition: ufo_geovals_mod.F90:7
ufo_geovals_mod::ufo_geovals_get_var
subroutine, public ufo_geovals_get_var(self, varname, geoval)
Definition: ufo_geovals_mod.F90:206
ufo_gnssro_bendmetoffice_tlad_mod
Fortran module for gnssro bending angle Met Office's tangent linear and adjoint.
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:9
ufo_gnssro_bendmetoffice_tlad_mod::jacobian_interface
subroutine jacobian_interface(nlevp, nlevq, za, zb, q, prs, pseudo_ops, vert_interp_ops, min_temp_grad, ro_rad_curv, latitude, ro_geoid_und, nobs, zobs, K)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:351
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_tlad_settraj
subroutine ufo_gnssro_bendmetoffice_tlad_settraj(self, geovals, obss)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:74
ufo_gnssro_bendmetoffice_tlad_mod::max_string
integer, parameter max_string
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:29
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_simobs_ad
subroutine ufo_gnssro_bendmetoffice_simobs_ad(self, geovals, hofx, obss)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:252
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_tlad_delete
subroutine ufo_gnssro_bendmetoffice_tlad_delete(self)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:320
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_setup
subroutine ufo_gnssro_bendmetoffice_setup(self, vert_interp_ops, pseudo_ops, min_temp_grad)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:54
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_simobs_tl
subroutine ufo_gnssro_bendmetoffice_simobs_tl(self, geovals, hofx, obss)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:187
ufo_gnssro_bendmetoffice_tlad_mod::ops_gpsro_getk
subroutine ops_gpsro_getk(nlevp, nRefLevels, nlevq, za, zb, model_heights, pseudo_ops, vert_interp_ops, min_temp_grad, pressure, humidity, ro_rad_curv, latitude, ro_geoid_und, ref_model, nobs, zobs, nr, K)
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:472
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_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
ufo_vars_mod
Definition: ufo_variables_mod.F90:8
ufo_vars_mod::var_prsi
character(len=maxvarlen), parameter, public var_prsi
Definition: ufo_variables_mod.F90:26
ufo_vars_mod::var_q
character(len=maxvarlen), parameter, public var_q
Definition: ufo_variables_mod.F90:22
ufo_vars_mod::var_zi
character(len=maxvarlen), parameter, public var_zi
Definition: ufo_variables_mod.F90:31
ufo_vars_mod::var_z
character(len=maxvarlen), parameter, public var_z
Definition: ufo_variables_mod.F90:29
vert_interp_mod
Fortran module to perform linear interpolation.
Definition: vert_interp.F90:8
ufo_basis_tlad_mod::ufo_basis_tlad
Definition: ufo_basis_tlad_mod.F90:12
ufo_geovals_mod::ufo_geoval
type to hold interpolated field for one variable, one observation
Definition: ufo_geovals_mod.F90:40
ufo_geovals_mod::ufo_geovals
type to hold interpolated fields required by the obs operators
Definition: ufo_geovals_mod.F90:47
ufo_gnssro_bendmetoffice_tlad_mod::ufo_gnssro_bendmetoffice_tlad
Fortran derived type for gnssro trajectory.
Definition: ufo_gnssro_bendmetoffice_tlad_mod.F90:32