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