UFO
ufo_sfcpcorrected_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 sfcpcorrected observation operator
7 
9 
10  use oops_variables_mod
11  use ufo_vars_mod
12  use missing_values_mod
13  use iso_c_binding
14  use kinds
15  use ufo_constants_mod, only : grav, rd, lclr, t2tv
16 
17  implicit none
18  private
19  integer, parameter :: max_string = 800
20 
21 !> Fortran derived type for the observation type
22  type, public :: ufo_sfcpcorrected
23  private
24  type(oops_variables), public :: obsvars
25  type(oops_variables), public :: geovars
26  character(len=MAXVARLEN) :: da_psfc_scheme
27  contains
28  procedure :: setup => ufo_sfcpcorrected_setup
29  procedure :: simobs => ufo_sfcpcorrected_simobs
30  end type ufo_sfcpcorrected
31 
32  character(len=MAXVARLEN), dimension(5), parameter :: geovars_default = (/ var_ps, var_geomz, var_sfc_geomz, var_tv, var_prs /)
33 
34 contains
35 
36 ! ------------------------------------------------------------------------------
37 subroutine ufo_sfcpcorrected_setup(self, f_conf)
38 use fckit_configuration_module, only: fckit_configuration
39 implicit none
40 class(ufo_sfcpcorrected), intent(inout) :: self
41 type(fckit_configuration), intent(in) :: f_conf
42 character(len=:), allocatable :: str
43 
44 call self%geovars%push_back(geovars_default)
45 
46 self%da_psfc_scheme = "UKMO"
47 if (f_conf%has("da_psfc_scheme")) then
48  call f_conf%get_or_die("da_psfc_scheme",str)
49  self%da_psfc_scheme = str
50 end if
51 
52 end subroutine ufo_sfcpcorrected_setup
53 
54 ! ------------------------------------------------------------------------------
55 subroutine ufo_sfcpcorrected_simobs(self, geovals, obss, nvars, nlocs, hofx)
57 use obsspace_mod
59 use fckit_log_module, only : fckit_log
60 implicit none
61 class(ufo_sfcpcorrected), intent(in) :: self
62 integer, intent(in) :: nvars, nlocs
63 type(ufo_geovals), intent(in) :: geovals
64 real(c_double), intent(inout) :: hofx(nvars, nlocs)
65 type(c_ptr), value, intent(in) :: obss
66 
67 ! Local variables
68 real(c_double) :: missing
69 real(kind_real) :: H2000 = 2000.0
70 integer :: nobs, iobs, ivar
71 real(kind_real), allocatable :: cor_psfc(:)
72 type(ufo_geoval), pointer :: model_ps, model_p, model_sfc_geomz, model_tv, model_geomz
73 character(len=*), parameter :: myname_="ufo_sfcpcorrected_simobs"
74 character(max_string) :: err_msg
75 real(kind_real) :: wf
76 integer :: wi
77 logical :: variable_present
78 real(kind_real), dimension(:), allocatable :: obs_height, obs_t, obs_q, obs_psfc
79 real(kind_real), dimension(:), allocatable :: model_tvs, model_zs, model_level1, model_p_2000, model_tv_2000, model_psfc
80 
81 missing = missing_value(missing)
82 nobs = obsspace_get_nlocs(obss)
83 
84 ! check if nobs is consistent in geovals & nlocs
85 if (geovals%nlocs /= nobs) then
86  write(err_msg,*) myname_, ' error: nlocs of model and obs is inconsistent!'
87  call abor1_ftn(err_msg)
88 endif
89 
90 ! cor_psfc: observed surface pressure at model surface height, corresponding to P_o2m in da_intpsfc_prs* subroutines
91 allocate(cor_psfc(nobs))
92 
93 ! get obs variables
94 allocate(obs_height(nobs))
95 allocate(obs_psfc(nobs))
96 call obsspace_get_db(obss, "MetaData", "station_elevation",obs_height)
97 call obsspace_get_db(obss, "ObsValue", "surface_pressure", obs_psfc)
98 
99 ! get model variables
100 call ufo_geovals_get_var(geovals, var_ps, model_ps)
101 call ufo_geovals_get_var(geovals, var_geomz, model_geomz)
102 call ufo_geovals_get_var(geovals, var_sfc_geomz, model_sfc_geomz)
103 call ufo_geovals_get_var(geovals, var_tv, model_tv)
104 call ufo_geovals_get_var(geovals, var_prs, model_p)
105 
106 if (model_geomz%vals(1,1) .gt. model_geomz%vals(model_geomz%nval,1) ) then
107  write(err_msg,'(a)') ' ufo_sfcpcorrected:'//new_line('a')// &
108  ' Model vertical height profile is from top to bottom'
109  call fckit_log%info(err_msg)
110 end if
111 
112 allocate(model_zs(nobs))
113 allocate(model_level1(nobs))
114 allocate(model_psfc(nobs))
115 
116 model_zs = model_sfc_geomz%vals(1,:)
117 model_level1 = model_geomz%vals(model_geomz%nval,:) !reverse
118 model_psfc = model_ps%vals(1,:)
119 
120 ! do terrain height correction, two optional schemes
121 select case (trim(self%da_psfc_scheme))
122 
123 case ("WRFDA")
124  ! get extra obs values
125  variable_present = obsspace_has(obss, "ObsValue", "air_temperature")
126  if (variable_present) then
127  allocate(obs_t(nobs))
128  call obsspace_get_db(obss, "ObsValue", "air_temperature", obs_t)
129  end if
130  variable_present = obsspace_has(obss, "ObsValue", "specific_humidity")
131  if (variable_present) then
132  allocate(obs_q(nobs))
133  call obsspace_get_db(obss, "ObsValue", "specific_humidity", obs_q)
134  end if
135 
136  ! get extra model values
137  allocate(model_tvs(nobs))
138  model_tvs = model_tv%vals(model_tv%nval,:) + lclr * ( model_level1 - model_zs ) !Lclr = 0.0065 K/m
139 
140  ! correction
141  call da_intpsfc_prs(nobs, missing, cor_psfc, obs_height, obs_psfc, model_zs, model_tvs, obs_t, obs_q)
142 
143  deallocate(obs_t)
144  deallocate(obs_q)
145  deallocate(model_tvs)
146 
147 case ("UKMO")
148 
149  allocate(model_p_2000(nobs))
150  allocate(model_tv_2000(nobs))
151  do iobs = 1, nobs
152  ! vertical interpolation for getting model P and tv at 2000 m
153  call vert_interp_weights(model_geomz%nval, h2000, model_geomz%vals(:,iobs), wi, wf)
154  call vert_interp_apply(model_p%nval, model_p%vals(:,iobs), model_p_2000(iobs), wi, wf)
155  call vert_interp_apply(model_tv%nval, model_tv%vals(:,iobs), model_tv_2000(iobs), wi, wf)
156  end do
157 
158  ! correction
159  call da_intpsfc_prs_ukmo(nobs, missing, cor_psfc, obs_height, obs_psfc, model_zs, model_psfc, model_tv_2000, model_p_2000)
160 
161  deallocate(model_p_2000)
162  deallocate(model_tv_2000)
163 
164 case default
165  write(err_msg,*) "ufo_sfcpcorrected_mod.F90: da_psfc_scheme must be WRFDA or UKMO"
166  call fckit_log%info(err_msg)
167  call abor1_ftn(err_msg)
168 end select
169 
170 ! update the obs surface pressure
171 do ivar = 1, nvars
172  do iobs = 1, nlocs
173  if ( cor_psfc(iobs) /= missing) then
174  hofx(ivar,iobs) = obs_psfc(iobs) - cor_psfc(iobs) + model_psfc(iobs)
175  else
176  hofx(ivar,iobs) = model_psfc(iobs)
177  end if
178  enddo
179 enddo
180 
181 deallocate(obs_height)
182 deallocate(obs_psfc)
183 
184 deallocate(model_zs)
185 deallocate(model_level1)
186 deallocate(model_psfc)
187 
188 end subroutine ufo_sfcpcorrected_simobs
189 
190 ! ------------------------------------------------------------------------------
191 !> \Conduct terrain height correction for surface pressure
192 !!
193 !! \This subroutine is based on a subroutine from WRFDA da_intpsfc_prs.inc file
194 !! corresponding to sfc_assi_options = 1 in WRFDA's namelist
195 !!
196 !! \Date: June 2019: Created
197 !! \Method: hydrosatic equation
198 !!
199 !! P_o2m = P_o * exp [-grav/rd * (H_m-H_o) / (TV_m + TV_o)/2)
200 !!
201 !! Where:
202 !! H_m = model surface height
203 !! H_o = station height
204 !! TV_m = virtual temperature at model surface height
205 !! TV_o = virtual temperature at station height
206 !! P_o2m = pressure interpolated from station height to model surface height
207 !! P_o = pressure at station height
208 !! grav = gravitational acceleration
209 !! rd = gas constant per mole
210 !!
211 
212 subroutine da_intpsfc_prs (nobs, missing, P_o2m, H_o, P_o, H_m, TV_m, T_o, Q_o)
213 implicit none
214 integer, intent (in) :: nobs !<total observation number
215 real(c_double), intent (in) :: missing
216 real(kind_real), dimension(nobs), intent (out) :: P_o2m !<observed PS at model sfc height
217 real(kind_real), dimension(nobs), intent (in) :: H_o, P_o !<observed Height and PS
218 real(kind_real), dimension(nobs), intent (in) :: H_m, TV_m !<model sfc height and TV
219 real(kind_real), dimension(nobs), intent (in), optional :: T_o, Q_o !<obserbed T and Q
220 real(kind_real), dimension(nobs) :: TV_o, TV
221 
222 ! 1. model and observation virtual temperature
223 ! ---------------------------------------------
224 
225 if (present(t_o) .and. present(q_o)) then
226  tv_o = t_o * (1.0 + t2tv * q_o)
227 else if (present(t_o) .and. .not.(present(q_o))) then
228  tv_o = t_o
229 else
230  tv_o = tv_m
231 end if
232 
233 tv = 0.5 * (tv_m + tv_o)
234 
235 ! 2. extrapolate pressure from station height to model surface height
236 ! --------------------------------------------------------------------
237 
238 where ( h_o /= missing .and. p_o /= missing )
239  p_o2m = p_o * exp( - (h_m - h_o) * grav / (rd * tv) )
240 elsewhere
241  p_o2m = p_o
242 end where
243 
244 end subroutine da_intpsfc_prs
245 
246 ! ------------------------------------------------------------------------------
247 !> \Conduct terrain height correction for surface pressure
248 !!
249 !! \Reference: Ingleby,2013. UKMO Technical Report No: 582. Appendix 1.
250 !!
251 !! \Method: integrate the hydrosatic equation dp/dz=-rho*g/RT to get P_m2o first, equation:
252 !!
253 !! (P_m2o/P_m)=(T_m2o/T_m)** (grav/rd*L)
254 !!
255 !! Where:
256 !! P_m2o = model surface pressure at station height
257 !! P_m = model surface pressure
258 !! T_m = temperature at model surface height; derived from TV_2000
259 !! T_m2o = model surface temperature at station height
260 !! grav = gravitational acceleration
261 !! rd = gas constant per mole
262 !! Lclr = constant lapse rate (0.0065 K/m)
263 !!
264 !! To avoid dirunal/local variations, use TV_2000 (2000 m above the model surface height) instead of direct T_m
265 !!
266 !! T_m = TV_2000 * (P_o / P_2000) ** (rd*L/grav)
267 !!
268 !! Where:
269 !! P_2000 = background pressure at 2000 m
270 !! TV_2000 = background virtual temperature at 2000 m
271 !! P_o = pressure at station height
272 !!
273 !! Finally, in practice, adjust P_o to the model surface height using
274 !!
275 !! P_o2m = P_o * (P_m / P_m2o)
276 !!
277 
278 subroutine da_intpsfc_prs_ukmo (nobs, missing, P_o2m, H_o, P_o, H_m, P_m, TV_2000, P_2000)
279 implicit none
280 integer, intent (in) :: nobs !<total observation number
281 real(c_double), intent (in) :: missing
282 real(kind_real), dimension(nobs), intent (out) :: P_o2m !<observed PS at model sfc height
283 real(kind_real), dimension(nobs), intent (in) :: H_o, P_o !<observed Height and PS
284 real(kind_real), dimension(nobs), intent (in) :: H_m, P_m, TV_2000, P_2000 !<model Height, PS, TV at 2000 m, and P at 2000 m
285 real(kind_real) :: ind !<local variable
286 real(kind_real), dimension(nobs) :: P_m2o, T_m, T_m2o !<local variables: model PS at observed height, model sfc T,
287  !! and model T at observed height
288 
289 ! define the constant power exponent
290 ind = rd * lclr / grav
291 
292 where ( h_o /= missing .and. p_o /= missing )
293  ! calculate T_m -- background temperature at model surface height
294  ! T_m2o -- background temperature at station height
295  t_m = tv_2000 * (p_m / p_2000) ** ind
296  t_m2o = t_m + lclr * ( h_m - h_o)
297  p_m2o = p_m * (t_m2o / t_m) ** (1.0 / ind)
298 
299  ! In practice, P_o is adjusted to the model surface height
300  p_o2m = p_o * p_m / p_m2o
301 elsewhere
302  p_o2m = p_o
303 end where
304 
305 end subroutine da_intpsfc_prs_ukmo
306 
307 ! ------------------------------------------------------------------------------
308 
309 end module ufo_sfcpcorrected_mod
ufo_avgkernel_mod::max_string
integer, parameter max_string
Definition: ufo_avgkernel_mod.F90:17
ufo_constants_mod::grav
real(kind_real), parameter, public grav
Definition: ufo_constants_mod.F90:9
ufo_sfcpcorrected_mod::ufo_sfcpcorrected_setup
subroutine ufo_sfcpcorrected_setup(self, f_conf)
Definition: ufo_sfcpcorrected_mod.F90:38
ufo_radarradialvelocity_mod::geovars_default
character(len=maxvarlen), dimension(3), parameter geovars_default
Definition: ufo_radarradialvelocity_mod.F90:27
ufo_constants_mod::t2tv
real(kind_real), parameter, public t2tv
Definition: ufo_constants_mod.F90:46
ufo_vars_mod::var_ps
character(len=maxvarlen), parameter, public var_ps
Definition: ufo_variables_mod.F90:28
ufo_sfcpcorrected_mod::da_intpsfc_prs_ukmo
subroutine da_intpsfc_prs_ukmo(nobs, missing, P_o2m, H_o, P_o, H_m, P_m, TV_2000, P_2000)
\Conduct terrain height correction for surface pressure
Definition: ufo_sfcpcorrected_mod.F90:279
ufo_vars_mod::var_geomz
character(len=maxvarlen), parameter, public var_geomz
Definition: ufo_variables_mod.F90:69
ufo_vars_mod::var_tv
character(len=maxvarlen), parameter, public var_tv
Definition: ufo_variables_mod.F90:18
vert_interp_mod
Fortran module to perform linear interpolation.
Definition: vert_interp.F90:8
vert_interp_mod::vert_interp_weights
subroutine vert_interp_weights(nlev, obl, vec, wi, wf)
Definition: vert_interp.F90:22
ufo_geovals_mod
Definition: ufo_geovals_mod.F90:7
ufo_sfcpcorrected_mod::ufo_sfcpcorrected
Fortran derived type for the observation type.
Definition: ufo_sfcpcorrected_mod.F90:22
ufo_sfcpcorrected_mod
Fortran module for sfcpcorrected observation operator.
Definition: ufo_sfcpcorrected_mod.F90:8
ufo_constants_mod
Definition: ufo_constants_mod.F90:2
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_constants_mod::lclr
real(kind_real), parameter, public lclr
Definition: ufo_constants_mod.F90:45
vert_interp_mod::vert_interp_apply
subroutine vert_interp_apply(nlev, fvec, f, wi, wf)
Definition: vert_interp.F90:73
ufo_geovals_mod::ufo_geovals_get_var
subroutine, public ufo_geovals_get_var(self, varname, geoval)
Definition: ufo_geovals_mod.F90:128
ufo_sfcpcorrected_mod::da_intpsfc_prs
subroutine da_intpsfc_prs(nobs, missing, P_o2m, H_o, P_o, H_m, TV_m, T_o, Q_o)
\Conduct terrain height correction for surface pressure
Definition: ufo_sfcpcorrected_mod.F90:213
ufo_sfcpcorrected_mod::ufo_sfcpcorrected_simobs
subroutine ufo_sfcpcorrected_simobs(self, geovals, obss, nvars, nlocs, hofx)
Definition: ufo_sfcpcorrected_mod.F90:56
ufo_geovals_mod::ufo_geovals
type to hold interpolated fields required by the obs operators
Definition: ufo_geovals_mod.F90:47
ufo_constants_mod::rd
real(kind_real), parameter, public rd
Definition: ufo_constants_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_vars_mod::var_prs
character(len=maxvarlen), parameter, public var_prs
Definition: ufo_variables_mod.F90:25