UFO
ufo_groundgnss_metoffice_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 ground based GNSS Met Office's tangent linear and adjoint
8 
9 module ufo_groundgnss_metoffice_tlad_mod
10 use iso_c_binding
11 
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 ufo_basis_tlad_mod, only: ufo_basis_tlad
17 use obsspace_mod
18 use missing_values_mod
19 use fckit_log_module, only : fckit_log
20 use ufo_utils_refractivity_calculator, only: &
21  ufo_calculate_refractivity, ufo_refractivity_kmat
22 
23 
24 
25 use ufo_constants_mod, only: &
26  rd, & ! Gas constant for dry air
27  grav, & ! Gravitational field strength
28  n_alpha, & ! Refractivity constant a
29  n_beta
30 
31 integer, parameter :: max_string=800
32 
33 !> Fortran derived type for groundgnss trajectory
34 type, extends(ufo_basis_tlad) :: ufo_groundgnss_metoffice_tlad
35  private
36 
37  integer :: nlevp, nlevq, nlocs, iflip
38  real(kind_real), allocatable :: k(:,:)
39  real(kind_real), allocatable :: dztd_dp(:,:)
40  real(kind_real), allocatable :: dztd_dq(:,:)
41  logical :: vert_interp_ops
42  logical :: pseudo_ops
43  real(kind_real) :: min_temp_grad
44  contains
45  procedure :: setup => ufo_groundgnss_metoffice_setup
46  procedure :: delete => ufo_groundgnss_metoffice_tlad_delete
47  procedure :: settraj => ufo_groundgnss_metoffice_tlad_settraj
48  procedure :: simobs_tl => ufo_groundgnss_metoffice_simobs_tl
49  procedure :: simobs_ad => ufo_groundgnss_metoffice_simobs_ad
50 end type ufo_groundgnss_metoffice_tlad
51 
52 contains
53 
54 ! ------------------------------------------------------------------------------
55 ! Get the optional settings for the forward model, and save them in the object
56 ! so that they can be used in the code.
57 ! ------------------------------------------------------------------------------
58 subroutine ufo_groundgnss_metoffice_setup(self, f_conf)
59 
60 use fckit_configuration_module, only: fckit_configuration
61 implicit none
62 class(ufo_groundgnss_metoffice_tlad), intent(inout) :: self
63 type(fckit_configuration), intent(in) :: f_conf
64 call f_conf%get_or_die("min_temp_grad", self % min_temp_grad)
65 
66 end subroutine ufo_groundgnss_metoffice_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_groundgnss_metoffice_tlad_settraj(self, geovals, obss)
74 
75  implicit none
76 ! Subroutine arguments
77  class(ufo_groundgnss_metoffice_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_groundgnss_metoffice_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 (q on theta)
90  integer :: nstate ! The size of the state vector
91  integer :: iobs ! Loop variable, observation number
92  integer :: nobs ! Number of observations
93  integer :: ilev ! Loop variable, vertical level number
94  real(kind_real), allocatable :: zStation(:) ! The station height
95  real(kind_real), allocatable :: pressure(:) ! Model background values of air pressure (monotonic order)
96  real(kind_real), allocatable :: humidity(:) ! Model background specific humidity (in pressure monotonic order)
97  real(kind_real), allocatable :: za(:) ! Model heights of rho levs (in pressure monotonic order)
98  real(kind_real), allocatable :: zb(:) ! Model heights of theta levs (in pressure monotonic order)
99 
100  integer, parameter :: max_string = 800
101  character(max_string) :: message ! General message for output
102 
103  write(err_msg,*) "TRACE: ufo_groundgnss_metoffice_tlad_settraj: begin"
104  call fckit_log%info(err_msg)
105 
106 ! Make sure that any previous values of geovals don't get carried over
107  call self%delete()
108 
109 ! get model state variables from geovals
110  call ufo_geovals_get_var(geovals, var_q, q) ! specific humidity
111  call ufo_geovals_get_var(geovals, var_prsi, prs) ! pressure
112  call ufo_geovals_get_var(geovals, var_z, theta_heights) ! Geopotential height of the normal model levels
113  call ufo_geovals_get_var(geovals, var_zi, rho_heights) ! Geopotential height of the pressure levels
114 
115 ! Keep copy of dimensions
116  self % nlevp = prs % nval
117  self % nlevq = q % nval
118  self % nlocs = obsspace_get_nlocs(obss)
119 
120 ! Check whether the pressure levels are in descending order
121  self%iflip = 0
122  IF (prs % vals(1,1)-prs % vals(self%nlevp,1) < 0.0) THEN
123  self%iflip = 1
124  WRITE(message, *) "Pressure is in ascending order. Reorder the variables in vertical direction"
125  CALL fckit_log % warning(message)
126  END IF
127 
128 
129 ! Get the meta-data from the observations
130  nobs = obsspace_get_nlocs(obss)
131  allocate(zstation(nobs))
132 
133  call obsspace_get_db(obss, "MetaData", "station_height", zstation)
134 
135  nstate = prs % nval + q % nval
136  ALLOCATE(self % K(1:self%nlocs, 1:nstate))
137  ALLOCATE(pressure(1:self%nlevp))
138  ALLOCATE(humidity(1:self%nlevq))
139  ALLOCATE(za(1:self%nlevp))
140  ALLOCATE(zb(1:self%nlevq))
141 
142 ! For each observation, calculate the K-matrix
143  obs_loop: do iobs = 1, self % nlocs
144  IF (self%iflip == 1) THEN
145  do ilev = 1, self%nlevp
146  pressure(ilev) = prs % vals(self%nlevp-ilev+1,iobs)
147  za(ilev) = rho_heights % vals(self%nlevp-ilev+1,iobs)
148  end do
149  do ilev = 1, self%nlevq
150  humidity(ilev) = q % vals(self%nlevq-ilev+1,iobs)
151  zb(ilev) = theta_heights % vals(self%nlevq-ilev+1,iobs)
152  end do
153  ELSE
154  pressure = prs % vals(:,iobs)
155  humidity = q % vals(:,iobs)
156  za = rho_heights % vals(:,iobs)
157  zb = theta_heights % vals(:,iobs)
158  END IF
159  CALL groundgnss_jacobian_interface(self % nlevp, & ! Number of pressure levels
160  self % nlevq, & ! Number of specific humidity levels
161  za(1:self%nlevp), & ! Heights of the pressure levels
162  zb(1:self%nlevq), & ! Heights of the specific humidity levels
163  humidity(1:self%nlevq), & ! Values of the specific humidity
164  pressure(1:self%nlevp), & ! Values of the pressure
165  zstation(iobs), & ! Station height
166  iobs, & ! Ob number
167  self % min_temp_grad, & ! Minimum temperature gradient allowed
168  self % K(:, 1:nstate)) ! K-matrix (Jacobian of the observation with respect to the inputs
169 
170  end do obs_loop
171 
172 ! Note that this routine has been run.
173  self%ltraj = .true.
174 
175  deallocate(zstation)
176 
177 end subroutine ufo_groundgnss_metoffice_tlad_settraj
178 
179 
180 ! ------------------------------------------------------------------------------
181 ! Given an increment to the model state, calculate an increment to the
182 ! observation
183 ! ------------------------------------------------------------------------------
184 subroutine ufo_groundgnss_metoffice_simobs_tl(self, geovals, hofx, obss)
185 
186  implicit none
187 
188 ! Subroutine arguments
189  class(ufo_groundgnss_metoffice_tlad), intent(in) :: self ! Object which is being used to transfer information
190  type(ufo_geovals), intent(in) :: geovals ! Model perturbations
191  real(kind_real), intent(inout) :: hofx(:) ! Increment to the observations
192  type(c_ptr), value, intent(in) :: obss ! Input - the observations
193 
194 ! Local parameters
195  character(len=*), parameter :: myname_="ufo_groundgnss_metoffice_simobs_tl"
196 
197 ! Local variables
198  integer :: iobs ! Loop variable, observation number
199  integer :: nlocs ! Number of observations
200  integer :: ilev ! Loop variable, pressure level number
201  integer :: iflip ! Index for vertical flip
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 :: pressure_d(:) ! Increment to the air pressure (monotonic order)
206  real(kind_real), allocatable :: humidity_d(:) ! Increment to the specific humidity (in pressure monotonic order)
207  real(kind_real), allocatable :: x_d(:) ! Increment to the complete state
208 
209  write(err_msg,*) "TRACE: ufo_groundgnss_metoffice_simobs_tl: begin"
210  call fckit_log%info(err_msg)
211 
212 ! Check if trajectory was set
213  if (.not. self%ltraj) then
214  write(err_msg,*) myname_, ' trajectory wasnt set!'
215  call abor1_ftn(err_msg)
216  endif
217 
218 ! Check if nlocs is consistent in geovals & hofx
219  if (geovals%nlocs /= size(hofx)) then
220  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
221  call abor1_ftn(err_msg)
222  endif
223 
224 ! Get variables from geovals
225  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
226  call ufo_geovals_get_var(geovals, var_prsi, prs_d) ! pressure on rho levels
227 
228  nlocs = self % nlocs ! number of observations
229 
230  allocate(x_d(1:prs_d%nval+q_d%nval))
231  allocate(pressure_d(1:self % nlevp))
232  allocate(humidity_d(1:self % nlevq))
233 
234  iflip = self%iflip
235 
236 ! Loop through the obs, calculating the increment to the observation
237  obs_loop: do iobs = 1, nlocs ! order of loop doesn't matter
238 
239  IF (iflip == 1) THEN
240  do ilev = 1, self % nlevp
241  pressure_d(ilev) = prs_d % vals(self % nlevp-ilev+1,iobs)
242  end do
243  do ilev = 1, self % nlevq
244  humidity_d(ilev) = q_d % vals(self % nlevq-ilev+1,iobs)
245  end do
246  ELSE
247  pressure_d(1:self % nlevp) = prs_d % vals(:,iobs)
248  humidity_d(1:self % nlevq) = q_d % vals(:,iobs)
249  END IF
250 
251  x_d(1:prs_d%nval) = pressure_d
252  x_d(prs_d%nval+1:prs_d%nval+q_d%nval) = humidity_d
253  hofx(iobs) = sum(self % K(iobs,:) * x_d)
254 
255  end do obs_loop
256 
257  deallocate(x_d)
258  deallocate(pressure_d)
259  deallocate(humidity_d)
260 
261  write(err_msg,*) "TRACE: ufo_groundgnss_metoffice_simobs_tl: complete"
262  call fckit_log%info(err_msg)
263 
264  return
265 
266 end subroutine ufo_groundgnss_metoffice_simobs_tl
267 
268 
269 
270 ! ------------------------------------------------------------------------------
271 ! Given an increment to the observation, find the equivalent increment to the
272 ! model state
273 ! ------------------------------------------------------------------------------
274 subroutine ufo_groundgnss_metoffice_simobs_ad(self, geovals, hofx, obss)
275 
276  use typesizes, only: wp => eightbytereal
277 
278  implicit none
279 
280 ! Subroutine arguments
281  class(ufo_groundgnss_metoffice_tlad), intent(in) :: self ! Object which is being used to transfer information
282  type(ufo_geovals), intent(inout) :: geovals ! Calculated perturbations to model state
283  real(kind_real), intent(in) :: hofx(:) ! Increment to the observations
284  type(c_ptr), value, intent(in) :: obss ! Input - the observations
285 
286 ! Local parameters
287  character(len=*), parameter :: myname_="ufo_groundgnss_metoffice_simobs_ad"
288 
289 ! Local variables
290  real(c_double) :: missing ! Missing data values
291  type(ufo_geoval), pointer :: q_d ! Pointer to the specific humidity perturbations
292  type(ufo_geoval), pointer :: prs_d ! Pointer to the pressure perturbations
293  integer :: iobs ! Loop variable, observation number
294  integer :: ilev ! Loop variable, pressure level number
295  integer :: iflip ! Index for vertical flip
296  real(kind_real), allocatable :: x_d(:) ! Perturbation to the full model state
297  real(kind_real), allocatable :: pressure_d(:) ! Perturbation to pressure (monotonic order)
298  real(kind_real), allocatable :: humidity_d(:) ! Perturbation to specific humidity (in pressure monotonic order)
299  character(max_string) :: err_msg ! Message to be output
300 
301  write(err_msg,*) "TRACE: ufo_groundgnss_metoffice_simobs_ad: begin"
302  call fckit_log%info(err_msg)
303 
304 ! Check if trajectory was set
305  if (.not. self%ltraj) then
306  write(err_msg,*) myname_, ' trajectory wasnt set!'
307  call abor1_ftn(err_msg)
308  endif
309 
310 ! Check if nlocs is consistent in geovals & hofx
311  if (geovals%nlocs /= size(hofx)) then
312  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
313  call abor1_ftn(err_msg)
314  endif
315 
316 ! Get variables from geovals
317  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
318  call ufo_geovals_get_var(geovals, var_prsi, prs_d) ! pressure
319 
320  missing = missing_value(missing)
321  allocate(x_d(1:prs_d%nval + q_d%nval))
322  allocate(pressure_d(1:prs_d%nval))
323  allocate(humidity_d(1:q_d%nval))
324 
325  iflip = self%iflip
326 
327 ! Loop through the obs, calculating the increment to the model state
328  obs_loop: do iobs = 1, self % nlocs
329 
330  if (hofx(iobs) /= missing) then
331  x_d = self % K(iobs,:) * hofx(iobs)
332  pressure_d(1:prs_d%nval) = x_d(1:prs_d%nval)
333  humidity_d(1:q_d%nval) = x_d(prs_d%nval+1:prs_d%nval+q_d%nval)
334  end if
335 
336  if (iflip == 1) then
337  do ilev = 1, self % nlevp
338  prs_d % vals(self % nlevp-ilev+1,iobs) = pressure_d(ilev)
339  end do
340  do ilev = 1, self % nlevq
341  q_d % vals(self % nlevq-ilev+1,iobs) = humidity_d(ilev)
342  end do
343  else
344  prs_d % vals(:,iobs) = pressure_d(1:self % nlevp)
345  q_d % vals(:,iobs) = humidity_d(1:self % nlevq)
346  end if
347 
348 
349  end do obs_loop
350 
351  deallocate(x_d)
352  deallocate(pressure_d)
353  deallocate(humidity_d)
354 
355  write(err_msg,*) "TRACE: ufo_groundgnss_metoffice_simobs_ad: complete"
356  call fckit_log%info(err_msg)
357 
358  return
359 
360 end subroutine ufo_groundgnss_metoffice_simobs_ad
361 
362 
363 
364 !-------------------------------------------------------------------------
365 ! Tidy up the variables that are used for passing information
366 !-------------------------------------------------------------------------
367 subroutine ufo_groundgnss_metoffice_tlad_delete(self)
368 
369  implicit none
370  class(ufo_groundgnss_metoffice_tlad), intent(inout) :: self
371  character(len=*), parameter :: myname_="ufo_groundgnss_metoffice_tlad_delete"
372 
373  self%nlocs = 0
374  self%nlevp = 0
375  self%nlevq = 0
376  if (allocated(self%K)) deallocate(self%K)
377  self%ltraj = .false.
378 
379 end subroutine ufo_groundgnss_metoffice_tlad_delete
380 
381 
382 !-------------------------------------------------------------------------
383 ! Interface for calculating the K-matrix for calculating TL/AD
384 !-------------------------------------------------------------------------
385 SUBROUTINE groundgnss_jacobian_interface(nlevp, &
386  nlevq, &
387  za, &
388  zb, &
389  q, &
390  prs, &
391  zStation, &
392  iobs, &
393  gbgnss_min_temp_grad, &
394  K)
395 
396 IMPLICIT NONE
397 
398 INTEGER, INTENT(IN) :: nlevP ! The number of model pressure levels
399 INTEGER, INTENT(IN) :: nlevq ! The number of model theta levels
400 REAL(kind_real), INTENT(IN) :: za(:) ! The geometric height of the model pressure levels
401 REAL(kind_real), INTENT(IN) :: zb(:) ! The geometric height of the model theta levels
402 REAL(kind_real), INTENT(IN) :: q(1:nlevq) ! The model values that are being perturbed
403 REAL(kind_real), INTENT(IN) :: prs(1:nlevP) ! The model values that are being perturbed
404 REAL(kind_real), INTENT(IN) :: zStation ! Station height
405 REAL(kind_real), INTENT(IN) :: gbgnss_min_temp_grad ! The minimum temperature gradient which is used
406 INTEGER, INTENT(IN) :: iobs ! Ob number
407 
408 REAL(kind_real), INTENT(INOUT) :: K(:,:) ! The calculated K matrix
409 !
410 ! Things that may need to be output, as they are used by the TL/AD calculation
411 !
412 
413 REAL(kind_real) :: T(1:nlevq) ! Temperature on model levels
414 REAL(kind_real), ALLOCATABLE :: refrac(:) ! model refractivity on theta levels
415 !
416 ! Local variables
417 !
418 INTEGER :: nstate ! Number of levels in state vector
419 REAL(kind_real) :: x(1:nlevp+nlevq) ! state vector
420 LOGICAL :: refracerr ! Whether we encountered an error in calculating the refractivity
421 CHARACTER(LEN=200) :: err_msg ! Output message
422 
423 REAL(kind_real) :: pN(nlevq) ! Presure on theta levels
424 
425 REAL(kind_real), ALLOCATABLE :: model_heights(:) ! Geopotential heights of the refractivity levels (not needed for this oper)
426 INTEGER :: nRefLevels ! Number of levels in refractivity calculation
427 
428 ! Set up the size of the state
429 nstate = nlevp + nlevq
430 
431 refracerr = .false.
432 
433 ! Calculate the refractivity
434 CALL ufo_calculate_refractivity(nlevp, &
435  nlevq, &
436  za, &
437  zb, &
438  prs, &
439  q, &
440  .true., & ! vert_interp_ops
441  .false., & ! pseudo_ops
442  gbgnss_min_temp_grad, &
443  refracerr, &
444  nreflevels, &
445  refrac, &
446  model_heights)
447 
448 IF (.NOT. refracerr) THEN
449  ! Calculate the K-matrix (Jacobian)
450  CALL groundgnss_getk(nstate, &
451  nlevp, &
452  nlevq, &
453  za, &
454  zb, &
455  prs, &
456  q, &
457  zstation, &
458  iobs, &
459  gbgnss_min_temp_grad, &
460  refracerr, &
461  refrac, &
462  k)
463 ELSE
464  k = 0
465  write(err_msg,*) "Error in refractivity calculation"
466  CALL fckit_log % warning(err_msg)
467 END IF
468 
469 
470 END SUBROUTINE groundgnss_jacobian_interface
471 
472 
473 
474 !-------------------------------------------------------------------------
475 ! Calculate the K-matrix (Jacobian)
476 !-------------------------------------------------------------------------
477 SUBROUTINE groundgnss_getk(nstate, &
478  nlevP, &
479  nlevq, &
480  za, &
481  zb, &
482  P, &
483  q, &
484  zStation, &
485  iobs, &
486  gbgnss_min_temp_grad, &
487  refracerr, &
488  refrac, &
489  K)
490 
491 !
492 ! Return the K-matrix for calculating TL/AD
493 !
494 IMPLICIT NONE
495 
496 INTEGER, INTENT(IN) :: nstate
497 INTEGER, INTENT(IN) :: nlevP ! The number of model pressure levels
498 INTEGER, INTENT(IN) :: nlevq ! The number of model theta levels
499 REAL(kind_real), INTENT(IN) :: za(:) ! The geometric height of the model pressure levels
500 REAL(kind_real), INTENT(IN) :: zb(:) ! The geometric height of the model theta levels
501 REAL(kind_real), INTENT(IN) :: P(:) ! The model pressure values
502 REAL(kind_real), INTENT(IN) :: q(:) ! The model humidity values
503 REAL(kind_real), INTENT(IN) :: zStation ! Station height
504 INTEGER, INTENT(IN) :: iobs ! Ob number
505 REAL(kind_real), INTENT(IN) :: gbgnss_min_temp_grad ! The minimum temperature gradient which is used
506 LOGICAL, INTENT(INOUT) :: refracerr ! Whether we encountered an error in calculating the refractivity
507 REAL(kind_real), INTENT(IN) :: refrac(:) ! Model refractivity on theta levels - returned from forward model
508 REAL(kind_real), INTENT(INOUT) :: K(:,:) ! The calculated K matrix
509 
510 REAL(kind_real), ALLOCATABLE :: dref_dP(:, :) ! Partial derivative of refractivity wrt. pressure
511 REAL(kind_real), ALLOCATABLE :: dref_dq(:, :) ! Partial derivative of refractivity wrt. specific humidity
512 
513 
514 ! Local constants
515 
516 CHARACTER (LEN=*), PARAMETER :: RoutineName = "Groundgnss_GetK"
517 REAL, PARAMETER :: refrac_scale = 1.0e-6 ! Conversion factor between refractivity and refractive index
518 
519 ! Local variables
520 
521 INTEGER :: Level ! Used for iteration over levels
522 INTEGER :: Lowest_Level ! Lowest height level
523 INTEGER :: FirstNeg ! First negative
524 
525 REAL(kind_real) :: p_local(nlevp) ! pressure on rho levels (with no negative pressures)
526 REAL(kind_real) :: pN(nlevq) ! pressure on theta levels
527 REAL(kind_real) :: LocalZenithDelay ! Zenith Total Delay
528 REAL(kind_real) :: h_diff, station_diff ! Height diff, station height diff
529 REAL(kind_real) :: c_rep ! 1/scale height
530 REAL(kind_real) :: const, c, term1, term2 ! constant, scale height,
531 REAL(kind_real) :: StationRefrac ! Refraction at station height
532 REAL(kind_real) :: z_weight1 ! For linear interpolation
533 REAL(kind_real) :: z_weight2 ! For linear interpolation
534 REAL(kind_real) :: dc_dref ! derivative of scale height wrt refrac term1
535 REAL(kind_real) :: dc_dref2 ! derivative of scale height wrt refrac term2
536 REAL(kind_real) :: dztd_dc ! derivative of ZTD wrt scale height
537 REAL(kind_real) :: drefsta_dref ! derivative of station refrac wrt refrac
538 REAL(kind_real) :: drefsta_dc ! derivative of station refrac wrt scale height
539 REAL(kind_real) :: dztd_drefsta ! derivative of ZTD wrt station refrac
540 REAL(kind_real), ALLOCATABLE :: dp_local_dPin(:,:) ! derivative of pressure rho levels wrt pressure
541 REAL(kind_real), ALLOCATABLE :: dztd_dpN(:) ! array for derivative w.r.t top theta level
542 REAL(kind_real) :: dztd_dq(nlevq) ! The calculated K matrix
543 REAL(kind_real) :: dztd_dp(nlevP) ! The calculated K matrix
544 REAL(kind_real), ALLOCATABLE :: dPb_dP(:,:) ! derivative of pressure theta wrt pressure
545 REAL(kind_real), ALLOCATABLE :: dztd_dref(:) ! derivative of ZTD wrt refrac
546 REAL(kind_real), ALLOCATABLE :: x1(:,:) ! Matrix placeholder
547 REAL(kind_real), ALLOCATABLE :: x2(:) ! Matrix placeholder
548 
549 REAL(kind_real), PARAMETER :: hpa_to_pa = 100.0 ! hPa to Pascal conversion
550 REAL(kind_real), PARAMETER :: PressScale = 6000.0 ! Pressure scale height
551 CHARACTER(LEN=200) :: err_msg ! Output message
552 integer, parameter :: max_string = 800
553 character(max_string) :: message ! General message for output
554 
555 !----------------------------------------------------------------------
556 
557 !-------------------------------------------------------
558 ! 0. Initialise variables
559 !-------------------------------------------------------
560 
561 ALLOCATE (dref_dp(nlevq,nlevp))
562 ALLOCATE (dref_dq(nlevq,nlevq))
563 ALLOCATE (dpb_dp(nlevq,nlevp))
564 ALLOCATE (dztd_dref(nlevq))
565 ALLOCATE (dztd_dpn(nlevq))
566 ALLOCATE (x2(nlevp))
567 ALLOCATE (x1(nlevq,nlevp))
568 ALLOCATE (dp_local_dpin(nlevp,nlevp))
569 
570 ! Initialise matrices
571 
572 dref_dq(:,:) = 0.0
573 dref_dp(:,:) = 0.0
574 dztd_dref(:) = 0.0
575 dztd_dpn(:) = 0.0
576 dpb_dp(:,:) = 0.0
577 dztd_dq(:) = 0.0
578 dztd_dp(:) = 0.0
579 x1(:,:) = 0.0
580 x2(:) = 0.0
581 dp_local_dpin(:,:) = 0.0
582 p_local(:) = 0.0
583 pn(:) = 0.0
584 
585 localzenithdelay = 0.0
586 stationrefrac = 0.0
587 
588 ! If negative pressures exist, replace these
589 ! and any above with values calculated as an exponential decay (scale height
590 ! = 6km) from the highest positive pressure.
591 
592 firstneg = 0
593 DO level=1, nlevp
594  IF (level==1) THEN
595  p_local(level) = p(level)
596  dp_local_dpin(level,level) = 1.0
597  ELSE
598  IF ((p(level) <= 0.0 .OR. (firstneg /= 0 .AND. level > firstneg)) .OR. &
599  (p(level) > p(level-1))) THEN
600 
601  IF (firstneg == 0) firstneg = level
602 
603  p_local(level) = p(firstneg-1) * exp(-(za(level) - za(firstneg-1)) / pressscale)
604  dp_local_dpin(level,firstneg-1) = p_local(level) / p(firstneg-1)
605 
606  ELSE
607  p_local(level) = p(level)
608  dp_local_dpin(level,level) = 1.0
609  END IF
610  END IF
611 END DO
612 
613 ! Calculate Pressure on theta
614 ! Assume ln(p) linear with height
615 
616 DO level = 1, nlevp-1
617 
618  z_weight1 = (za(level+1) - zb(level)) / (za(level+1) - za(level))
619  z_weight2 = 1.0 - z_weight1
620 
621  pn(level) = exp(z_weight1 * log(p_local(level)) + z_weight2 * log(p_local(level+1)))
622 
623  dpb_dp(level,level) = pn(level) * z_weight1 / p_local(level)
624  dpb_dp(level,level+1) = pn(level) * z_weight2 / p_local(level+1)
625 
626 END DO
627 
628 ! Calculate the gradient of ref wrt p (on rho levels) and q (on theta levels)
629 
630 CALL ufo_refractivity_kmat (nlevp, &
631  nlevq, &
632  nlevq, &
633  za, &
634  zb, &
635  p, &
636  q, &
637  .false., & !pseudo_ops
638  .true., & ! vert_interp_ops
639  gbgnss_min_temp_grad, &
640  dref_dp, &
641  dref_dq, &
642  dpb_dp)
643 
644 
645 
646 ! In Layer where station height lies, define lowest level required for
647 ! iteration and integration
648 
649 DO level = 1, nlevp-1
650  IF (zb(level) > zstation) THEN
651  lowest_level = level
652  EXIT
653  END IF
654 END DO
655 
656 DO level = lowest_level, nlevq-1
657  IF (refrac(level+1) / refrac(level) <= 0.0) THEN
658 
659  write(err_msg,*) "Refractivity error. Refractivity < 0.0"
660  CALL fckit_log % warning(err_msg)
661 
662  RETURN
663 
664  END IF
665 END DO
666 
667 !---------------------------
668 ! 3. Calculate Zenith delays
669 !---------------------------
670 
671 ! Start at bottom level
672 
673 DO level = lowest_level, nlevq
674 
675  localzenithdelay = 0.0
676 
677  IF (level == lowest_level .AND. level /= 1) THEN
678 
679  ! If station lies above the lowest model level, interpolate refractivity
680  ! to station height
681  h_diff = zb(level - 1) - zb(level)
682  station_diff = zstation - zb(level)
683  c = (log(refrac(level) / refrac(level - 1))) / h_diff
684  stationrefrac = refrac(level - 1) * exp(-c * (zstation-zb(level - 1)))
685  const = (-stationrefrac / c) * exp(c * zstation)
686  term1 = exp(-c * (zb(level)))
687  term2 = exp(-c * zstation)
688  localzenithdelay = refrac_scale * const * (term1 - term2)
689 
690  c_rep = 1.0 / c
691  dztd_dc = (-refrac_scale * stationrefrac / c) * &
692  (c_rep + exp(c * station_diff) * (station_diff - c_rep))
693  dztd_dc = dztd_dc - (station_diff * localzenithdelay)
694  dc_dref = -1.0 / (refrac(level - 1) * h_diff)
695  dc_dref2 = 1.0 / (refrac(level) * h_diff)
696  drefsta_dref = stationrefrac / refrac(level - 1)
697  drefsta_dc = -(zstation - zb(level - 1)) * stationrefrac
698  drefsta_dref = drefsta_dref + drefsta_dc * dc_dref
699  dztd_drefsta = localzenithdelay / stationrefrac
700 
701  dztd_dref(level-1) = dztd_drefsta * drefsta_dref + dztd_dc * dc_dref
702  dztd_dref(level) = dztd_dc * dc_dref2
703 
704  ELSE IF (level == 1) THEN
705 
706  ! If station lies below model level 1 (ie. the lowest level for which refractivity is
707  ! calculated), then use c from the first full layer, but integrate down to height of
708  ! station
709 
710  h_diff = zb(level) - zb(level + 1)
711  c = (log(refrac(level + 1) / refrac(level))) / h_diff
712  const = (-refrac(level) / c) * exp(c * (zb(level)))
713  term1 = exp(-c * (zb(level + 1)))
714  term2 = exp(-c * zstation)
715  localzenithdelay = refrac_scale * const * (term1 - term2)
716 
717  c_rep = 1.0 / c
718  dztd_dc = (-refrac_scale * refrac(level) / c) * &
719  (c_rep + exp(c * h_diff) * (h_diff - c_rep))
720  dc_dref = -1.0 / (refrac(level) * h_diff)
721  dc_dref2 = 1.0 / (refrac(level + 1) * h_diff)
722 
723  dztd_dref(level) = localzenithdelay / refrac(level)
724  dztd_dref(level) = dztd_dref(level) + dztd_dc * dc_dref
725  dztd_dref(level+1)= dztd_dc * dc_dref2
726 
727  ELSE IF (level <= nlevq .AND. level > 2 .AND. level > lowest_level) THEN
728 
729  ! For the other Levels in the column above the station.
730 
731  h_diff = zb(level - 1) - zb(level)
732  c = (log(refrac(level) / refrac(level-1))) / h_diff
733  const = (-refrac(level - 1) / c) * exp(c * (zb(level - 1)))
734  term1 = exp(-c * (zb(level)))
735  term2 = exp(-c * (zb(level - 1)))
736  localzenithdelay = refrac_scale * const * (term1 - term2)
737 
738  c_rep = 1.0 / c
739  dztd_dc = (-refrac_scale * refrac(level - 1) / c) * (c_rep + exp(c * h_diff) * (h_diff - c_rep))
740  dc_dref = -1.0 / (refrac(level - 1) * h_diff)
741  dc_dref2 = 1.0 / (refrac(level) * h_diff)
742 
743  dztd_dref(level-1) = dztd_dref(level - 1) + localzenithdelay / refrac(level - 1) + dztd_dc * dc_dref
744  dztd_dref(level) = dztd_dc * dc_dref2
745 
746  END IF
747 
748 END DO
749 
750 !-------------------------------------------
751 ! Construct K Matrices
752 !-------------------------------------------
753 
754 ! Account for negative pressure
755 dref_dp = matmul(dref_dp,dp_local_dpin)
756 
757 dztd_dq(:) = matmul(dztd_dref,dref_dq)
758 dztd_dp(:) = matmul(dztd_dref,dref_dp)
759 
760 ! Account for negative pressure
761 dztd_dp(:) = matmul(dztd_dp,dp_local_dpin)
762 
763 ! First add in dZTD/dp for the top correction, which only depends on top level theta pressure
764 
765 dztd_dpn(nlevq) = refrac_scale * n_alpha * rd / (hpa_to_pa * grav)
766 x1 = matmul(dpb_dp, dp_local_dpin)
767 x2 = matmul(dztd_dpn, x1)
768 dztd_dp = x2 + dztd_dp
769 
770 k(iobs, 1:nlevp) = dztd_dp
771 k(iobs, nlevp+1:nstate) = dztd_dq
772 
773 DEALLOCATE (dp_local_dpin)
774 DEALLOCATE (x1)
775 DEALLOCATE (x2)
776 DEALLOCATE (dztd_dpn)
777 DEALLOCATE (dztd_dref)
778 DEALLOCATE (dpb_dp)
779 DEALLOCATE (dref_dq)
780 DEALLOCATE (dref_dp)
781 
782 END SUBROUTINE groundgnss_getk
783 
784 end module ufo_groundgnss_metoffice_tlad_mod
785 
ufo_basis_tlad_mod
Definition: ufo_basis_tlad_mod.F90:6
ufo_constants_mod
Definition: ufo_constants_mod.F90:2
ufo_constants_mod::rd
real(kind_real), parameter, public rd
Definition: ufo_constants_mod.F90:11
ufo_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_groundgnss_metoffice_tlad_mod
Fortran module for ground based GNSS Met Office's tangent linear and adjoint.
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:9
ufo_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_tlad_delete
subroutine ufo_groundgnss_metoffice_tlad_delete(self)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:368
ufo_groundgnss_metoffice_tlad_mod::groundgnss_getk
subroutine groundgnss_getk(nstate, nlevP, nlevq, za, zb, P, q, zStation, iobs, gbgnss_min_temp_grad, refracerr, refrac, K)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:490
ufo_groundgnss_metoffice_tlad_mod::max_string
integer, parameter max_string
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:31
ufo_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_tlad_settraj
subroutine ufo_groundgnss_metoffice_tlad_settraj(self, geovals, obss)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:74
ufo_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_simobs_ad
subroutine ufo_groundgnss_metoffice_simobs_ad(self, geovals, hofx, obss)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:275
ufo_groundgnss_metoffice_tlad_mod::groundgnss_jacobian_interface
subroutine groundgnss_jacobian_interface(nlevp, nlevq, za, zb, q, prs, zStation, iobs, gbgnss_min_temp_grad, K)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:395
ufo_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_simobs_tl
subroutine ufo_groundgnss_metoffice_simobs_tl(self, geovals, hofx, obss)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:185
ufo_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_setup
subroutine ufo_groundgnss_metoffice_setup(self, f_conf)
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:59
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
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_groundgnss_metoffice_tlad_mod::ufo_groundgnss_metoffice_tlad
Fortran derived type for groundgnss trajectory.
Definition: ufo_groundgnss_metoffice_tlad_mod.F90:34