UFO
ufo_gnssro_bndropp1d_tlad_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2017-2018 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 !> Fortran module for gnssro bending angle ropp1d tangent linear and adjoint
7 !> following the ROPP (2018 Aug) implementation
8 
9 module ufo_gnssro_bndropp1d_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 ufo_gnssro_ropp1d_utils_mod
22 use fckit_log_module, only : fckit_log
23 
24 integer, parameter :: max_string=800
25 
26 !> Fortran derived type for gnssro trajectory
27 type, extends(ufo_basis_tlad) :: ufo_gnssro_bndropp1d_tlad
28  private
29  integer :: nval, nlocs, iflip
30  real(kind_real), allocatable :: prs(:,:), t(:,:), q(:,:), gph(:,:), gph_sfc(:,:)
31  contains
32  procedure :: delete => ufo_gnssro_bndropp1d_tlad_delete
33  procedure :: settraj => ufo_gnssro_bndropp1d_tlad_settraj
34  procedure :: simobs_tl => ufo_gnssro_bndropp1d_simobs_tl
35  procedure :: simobs_ad => ufo_gnssro_bndropp1d_simobs_ad
36 end type ufo_gnssro_bndropp1d_tlad
37 
38 contains
39 
40 ! ------------------------------------------------------------------------------
41 ! ------------------------------------------------------------------------------
42 subroutine ufo_gnssro_bndropp1d_tlad_settraj(self, geovals, obss)
43 
44  implicit none
45  class(ufo_gnssro_bndropp1d_tlad), intent(inout) :: self
46  type(ufo_geovals), intent(in) :: geovals
47  type(c_ptr), value, intent(in) :: obss
48  character(len=*), parameter :: myname_="ufo_gnssro_bndropp1d_tlad_settraj"
49  character(max_string) :: err_msg
50  type(ufo_geoval), pointer :: t, q, prs, gph, gph_sfc
51 
52  write(err_msg,*) "TRACE: ufo_gnssro_bndropp1d_tlad_settraj: begin"
53  call fckit_log%info(err_msg)
54 
55 ! get model state variables from geovals
56  call ufo_geovals_get_var(geovals, var_ts, t) ! temperature
57  call ufo_geovals_get_var(geovals, var_q, q) ! specific humidity
58  call ufo_geovals_get_var(geovals, var_prs, prs) ! pressure
59  call ufo_geovals_get_var(geovals, var_z, gph) ! geopotential height
60  call ufo_geovals_get_var(geovals, var_sfc_geomz, gph_sfc) ! surface geopotential height
61 
62  call self%delete()
63 
64 ! Keep copy of dimensions
65  self%nval = prs%nval
66  self%nlocs = obsspace_get_nlocs(obss)
67 
68  if (self%nlocs > 0) then
69  self%iflip = 0
70  if (prs%vals(1,1) .lt. prs%vals(prs%nval,1) ) then
71  self%iflip = 1
72  write(err_msg,'(a)') ' ufo_gnssro_bndropp1d_tlad_settraj:'//new_line('a')// &
73  ' Model vertical height profile is in descending order,'//new_line('a')// &
74  ' but ROPP requires it to be ascending order, need flip'
75  call fckit_log%info(err_msg)
76  end if
77 
78  allocate(self%t(self%nval,self%nlocs))
79  allocate(self%q(self%nval,self%nlocs))
80  allocate(self%prs(self%nval,self%nlocs))
81  allocate(self%gph(self%nval,self%nlocs))
82  allocate(self%gph_sfc(1,self%nlocs))
83  self%gph = gph%vals
84  self%t = t%vals
85  self%q = q%vals
86  self%prs = prs%vals
87  self%gph_sfc = gph_sfc%vals
88  end if
89  self%ltraj = .true.
90 
91 end subroutine ufo_gnssro_bndropp1d_tlad_settraj
92 
93 ! ------------------------------------------------------------------------------
94 ! ------------------------------------------------------------------------------
95 subroutine ufo_gnssro_bndropp1d_simobs_tl(self, geovals, hofx, obss)
96 
97  use ropp_fm_types, only: state1dfm
98  use ropp_fm_types, only: obs1dbangle
99  use datetimetypes, only: dp
100  implicit none
101  class(ufo_gnssro_bndropp1d_tlad), intent(in) :: self
102  type(ufo_geovals), intent(in) :: geovals ! perturbed quantities
103  real(kind_real), intent(inout) :: hofx(:)
104  type(c_ptr), value, intent(in) :: obss
105 
106  type(state1dfm) :: x,x_tl
107  type(obs1dbangle) :: y,y_tl
108 
109  integer :: iobs,nlev, nlocs
110  integer :: nvprof
111 
112  character(len=*), parameter :: myname_="ufo_gnssro_bndropp1d_simobs_tl"
113  character(max_string) :: err_msg
114  real(kind=dp) :: ob_time
115  type(ufo_geoval), pointer :: t_d, q_d, prs_d
116 
117 ! hack - set local geopotential height to zero for ropp routines
118  real(kind_real), allocatable :: gph_d_zero(:)
119  real(kind_real) :: gph_sfc_d_zero
120  real(kind_real), allocatable :: obsLat(:), obsLon(:), obsImpP(:), obsLocR(:), obsGeoid(:)
121 ! hack - set local geopotential height to zero for ropp routines
122 
123  write(err_msg,*) "TRACE: ufo_gnssro_bndropp1d_simobs_tl: begin"
124  call fckit_log%info(err_msg)
125 
126 ! check if trajectory was set
127  if (.not. self%ltraj) then
128  write(err_msg,*) myname_, ' trajectory wasnt set!'
129  call abor1_ftn(err_msg)
130  endif
131 
132 ! check if nlocs is consistent in geovals & hofx
133  if (geovals%nlocs /= size(hofx)) then
134  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
135  call abor1_ftn(err_msg)
136  endif
137 
138 ! get variables from geovals
139  call ufo_geovals_get_var(geovals, var_ts, t_d) ! temperature
140  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
141  call ufo_geovals_get_var(geovals, var_prs, prs_d) ! pressure
142 
143  nlev = self%nval
144  nlocs = self%nlocs ! number of observations
145  if (nlocs > 0) then
146  allocate(gph_d_zero(nlev))
147  gph_d_zero = 0.0
148  gph_sfc_d_zero = 0.0
149 
150  ! set obs space struture
151  allocate(obslon(nlocs))
152  allocate(obslat(nlocs))
153  allocate(obsimpp(nlocs))
154  allocate(obslocr(nlocs))
155  allocate(obsgeoid(nlocs))
156  call obsspace_get_db(obss, "MetaData", "longitude", obslon)
157  call obsspace_get_db(obss, "MetaData", "latitude", obslat)
158  call obsspace_get_db(obss, "MetaData", "impact_parameter", obsimpp)
159  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", obslocr)
160  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", obsgeoid)
161 
162  nvprof = 1 ! no. of bending angles in profile
163 
164  ! loop through the obs
165  obs_loop: do iobs = 1, nlocs ! order of loop doesn't matter
166 
167  ob_time = 0.0
168  ! map the trajectory to ROPP structure x
169  call init_ropp_1d_statevec( ob_time, &
170  obslon(iobs), &
171  obslat(iobs), &
172  self%t(:,iobs), &
173  self%q(:,iobs), &
174  self%prs(:,iobs), &
175  self%gph(:,iobs), &
176  nlev, &
177  self%gph_sfc(1,iobs), &
178  x, self%iflip)
179  ! hack -- make non zero humidity to avoid zero denominator in tangent linear
180  ! see ropp_fm/bangle_1d/ropp_fm_bangle_1d_tl.f90
181  where(x%shum .le. 1e-8) x%shum = 1e-8
182  ! hack -- make non zero humidity to avoid zero denominator in tangent linear
183 
184  call init_ropp_1d_statevec( ob_time, &
185  obslon(iobs), &
186  obslat(iobs), &
187  t_d%vals(:,iobs), &
188  q_d%vals(:,iobs), &
189  prs_d%vals(:,iobs), &
190  gph_d_zero(:), &
191  nlev, &
192  gph_sfc_d_zero, &
193  x_tl, self%iflip)
194  ! set both y and y_tl structures
195  call init_ropp_1d_obvec_tlad(iobs, nvprof, &
196  obsimpp(iobs), &
197  obslat(iobs), &
198  obslon(iobs), &
199  obslocr(iobs), &
200  obsgeoid(iobs), &
201  y,y_tl)
202 
203  ! now call TL of forward model
204  call ropp_fm_bangle_1d_tl(x,x_tl,y, y_tl%bangle(nvprof))
205  hofx(iobs) = y_tl%bangle(nvprof) ! this will need to change if profile is passed
206 
207  ! tidy up -deallocate ropp structures
208  call ropp_tidy_up_tlad_1d(x,x_tl,y,y_tl)
209  end do obs_loop
210 
211  ! tidy up - deallocate obsspace structures
212  deallocate(obslat)
213  deallocate(obslon)
214  deallocate(obsimpp)
215  deallocate(obslocr)
216  deallocate(obsgeoid)
217  end if ! nlocs > 0
218 
219  write(err_msg,*) "TRACE: ufo_gnssro_bndropp1d_simobs_tl: complete"
220  call fckit_log%info(err_msg)
221 
222  return
223 
224 end subroutine ufo_gnssro_bndropp1d_simobs_tl
225 
226 ! ------------------------------------------------------------------------------
227 ! ------------------------------------------------------------------------------
228 subroutine ufo_gnssro_bndropp1d_simobs_ad(self, geovals, hofx, obss)
229 
230  use ropp_fm_types, only: state1dfm
231  use ropp_fm_types, only: obs1dbangle
232  use typesizes, only: wp => eightbytereal
233  use datetimetypes, only: dp
234 
235  implicit none
236  class(ufo_gnssro_bndropp1d_tlad), intent(in) :: self
237  type(ufo_geovals), intent(inout) :: geovals ! perturbed quantities
238  real(kind_real), intent(in) :: hofx(:)
239  type(c_ptr), value, intent(in) :: obss
240  real(c_double) :: missing
241 
242  type(ufo_geoval), pointer :: t_d, q_d, prs_d
243 ! set local geopotential height to zero for ropp routines
244  real(kind_real), parameter :: gph_sfc_d_zero = 0.0
245  real(kind_real), allocatable :: gph_d_zero(:)
246 
247  real(kind_real), allocatable :: obsLat(:), obsLon(:), obsImpP(:), obsLocR(:), obsGeoid(:)
248  type(state1dfm) :: x,x_ad
249  type(obs1dbangle) :: y,y_ad
250  integer :: iobs,nlev, nlocs
251  integer :: nvprof
252  real(kind=dp) :: ob_time
253  character(len=*), parameter :: myname_="ufo_gnssro_bndropp1d_simobs_ad"
254  character(max_string) :: err_msg
255 
256  write(err_msg,*) "TRACE: ufo_gnssro_bndropp1d_simobs_ad: begin"
257  call fckit_log%info(err_msg)
258  if (self%nlocs > 0) then
259  ! check if trajectory was set
260  if (.not. self%ltraj) then
261  write(err_msg,*) myname_, ' trajectory wasnt set!'
262  call abor1_ftn(err_msg)
263  endif
264  ! check if nlocs is consistent in geovals & hofx
265  if (geovals%nlocs /= size(hofx)) then
266  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
267  call abor1_ftn(err_msg)
268  endif
269 
270  ! get variables from geovals
271  call ufo_geovals_get_var(geovals, var_ts, t_d) ! temperature
272  call ufo_geovals_get_var(geovals, var_q, q_d) ! specific humidity
273  call ufo_geovals_get_var(geovals, var_prs, prs_d) ! pressure
274 
275  ! allocate if not yet allocated
276  if (.not. allocated(t_d%vals)) then
277  t_d%nlocs = self%nlocs
278  t_d%nval = self%nval
279  allocate(t_d%vals(t_d%nval,t_d%nlocs))
280  t_d%vals = 0.0_kind_real
281  endif
282 
283  if (.not. allocated(prs_d%vals)) then
284  prs_d%nlocs = self%nlocs
285  prs_d%nval = self%nval
286  allocate(prs_d%vals(prs_d%nval,prs_d%nlocs))
287  prs_d%vals = 0.0_kind_real
288  endif
289 
290  if (.not. allocated(q_d%vals)) then
291  q_d%nlocs = self%nlocs
292  q_d%nval = self%nval
293  allocate(q_d%vals(q_d%nval,q_d%nlocs))
294  q_d%vals = 0.0_kind_real
295  endif
296 
297  if (.not. geovals%linit ) geovals%linit=.true.
298 
299  nlev = self%nval
300  nlocs = self%nlocs
301 
302  allocate(gph_d_zero(nlev))
303  gph_d_zero = 0.0
304 
305  ! set obs space struture
306  allocate(obslon(nlocs))
307  allocate(obslat(nlocs))
308  allocate(obsimpp(nlocs))
309  allocate(obslocr(nlocs))
310  allocate(obsgeoid(nlocs))
311 
312  call obsspace_get_db(obss, "MetaData", "longitude", obslon)
313  call obsspace_get_db(obss, "MetaData", "latitude", obslat)
314  call obsspace_get_db(obss, "MetaData", "impact_parameter", obsimpp)
315  call obsspace_get_db(obss, "MetaData", "earth_radius_of_curvature", obslocr)
316  call obsspace_get_db(obss, "MetaData", "geoid_height_above_reference_ellipsoid", obsgeoid)
317 
318  missing = missing_value(missing)
319 
320  ! loop through the obs
321  nvprof=1 ! no. of bending angles in profile
322  obs_loop: do iobs = 1, nlocs
323 
324  if (hofx(iobs) .gt. missing) then
325  ob_time = 0.0
326 
327  ! map the trajectory to ROPP structure x
328  call init_ropp_1d_statevec( ob_time, &
329  obslon(iobs), &
330  obslat(iobs), &
331  self%t(:,iobs), &
332  self%q(:,iobs), &
333  self%prs(:,iobs), &
334  self%gph(:,iobs), &
335  nlev, &
336  self%gph_sfc(1,iobs),&
337  x, self%iflip)
338 
339  call init_ropp_1d_statevec( ob_time, &
340  obslon(iobs), &
341  obslat(iobs), &
342  t_d%vals(:,iobs), &
343  q_d%vals(:,iobs), &
344  prs_d%vals(:,iobs), &
345  gph_d_zero(:), &
346  nlev, &
347  gph_sfc_d_zero, &
348  x_ad, self%iflip)
349 
350  ! x_ad is local so initialise to 0.0
351  x_ad%temp(:) = 0.0_wp
352  x_ad%pres(:) = 0.0_wp
353  x_ad%shum(:) = 0.0_wp
354  x_ad%geop(:) = 0.0_wp
355 
356  ! set both y and y_ad structures
357  call init_ropp_1d_obvec_tlad(iobs, nvprof, &
358  obsimpp(iobs), &
359  obslat(iobs), &
360  obslon(iobs), &
361  obslocr(iobs), &
362  obsgeoid(iobs), &
363  y,y_ad)
364 
365  ! local variable initialise
366  y_ad%bangle(:) = 0.0_wp
367 
368  ! now call AD of forward model
369  y_ad%bangle(nvprof) = y_ad%bangle(nvprof) + hofx(iobs)
370  call ropp_fm_bangle_1d_ad(x,x_ad,y,y_ad)
371  call init_ropp_1d_statevec_ad( &
372  t_d%vals(:,iobs), &
373  q_d%vals(:,iobs), &
374  prs_d%vals(:,iobs), &
375  gph_d_zero(:), &
376  nlev, x_ad, self%iflip)
377 
378  ! tidy up - deallocate ropp structures
379  call ropp_tidy_up_tlad_1d(x,x_ad,y,y_ad)
380 
381  end if ! end missing value check
382 
383  end do obs_loop
384 
385  ! tidy up - deallocate obsspace structures
386  deallocate(obslat)
387  deallocate(obslon)
388  deallocate(obsimpp)
389  deallocate(obslocr)
390  deallocate(obsgeoid)
391  deallocate(gph_d_zero)
392  end if ! nlocs > 0
393 
394  write(err_msg,*) "TRACE: ufo_gnssro_bndropp1d_simobs_ad: complete"
395  call fckit_log%info(err_msg)
396 
397  return
398 
399 end subroutine ufo_gnssro_bndropp1d_simobs_ad
400 
401 !-------------------------------------------------------------------------
402 !-------------------------------------------------------------------------
403 subroutine ufo_gnssro_bndropp1d_tlad_delete(self)
404 
405  implicit none
406  class(ufo_gnssro_bndropp1d_tlad), intent(inout) :: self
407  character(len=*), parameter :: myname_="ufo_gnssro_bndropp_tlad_delete"
408 
409  self%nval = 0
410  if (allocated(self%prs)) deallocate(self%prs)
411  if (allocated(self%t)) deallocate(self%t)
412  if (allocated(self%q)) deallocate(self%q)
413  if (allocated(self%gph)) deallocate(self%gph)
414  if (allocated(self%gph_sfc)) deallocate(self%gph_sfc)
415  self%ltraj = .false.
416 
417 end subroutine ufo_gnssro_bndropp1d_tlad_delete
418 
419 !-------------------------------------------------------------------------
420 
421 end module ufo_gnssro_bndropp1d_tlad_mod
ufo_avgkernel_mod::max_string
integer, parameter max_string
Definition: ufo_avgkernel_mod.F90:17
ufo_basis_tlad_mod
Definition: ufo_basis_tlad_mod.F90:6
ufo_gnssro_ropp1d_utils_mod::init_ropp_1d_statevec
subroutine, public init_ropp_1d_statevec(step_time, rlon, rlat, temp, shum, pres, phi, lm, phi_sfc, x, iflip)
Definition: ufo_gnssro_ropp1d_utils_mod.F90:37
ufo_gnssro_ropp1d_utils_mod::init_ropp_1d_obvec_tlad
subroutine, public init_ropp_1d_obvec_tlad(iloop, nvprof, obs_impact, rlat, rlon, roc, undulat, y, y_p)
Definition: ufo_gnssro_ropp1d_utils_mod.F90:326
ufo_gnssro_bndropp1d_tlad_mod::ufo_gnssro_bndropp1d_tlad_delete
subroutine ufo_gnssro_bndropp1d_tlad_delete(self)
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:404
ufo_gnssro_ropp1d_utils_mod::ropp_tidy_up_tlad_1d
subroutine, public ropp_tidy_up_tlad_1d(x, x_p, y, y_p)
Definition: ufo_gnssro_ropp1d_utils_mod.F90:399
ufo_gnssro_ropp1d_utils_mod::init_ropp_1d_statevec_ad
subroutine, public init_ropp_1d_statevec_ad(temp_d, shum_d, pres_d, phi_d, lm, x_ad, iflip)
Definition: ufo_gnssro_ropp1d_utils_mod.F90:166
ufo_gnssro_bndropp1d_tlad_mod::ufo_gnssro_bndropp1d_tlad_settraj
subroutine ufo_gnssro_bndropp1d_tlad_settraj(self, geovals, obss)
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:43
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_gnssro_bndropp1d_tlad_mod
Fortran module for gnssro bending angle ropp1d tangent linear and adjoint following the ROPP (2018 Au...
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:9
ufo_vars_mod::var_sfc_geomz
character(len=maxvarlen), parameter, public var_sfc_geomz
Definition: ufo_variables_mod.F90:70
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_bndropp1d_tlad_mod::ufo_gnssro_bndropp1d_tlad
Fortran derived type for gnssro trajectory.
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:27
ufo_vars_mod::var_ts
character(len=maxvarlen), parameter, public var_ts
Definition: ufo_variables_mod.F90:19
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_gnssro_ropp1d_utils_mod
Fortran module to handle gnssro bending angle observations following the ROPP (2018 Aug) implementati...
Definition: ufo_gnssro_ropp1d_utils_mod.F90:9
ufo_geovals_mod::ufo_geoval
type to hold interpolated field for one variable, one observation
Definition: ufo_geovals_mod.F90:40
ufo_gnssro_bndropp1d_tlad_mod::ufo_gnssro_bndropp1d_simobs_tl
subroutine ufo_gnssro_bndropp1d_simobs_tl(self, geovals, hofx, obss)
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:96
ufo_vars_mod::var_prs
character(len=maxvarlen), parameter, public var_prs
Definition: ufo_variables_mod.F90:25
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
ufo_gnssro_bndropp1d_tlad_mod::ufo_gnssro_bndropp1d_simobs_ad
subroutine ufo_gnssro_bndropp1d_simobs_ad(self, geovals, hofx, obss)
Definition: ufo_gnssro_bndropp1d_tlad_mod.F90:229