UFO
ufo_gnssro_ropp2d_utils_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 to handle gnssro bending angle observations following
7 !> the ROPP (2018 Aug) implementation
8 
9 module ufo_gnssro_ropp2d_utils_mod
10 
11 !use iso_c_binding
12 use fckit_log_module, only: fckit_log
13 use kinds, only: kind_real
14 
15 ! ROPP data type and library subroutines
16 use typesizes, only: wp => eightbytereal
17 use datetimetypes, only: dp
18 use ropp_fm_types, only: state2dfm, obs1dbangle
19 use geodesy, only: gravity, r_eff, geometric2geopotential
20 use arrays, only: callocate
21 
22 implicit none
23 public :: init_ropp_2d_statevec
24 public :: init_ropp_2d_statevec_ad
25 public :: init_ropp_2d_obvec
26 public :: init_ropp_2d_obvec_tlad
27 public :: ropp_tidy_up_2d
28 public :: ropp_tidy_up_tlad_2d
29 private
30 
31 contains
32 
33 ! ------------------------------------------------------------------------------------
34 ! ------------------------------------------------------------------------------------
35 subroutine init_ropp_2d_statevec(rlon,rlat,temp,shum,pres,phi,lm, x, n_horiz, dtheta, iflip)
36 
37 ! Description:
38 ! subroutine to fill a ROPP state vector structure with
39 ! model background fields: Temperature, pressure, specific
40 ! humidity at full-levels, and surface geopotential height.
41 !
42 ! Inputs:
43 ! temp background temperature
44 ! shum background specific humidity
45 ! pres background pressure
46 ! phi geopotential height
47 !
48 ! phi_sfc terrain geopotential of background field
49 ! lm number of vertical levels in the background
50 !
51 ! Outputs:
52 ! x Forward model state vector
53 !
54 ! ###############################################################
55  implicit none
56 ! Output state vector
57  type(state2dfm), intent(out) :: x
58  integer, intent(in) :: lm, n_horiz
59  real(kind=kind_real), intent(in) :: dtheta
60  real(kind=kind_real), dimension(n_horiz), intent(in) :: rlon, rlat
61  real(kind=kind_real), dimension(lm,n_horiz), intent(in) :: temp,shum,pres,phi
62 ! Local variables
63  integer :: n,k
64  integer, optional, intent(in) :: iflip
65 !-------------------------------------------------------------------------
66 ! number of profiles in plane
67  x%n_horiz = n_horiz
68  x%dtheta = dtheta
69 
70 ! Number of levels in background profile. What about (lm+1) field ?
71  x%n_lev=lm
72 
73 !-------------------------------------------------------------
74 ! allocate arrays for temperature, specific humidity, pressure
75 ! and geopotential height data and additional 2d fields
76 !-------------------------------------------------------------
77  allocate(x%temp(x%n_lev,x%n_horiz))
78  allocate(x%shum(x%n_lev,x%n_horiz))
79  allocate(x%pres(x%n_lev,x%n_horiz))
80  allocate(x%geop(x%n_lev,x%n_horiz))
81 
82 ! needed in 2D
83  allocate(x%refrac(x%n_lev,x%n_horiz))
84  allocate(x%nr(x%n_lev,x%n_horiz))
85  allocate(x%lat(x%n_horiz))
86  allocate(x%lon(x%n_horiz))
87 
88  x%lat(:) = real(rlat(:),kind=wp)
89  x%lon(:) = real(rlon(:),kind=wp)
90  where (x%lon .gt. 180.0) x%lon = x%lon -360.0
91 
92 
93 ! TEMPORARY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
94 !----------------------------------------------------
95 ! ROPP FM requires vertical height profile to be of the ascending order.
96 ! (see ropp_io_ascend ( ROdata )). So we need to flip the data.
97 !----------------------------------------------------
98 
99  n = lm
100  if ( present(iflip) .and. iflip .eq. 1) then
101 
102  do k = 1, lm
103  x%temp(n,:) = real(temp(k,:),kind=wp)
104  x%shum(n,:) = real(shum(k,:),kind=wp)
105  x%pres(n,:) = real(pres(k,:),kind=wp)
106  x%geop(n,:) = real(phi(k,:),kind=wp)
107  n = n - 1
108  end do
109 
110  else
111  do k = 1, lm
112  x%temp(k,:) = real(temp(k,:),kind=wp)
113  x%shum(k,:) = real(shum(k,:),kind=wp)
114  x%pres(k,:) = real(pres(k,:),kind=wp)
115  x%geop(k,:) = real(phi(k,:),kind=wp)
116  end do
117 
118  end if
119  return
120 end subroutine init_ropp_2d_statevec
121 
122 ! ------------------------------------------------------------------------------
123 
124 subroutine init_ropp_2d_statevec_ad(temp_d,shum_d,pres_d,phi_d,lm,x_ad,n_horiz,iflip)
125 
126 ! Description:
127 ! subroutine to fill a ROPP state vector structure with
128 ! model background fields: Temperature, pressure, specific
129 ! humidity at full-levels, and surface geopotential height.
130 !
131 ! Inputs:
132 ! temp background temperature
133 ! shum background specific humidity
134 ! pres background pressure
135 ! phi geopotential height
136 !
137 ! phi_sfc terrain geopotential of background field
138 ! lm number of vertical levels in the background
139 !
140 ! Outputs:
141 ! x Forward model state vector
142 !
143 ! ###############################################################
144  implicit none
145 
146 ! Function arguments
147 
148 ! Output state vector
149 
150  type(state2dfm), intent(inout) :: x_ad
151  integer, intent(in) :: lm, n_horiz
152  real(kind=kind_real), &
153  dimension(lm,n_horiz), intent(inout) :: temp_d,shum_d,pres_d,phi_d
154 
155 ! Local variables
156  integer :: n,j,k
157  integer, optional, intent(in) :: iflip
158 !-------------------------------------------------------------------------
159  n = lm
160  x_ad%n_horiz = n_horiz
161 
162  if ( present(iflip) .and. iflip .eq. 1) then
163 
164  do k = 1, lm
165  do j = 1, x_ad%n_horiz
166 !!! x_tl%temp(n,:) = real(temp_d(k),kind=wp)
167  temp_d(k,j) = temp_d(k,j) + real(x_ad%temp(n,j),kind=kind_real)
168  x_ad%temp(n,j) = 0.0_wp
169 
170 !!! x_tl%shum(n,:) = real(shum_d(k),kind=wp)
171  shum_d(k,j) = shum_d(k,j) + real(x_ad%shum(n,j),kind=kind_real)
172  x_ad%shum(n,j) = 0.0_wp
173 
174 !!! x_tl%pres(n,:) = real(pres_d(k),kind=wp)
175  pres_d(k,j) = pres_d(k,j) + real(x_ad%pres(n,j),kind=kind_real)
176  x_ad%pres(n,j) = 0.0_wp
177 
178 !!! x_tl%geop(n,:) = real(phi_d(k),kind=wp)
179  phi_d(k,j) = phi_d(k,j) + real(x_ad%geop(n,j),kind=kind_real)
180  x_ad%geop(n,j) = 0.0_wp
181  enddo
182 
183  n = n - 1
184  end do
185  else
186  do k = 1, lm
187  do j = 1, x_ad%n_horiz
188  temp_d(k,j) = temp_d(k,j) + real(x_ad%temp(k,j),kind=kind_real)
189  x_ad%temp(k,j) = 0.0_wp
190  shum_d(k,j) = shum_d(k,j) + real(x_ad%shum(k,j),kind=kind_real)
191  x_ad%shum(k,j) = 0.0_wp
192  pres_d(k,j) = pres_d(k,j) + real(x_ad%pres(k,j),kind=kind_real)
193  x_ad%pres(k,j) = 0.0_wp
194  phi_d(k,j) = phi_d(k,j) + real(x_ad%geop(k,j),kind=kind_real)
195  x_ad%geop(k,j) = 0.0_wp
196  enddo
197  end do
198  end if
199  return
200  end subroutine init_ropp_2d_statevec_ad
201 !-------------------------------------------------------------------------
202 !---------------------------------------------------------------------------------------------
203 subroutine init_ropp_2d_obvec(nvprof,obs_impact,rlat,rlon,roc,undulat,y)
204 
205 ! Description:
206 ! subroutine to fill a ROPP observation vector structure
207 ! observation provides the inputs of:
208 ! impact parameter, lat, lon, time, radius of curvature
209 !
210 ! forward model will provide the simulation of bending angle
211 !
212 ! Inputs:
213 ! obs_impact impact parameter
214 ! rlat latitude
215 ! rlon longitude
216 ! roc radius of curvature
217 ! undulat undulation correction for radius of curvature
218 !
219 ! Outputs:
220 ! y: Partially filled Forward model observation vector
221 !
222 ! ###############################################################
223  implicit none
224 ! Output state vector
225  type(obs1dbangle), intent(out) :: y
226 
227  integer, intent(in) :: nvprof
228  real(kind=kind_real), dimension(nvprof), intent(in) :: obs_impact
229  real(kind=kind_real), intent(in) :: rlat,rlon
230  real(kind=kind_real), intent(in) :: roc, undulat
231 
232  real(kind=wp) :: r8lat
233  integer :: i
234  real(kind=kind_real) :: rlon_local
235 
236 !-------------------------------------------------------------------------
237  r8lat = real(rlat,kind=wp)
238  y%lat = r8lat
239  rlon_local = rlon
240  if (rlon_local .gt. 180.) rlon_local = rlon_local - 360.0 ! ROPP Longitude value -180 to 180
241  y%lon = real(rlon_local,kind=wp)
242  y%g_sfc = gravity(r8lat, 0.0_wp) ! 2nd argument is height(m) above the sfc
243  y%r_curve = real(roc,kind=wp) ! ROPP rad of curve (m)
244  y%undulation = real(undulat,kind=wp) ! ROPP undulation corr for rad of curve (m)
245  y%r_earth = r_eff(r8lat)
246 
247 !---------------------------------------------------
248 ! allocate bending angle, impact parameter & weights
249 !----------------------------------------------------
250  allocate(y%bangle(1:nvprof)) ! value computed in fwd model
251  allocate(y%impact(1:nvprof))
252 ! needed for 2D
253  allocate(y%a_path(1:nvprof,2))
254  allocate(y%rtan(1:nvprof))
255 
256 ! number of points in profile
257  y%nobs = nvprof
258 
259  do i=1,nvprof
260  y%impact(i) = real(obs_impact(i),kind=wp) ! ROPP expects impact parameter in meters
261  end do
262 
263  return
264 
265 end subroutine init_ropp_2d_obvec
266 
267 !-------------------------------------------------------------------------
268 !-------------------------------------------------------------------------
269 subroutine init_ropp_2d_obvec_tlad(iloop,nvprof,obs_impact, &
270  rlat,rlon,roc,undulat,y,y_p)
271 
272 ! Description:
273 ! subroutine to fill a ROPP observation vector structure
274 ! observation provides the inputs of:
275 ! impact parameter, lat, lon, time, radius of curvature
276 !
277 ! forward model will provide the simulation of bending angle
278 !
279 ! Inputs:
280 ! obs_impact impact parameter
281 ! rlat latitude
282 ! rlon longitude
283 ! roc radius of curvature
284 ! undulat undulation correction for radius of curvature
285 !
286 ! Outputs:
287 ! y: Partially filled Forward model observation vector
288 !-------------------------------------------------------------------------
289  implicit none
290 ! Output state vector
291  type(obs1dbangle), intent(out) :: y,y_p
292 
293  integer, intent(in) :: iloop
294  integer, intent(in) :: nvprof
295  real(kind=kind_real), dimension(nvprof), intent(in) :: obs_impact
296  real(kind=kind_real), intent(in) :: rlat,rlon
297  real(kind=kind_real), intent(in) :: roc, undulat
298 
299  real(kind=wp) :: r8lat
300  integer :: i
301  real(kind=kind_real) :: rlon_local
302 
303 !-------------------------------------------------------------------------
304  r8lat = real(rlat,kind=wp)
305  y%lat = r8lat
306  rlon_local = rlon
307  if (rlon_local .gt. 180.) rlon_local = rlon_local - 360.0 ! ROPP Longitude value -180 to 180
308  y%lon = real(rlon_local,kind=wp)
309  y%g_sfc = gravity(r8lat, 0.0_wp) ! 2nd argument is height(m) above the sfc
310  y%r_curve = real(roc,kind=wp) ! ROPP rad of curve (m)
311  y%undulation = real(undulat,kind=wp) ! ROPP undulation corr for rad of curve (m)
312  y%r_earth = r_eff(r8lat)
313 
314 !--------------------------------------------------------
315 ! allocate bending angle, impact parameter & weights
316 !---------------------------------------------------------
317 ! if (iloop ==1) then
318 
319  allocate(y%bangle(1:nvprof)) ! value computed in fwd model
320  allocate(y%impact(1:nvprof))
321 ! needed for 2D
322  allocate(y%a_path(1:nvprof,2))
323  allocate(y%rtan(1:nvprof))
324 ! TL code
325  allocate(y_p%bangle(1:nvprof)) ! value computed in fwd model
326  allocate(y_p%impact(1:nvprof))
327 ! needed for 2D
328  allocate(y_p%a_path(1:nvprof,2))
329  allocate(y_p%rtan(1:nvprof))
330 
331 ! endif
332 
333 
334 ! number of points in profile
335  y%nobs = nvprof
336  y_p%nobs = nvprof
337 
338  do i = 1 ,nvprof
339  y%impact(i) = real(obs_impact(i),kind=wp) ! ROPP expects impact parameter in meters
340  y_p%impact(i) = real(obs_impact(i),kind=wp) ! ROPP expects impact parameter in meters
341  end do
342 
343  return
344 
345 end subroutine init_ropp_2d_obvec_tlad
346 
347 
348 !-------------------------------------------------------------------------
349 !-------------------------------------------------------------------------
350 subroutine ropp_tidy_up_2d(x,y)
351  implicit none
352  type(state2dfm), intent(inout) :: x
353  type(obs1dbangle), intent(inout) :: y
354 ! x
355  if (associated(x%temp)) deallocate(x%temp)
356  if (associated(x%shum)) deallocate(x%shum)
357  if (associated(x%pres)) deallocate(x%pres)
358  if (associated(x%geop)) deallocate(x%geop)
359  if (associated(x%nr)) deallocate(x%nr)
360  if (associated(x%refrac)) deallocate(x%refrac)
361  if (associated(x%lat)) deallocate(x%lat)
362  if (associated(x%lon)) deallocate(x%lon)
363 ! y
364  if (associated(y%impact)) deallocate(y%impact)
365  if (associated(y%bangle)) deallocate(y%bangle)
366  if (associated(y%a_path)) deallocate(y%a_path)
367  if (associated(y%rtan)) deallocate(y%rtan)
368 
369 end subroutine ropp_tidy_up_2d
370 
371 !-------------------------------------------------------------------------
372 !-------------------------------------------------------------------------
373 subroutine ropp_tidy_up_tlad_2d(x,x_p,y,y_p)
374  implicit none
375  type(state2dfm), intent(inout) :: x,x_p ! x_p can be either x_tl or x_ad
376  type(obs1dbangle), intent(inout) :: y,y_p ! y_p can be either y_tl or y_ad
377 ! x
378  deallocate(x%temp)
379  deallocate(x%shum)
380  deallocate(x%pres)
381  deallocate(x%geop)
382  deallocate(x%nr)
383  deallocate(x%refrac)
384  deallocate(x%lat)
385  deallocate(x%lon)
386  deallocate(x_p%temp)
387  deallocate(x_p%shum)
388  deallocate(x_p%pres)
389  deallocate(x_p%geop)
390  deallocate(x_p%nr)
391  deallocate(x_p%refrac)
392  deallocate(x_p%lat)
393  deallocate(x_p%lon)
394 ! y
395  deallocate(y%impact)
396  deallocate(y%bangle)
397  deallocate(y%a_path)
398  deallocate(y%rtan)
399  deallocate(y_p%impact)
400  deallocate(y_p%bangle)
401  deallocate(y_p%a_path)
402  deallocate(y_p%rtan)
403 
404  return
405 
406 end subroutine ropp_tidy_up_tlad_2d
407 
408 end module ufo_gnssro_ropp2d_utils_mod
ufo_gnssro_ropp2d_utils_mod
Fortran module to handle gnssro bending angle observations following the ROPP (2018 Aug) implementati...
Definition: ufo_gnssro_ropp2d_utils_mod.F90:9
ufo_gnssro_ropp2d_utils_mod::ropp_tidy_up_2d
subroutine, public ropp_tidy_up_2d(x, y)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:351
ufo_gnssro_ropp2d_utils_mod::init_ropp_2d_statevec
subroutine, public init_ropp_2d_statevec(rlon, rlat, temp, shum, pres, phi, lm, x, n_horiz, dtheta, iflip)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:36
ufo_gnssro_ropp2d_utils_mod::ropp_tidy_up_tlad_2d
subroutine, public ropp_tidy_up_tlad_2d(x, x_p, y, y_p)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:374
ufo_gnssro_ropp2d_utils_mod::init_ropp_2d_obvec
subroutine, public init_ropp_2d_obvec(nvprof, obs_impact, rlat, rlon, roc, undulat, y)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:204
ufo_gnssro_ropp2d_utils_mod::init_ropp_2d_obvec_tlad
subroutine, public init_ropp_2d_obvec_tlad(iloop, nvprof, obs_impact, rlat, rlon, roc, undulat, y, y_p)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:271
ufo_gnssro_ropp2d_utils_mod::init_ropp_2d_statevec_ad
subroutine, public init_ropp_2d_statevec_ad(temp_d, shum_d, pres_d, phi_d, lm, x_ad, n_horiz, iflip)
Definition: ufo_gnssro_ropp2d_utils_mod.F90:125