10 use oops_variables_mod
12 use missing_values_mod
24 type(oops_variables),
public :: obsvars
25 type(oops_variables),
public :: geovars
26 character(len=MAXVARLEN) :: da_psfc_scheme
38 use fckit_configuration_module,
only: fckit_configuration
41 type(fckit_configuration),
intent(in) :: f_conf
42 character(len=:),
allocatable :: str
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
59 use fckit_log_module,
only : fckit_log
62 integer,
intent(in) :: nvars, nlocs
64 real(c_double),
intent(inout) :: hofx(nvars, nlocs)
65 type(c_ptr),
value,
intent(in) :: obss
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
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
81 missing = missing_value(missing)
82 nobs = obsspace_get_nlocs(obss)
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)
91 allocate(cor_psfc(nobs))
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)
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)
112 allocate(model_zs(nobs))
113 allocate(model_level1(nobs))
114 allocate(model_psfc(nobs))
116 model_zs = model_sfc_geomz%vals(1,:)
117 model_level1 = model_geomz%vals(model_geomz%nval,:)
118 model_psfc = model_ps%vals(1,:)
121 select case (trim(self%da_psfc_scheme))
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)
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)
137 allocate(model_tvs(nobs))
138 model_tvs = model_tv%vals(model_tv%nval,:) +
lclr * ( model_level1 - model_zs )
141 call da_intpsfc_prs(nobs, missing, cor_psfc, obs_height, obs_psfc, model_zs, model_tvs, obs_t, obs_q)
145 deallocate(model_tvs)
149 allocate(model_p_2000(nobs))
150 allocate(model_tv_2000(nobs))
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)
159 call da_intpsfc_prs_ukmo(nobs, missing, cor_psfc, obs_height, obs_psfc, model_zs, model_psfc, model_tv_2000, model_p_2000)
161 deallocate(model_p_2000)
162 deallocate(model_tv_2000)
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)
173 if ( cor_psfc(iobs) /= missing)
then
174 hofx(ivar,iobs) = obs_psfc(iobs) - cor_psfc(iobs) + model_psfc(iobs)
176 hofx(ivar,iobs) = model_psfc(iobs)
181 deallocate(obs_height)
185 deallocate(model_level1)
186 deallocate(model_psfc)
214 integer,
intent (in) :: nobs
215 real(c_double),
intent (in) :: missing
216 real(kind_real),
dimension(nobs),
intent (out) :: P_o2m
217 real(kind_real),
dimension(nobs),
intent (in) :: H_o, P_o
218 real(kind_real),
dimension(nobs),
intent (in) :: H_m, TV_m
219 real(kind_real),
dimension(nobs),
intent (in),
optional :: T_o, Q_o
220 real(kind_real),
dimension(nobs) :: TV_o, TV
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
233 tv = 0.5 * (tv_m + tv_o)
238 where ( h_o /= missing .and. p_o /= missing )
239 p_o2m = p_o * exp( - (h_m - h_o) *
grav / (
rd * tv) )
280 integer,
intent (in) :: nobs
281 real(c_double),
intent (in) :: missing
282 real(kind_real),
dimension(nobs),
intent (out) :: P_o2m
283 real(kind_real),
dimension(nobs),
intent (in) :: H_o, P_o
284 real(kind_real),
dimension(nobs),
intent (in) :: H_m, P_m, TV_2000, P_2000
285 real(kind_real) :: ind
286 real(kind_real),
dimension(nobs) :: P_m2o, T_m, T_m2o
292 where ( h_o /= missing .and. p_o /= missing )
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)
300 p_o2m = p_o * p_m / p_m2o