UFO
ufo_gnssro_bendmetoffice_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_gnssro_bendmetoffice_mod
11 
12 use iso_c_binding
13 use kinds
14 use ufo_vars_mod
15 use ufo_geovals_mod
16 use ufo_geovals_mod_c, only: ufo_geovals_registry
17 use ufo_basis_mod, only: ufo_basis
18 use vert_interp_mod
19 use lag_interp_mod, only: lag_interp_const, lag_interp_smthweights
20 use obsspace_mod
21 use missing_values_mod
22 use ufo_gnssro_ukmo1d_utils_mod
23 use fckit_log_module, only : fckit_log
24 
25 implicit none
26 public :: ufo_gnssro_bendmetoffice
27 private
28 
29  !> Fortran derived type for gnssro trajectory
30 type, extends(ufo_basis) :: ufo_gnssro_bendmetoffice
31  logical :: vert_interp_ops
32  logical :: pseudo_ops
33  contains
34  procedure :: setup => ufo_gnssro_bendmetoffice_setup
35  procedure :: simobs => ufo_gnssro_bendmetoffice_simobs
36 end type ufo_gnssro_bendmetoffice
37 
38 contains
39 
40 ! ------------------------------------------------------------------------------
41 ! Get the optional settings for the forward model, and save them in the object
42 ! so that they can be used in the code.
43 ! ------------------------------------------------------------------------------
44 subroutine ufo_gnssro_bendmetoffice_setup(self, f_conf)
45 
46 use fckit_configuration_module, only: fckit_configuration
47 implicit none
48 class(ufo_gnssro_bendmetoffice), intent(inout) :: self
49 type(fckit_configuration), intent(in) :: f_conf
50 
51 call f_conf%get_or_die("vert_interp_ops", self % vert_interp_ops)
52 call f_conf%get_or_die("pseudo_ops", self % pseudo_ops)
53 
54 end subroutine ufo_gnssro_bendmetoffice_setup
55 
56 ! ------------------------------------------------------------------------------
57 ! 1-dimensional GNSS-RO forward operator for the Met Office system
58 ! ------------------------------------------------------------------------------
59 subroutine ufo_gnssro_bendmetoffice_simobs(self, geovals, hofx, obss)
60 
61  implicit none
62 
63  ! Arguments to this routine
64  class(ufo_gnssro_bendmetoffice), intent(in) :: self ! The object in which this operator is contained
65  type(ufo_geovals), intent(in) :: geovals ! The model values, interpolated to the obsevation locations
66  real(kind_real), intent(inout) :: hofx(:) ! The model forecast of the observations
67  type(c_ptr), value, intent(in) :: obss ! The observations, and meta-data for those observations
68 
69  character(len=*), parameter :: myname_ = "ufo_gnssro_bendmetoffice_simobs"
70  integer, parameter :: max_string = 800
71 
72  character(max_string) :: err_msg ! Error message for output
73  character(max_string) :: message ! General message for output
74  integer :: nobs ! Number of observations
75  integer :: ilev ! Loop variable, level number
76  integer :: iobs ! Loop variable, observation number
77  type(ufo_geoval), pointer :: q ! Model background values of specific humidity
78  type(ufo_geoval), pointer :: prs ! Model background values of air pressure
79  type(ufo_geoval), pointer :: theta_heights ! Model heights of levels containing specific humidity
80  type(ufo_geoval), pointer :: rho_heights ! Model heights of levels containing air pressure
81  real(kind_real), allocatable :: obsLat(:) ! Latitude of the observation
82  real(kind_real), allocatable :: obsLon(:) ! Longitude of the observation
83  real(kind_real), allocatable :: impact_param(:) ! Impact parameter of the observation
84  real(kind_real), allocatable :: radius_curv(:) ! Earth's radius of curvature at the observation tangent point
85  real(kind_real), allocatable :: undulation(:) ! Undulation - height of the geoid above the ellipsoid
86  logical :: flip_data ! Whether to reverse the order of the model data
87  logical :: BAErr ! Was there an error in the calculation?
88 
89  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs: begin"
90  call fckit_log%info(err_msg)
91 
92 ! check if nlocs is consistent in geovals & hofx
93  if (geovals%nlocs /= size(hofx)) then
94  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
95  call abor1_ftn(err_msg)
96  endif
97 
98  write(message, *) myname_, ' Running Met Office GNSS-RO forward operator with'
99  call fckit_log%info(message)
100  write(message, *) 'vert_interp_ops =', self % vert_interp_ops, &
101  'pseudo_ops =', self % pseudo_ops
102  call fckit_log%info(message)
103 
104 ! get variables from geovals
105  call ufo_geovals_get_var(geovals, var_q, q) ! specific humidity
106  call ufo_geovals_get_var(geovals, var_prsi, prs) ! pressure
107  call ufo_geovals_get_var(geovals, var_z, theta_heights) ! Geopotential height of the normal model levels
108  call ufo_geovals_get_var(geovals, var_zi, rho_heights) ! Geopotential height of the pressure levels
109 
110  write(message, '(A,10I6)') 'Q: ', q%nval, q%nlocs, shape(q%vals)
111  call fckit_log%info(message)
112  write(message, '(A,10I6)') 'Pressure: ', prs%nval, prs%nlocs, shape(prs%vals)
113  call fckit_log%info(message)
114 
115  nobs = obsspace_get_nlocs(obss)
116 
117  flip_data = .false.
118  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
119  write(err_msg,'(a)') ' ufo_gnssro_bendmetoffice_simobs:'//new_line('a')// &
120  ' Model vertical height profile is in descending order,'//new_line('a')// &
121  ' The data will be flipped for processing'
122  call fckit_log%info(err_msg)
123  flip_data = .true.
124  end if
125 
126 ! set obs space struture
127  allocate(obslon(nobs))
128  allocate(obslat(nobs))
129  allocate(impact_param(nobs))
130  allocate(radius_curv(nobs))
131  allocate(undulation(nobs))
132 
133  call obsspace_get_db(obss, "MetaData", "longitude", obslon)
134  call obsspace_get_db(obss, "MetaData", "latitude", obslat)
135  call obsspace_get_db(obss, "MetaData", "impact_parameter", impact_param)
136  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", radius_curv)
137  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", undulation)
138 
139  call fckit_log%info(err_msg)
140 
141  obs_loop: do iobs = 1, nobs
142 
143  call fckit_log%info(err_msg)
144 
145  if (flip_data) then
146  call ops_gpsro_forwardmodel(prs % nval, &
147  q % nval, &
148  rho_heights % vals(rho_heights%nval:1:-1, iobs), &
149  theta_heights % vals(theta_heights%nval:1:-1, iobs), &
150  prs % vals(prs%nval:1:-1, iobs), &
151  q % vals(q%nval:1:-1, iobs), &
152  self % pseudo_ops, &
153  self % vert_interp_ops, &
154  1, &
155  impact_param(iobs:iobs), &
156  radius_curv(iobs), &
157  obslat(iobs), &
158  undulation(iobs), &
159  hofx(iobs:iobs), &
160  baerr)
161  else
162  call ops_gpsro_forwardmodel(prs % nval, &
163  q % nval, &
164  rho_heights % vals(:,iobs), &
165  theta_heights % vals(:,iobs), &
166  prs % vals(:,iobs), &
167  q % vals(:,iobs), &
168  self % pseudo_ops, &
169  self % vert_interp_ops, &
170  1, &
171  impact_param(iobs:iobs), &
172  radius_curv(iobs), &
173  obslat(iobs), &
174  undulation(iobs), &
175  hofx(iobs:iobs), &
176  baerr)
177  end if
178 
179  if (baerr) then
180  write(err_msg,*) "Error with observation processing ", iobs
181  call fckit_log % info(err_msg)
182  end if
183 
184  end do obs_loop
185 
186  deallocate(obslat)
187  deallocate(obslon)
188  deallocate(impact_param)
189  deallocate(radius_curv)
190  deallocate(undulation)
191 
192  write(err_msg,*) "TRACE: ufo_gnssro_bendmetoffice_simobs: completed"
193  call fckit_log%info(err_msg)
194 
195 end subroutine ufo_gnssro_bendmetoffice_simobs
196 ! ------------------------------------------------------------------------------
197 
198 
199 SUBROUTINE ops_gpsro_forwardmodel(nlevP, &
200  nlevq, &
201  za, &
202  zb, &
203  pressure, &
204  humidity, &
205  GPSRO_pseudo_ops, &
206  GPSRO_vert_interp_ops, &
207  nobs, &
208  zobs, &
209  RO_Rad_Curv, &
210  Latitude, &
211  RO_geoid_und, &
212  ycalc, &
213  BAErr)
214 
215 INTEGER, INTENT(IN) :: nlevP ! no. of p levels in state vec.
216 INTEGER, INTENT(IN) :: nlevq ! no. of theta levels
217 REAL(kind_real), INTENT(IN) :: za(1:nlevP) ! heights of rho levs
218 REAL(kind_real), INTENT(IN) :: zb(1:nlevQ) ! heights of theta levs
219 REAL(kind_real), INTENT(IN) :: pressure(1:nlevP) ! Model background pressure
220 REAL(kind_real), INTENT(IN) :: humidity(1:nlevQ) ! Model background specific humidity
221 LOGICAL, INTENT(IN) :: GPSRO_pseudo_ops ! Option: Use pseudo-levels in vertical interpolation?
222 LOGICAL, INTENT(IN) :: GPSRO_vert_interp_ops ! Option: Use ln(p) for vertical interpolation? (rather than exner)
223 INTEGER, INTENT(IN) :: nobs ! Number of observations in the profile
224 REAL(kind_real), INTENT(IN) :: zobs(1:nobs) ! Impact parameter for the obs
225 REAL(kind_real), INTENT(IN) :: RO_Rad_Curv ! Earth's radius of curvature for these observations
226 REAL(kind_real), INTENT(IN) :: Latitude ! Latitude of this profile
227 REAL(kind_real), INTENT(IN) :: RO_geoid_und ! Undulation - difference between the geoid and the ellipsoid
228 REAL(kind_real), INTENT(INOUT) :: ycalc(1:nobs) ! Model forecast of the observations
229 LOGICAL, INTENT(OUT) :: BAErr ! Was an error encountered during the calculation?
230 !
231 ! Things that may need to be output, as they are used by the TL/AD calculation
232 !
233 REAL(kind_real), ALLOCATABLE :: z_pseudo(:) ! Heights of the pseudo levels | Allocated by
234 REAL(kind_real), ALLOCATABLE :: N_pseudo(:) ! Refractivity on the pseudo levels | Ops_GPSRO_refrac
235 INTEGER :: nb_pseudo ! Number of pseudo levels
236 REAL(kind_real) :: T(1:nlevq) ! Temperature on model levels
237 REAL(kind_real), ALLOCATABLE :: nr(:) ! Model calculation of impact parameters
238 REAL(kind_real) :: Refmodel(1:nlevq) ! model refractivity on theta levels
239 !
240 ! Local parameters
241 !
242 integer, parameter :: max_string = 800 ! Length of strings
243 character(len=*), parameter :: myname_ = "Ops_GPSRO_ForwardModel"
244 !
245 ! Local variables
246 !
247 INTEGER :: nstate ! no. of levels in state vec.
248 INTEGER :: num_pseudo ! Number of levels, including pseudo levels
249 INTEGER :: nb ! no. of non-pseudo levs
250 REAL(kind_real) :: x(1:nlevP+nlevQ) ! state vector
251 character(max_string) :: err_msg ! Error message to be output
252 
253 ! The model data must be on a staggered grid, with nlevp = nlevq+1
254 IF (nlevp /= nlevq + 1) THEN
255  write(err_msg,*) myname_ // ':' // ' Data must be on a staggered grid nlevp, nlevq = ', nlevp, nlevq
256  call fckit_log % warning(err_msg)
257  write(err_msg,*) myname_ // ':' // ' error: number of levels inconsistent!'
258  call abor1_ftn(err_msg)
259 END IF
260 
261 nstate = nlevp + nlevq
262 nb = nlevq
263 x(1:nlevp) = pressure
264 x(nlevp+1:nstate) = humidity
265 
266 IF (gpsro_pseudo_ops) THEN
267  num_pseudo = 2 * nlevq - 1
268 ELSE
269  num_pseudo = nlevq
270 END IF
271 ALLOCATE(nr(1:num_pseudo))
272 
273 baerr = .false.
274 
275 CALL ops_gpsro_refrac (nstate, &
276  nlevp, &
277  nb, &
278  nlevq, &
279  za, &
280  zb, &
281  x, &
282  gpsro_pseudo_ops, &
283  gpsro_vert_interp_ops, &
284  baerr, &
285  refmodel, &
286  t, &
287  z_pseudo, &
288  n_pseudo, &
289  nb_pseudo)
290 
291 ! no point proceeding further if ...
292 IF (.NOT. baerr) THEN
293  ! Pseudo levels
294  IF (gpsro_pseudo_ops) THEN
295  ! 2. Calculate the refractive index * radius on theta model levels (or model impact parameter)
296  CALL ops_gpsrocalc_nr (z_pseudo, & ! geopotential heights of pseudo levels
297  nb_pseudo, & ! number of model+pseudo-levels
298  ro_rad_curv, & ! radius of curvature of earth at observation
299  latitude, & ! latitude at observation
300  ro_geoid_und, & ! geoid undulation above WGS-84
301  n_pseudo, & ! refractivity of model on model+pseudo levels
302  nr) ! Calculated model impact parameters
303 
304  ! 3. Calculate model bending angle on observation impact parameters
305  CALL ops_gpsrocalc_alpha (nobs, & ! size of ob. vector
306  nb_pseudo, & ! no. of refractivity levels
307  zobs, & ! obs impact parameters
308  n_pseudo, & ! refractivity values on model+pseudo levels
309  nr, & ! index * radius product
310  ycalc) ! forward modelled bending angle
311  ! Model levels only
312  ELSE
313  ! 2. Calculate the refractive index * radius on theta model levels (or model impact parameter)
314  CALL ops_gpsrocalc_nr (zb, & ! geopotential heights of model levels
315  nb, & ! number of levels in zb
316  ro_rad_curv, & ! radius of curvature of earth at observation
317  latitude, & ! latitude at observation
318  ro_geoid_und, & ! geoid undulation above WGS-84
319  refmodel, & ! refractivity of model on model levels
320  nr) ! Calculated model impact parameters
321 
322  ! 3. Calculate model bending angle on observation impact parameters
323  CALL ops_gpsrocalc_alpha (nobs, & ! size of ob. vector
324  nb, & ! no. of refractivity levels
325  zobs, & ! obs impact parameters
326  refmodel, & ! refractivity values on model levels
327  nr, & ! index * radius product
328  ycalc) ! forward modelled bending angle
329  END IF
330 END IF
331 
332 END SUBROUTINE ops_gpsro_forwardmodel
333 
334 end module ufo_gnssro_bendmetoffice_mod
ufo_gnssro_bendmetoffice_mod
Fortran module for gnssro bending angle Met Office forward operator.
Definition: ufo_gnssro_bendmetoffice_mod.F90:10
ufo_gnssro_bendmetoffice_mod::ops_gpsro_forwardmodel
subroutine ops_gpsro_forwardmodel(nlevP, nlevq, za, zb, pressure, humidity, GPSRO_pseudo_ops, GPSRO_vert_interp_ops, nobs, zobs, RO_Rad_Curv, Latitude, RO_geoid_und, ycalc, BAErr)
Definition: ufo_gnssro_bendmetoffice_mod.F90:214
lag_interp_mod
Fortran module to prepare for Lagrange polynomial interpolation. based on GSI: lagmod....
Definition: lag_interp.F90:4
ufo_vars_mod::var_zi
character(len=maxvarlen), parameter, public var_zi
Definition: ufo_variables_mod.F90:30
ufo_basis_mod
Definition: ufo_basis_mod.F90:6
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
lag_interp_mod::lag_interp_smthweights
subroutine, public lag_interp_smthweights(x, xt, aq, bq, w, dw, n)
Definition: lag_interp.F90:174
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_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_basis_mod::ufo_basis
Definition: ufo_basis_mod.F90:12
lag_interp_mod::lag_interp_const
subroutine, public lag_interp_const(q, x, n)
Definition: lag_interp.F90:24
ufo_vars_mod
Definition: ufo_variables_mod.F90:8
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:51
ufo_gnssro_bendmetoffice_mod::ufo_gnssro_bendmetoffice_simobs
subroutine ufo_gnssro_bendmetoffice_simobs(self, geovals, hofx, obss)
Definition: ufo_gnssro_bendmetoffice_mod.F90:60
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_geovals_mod::ufo_geovals
type to hold interpolated fields required by the obs operators
Definition: ufo_geovals_mod.F90:47
ufo_geovals_mod::ufo_geoval
type to hold interpolated field for one variable, one observation
Definition: ufo_geovals_mod.F90:40
ufo_gnssro_bendmetoffice_mod::ufo_gnssro_bendmetoffice
Fortran derived type for gnssro trajectory.
Definition: ufo_gnssro_bendmetoffice_mod.F90:30
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