OOPS
qg_obsvec_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2009-2016 ECMWF.
2 ! (C) Copyright 2017-2019 UCAR.
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 ! In applying this licence, ECMWF does not waive the privileges and immunities
7 ! granted to it by virtue of its status as an intergovernmental organisation nor
8 ! does it submit to any jurisdiction.
9 
11 
12 use iso_c_binding
13 use kinds
14 use random_mod
15 
16 implicit none
17 
18 private
19 public :: qg_obsvec
20 public :: qg_obsvec_registry
24 ! ------------------------------------------------------------------------------
25 interface
26  subroutine qg_obsvec_random_i(odb,nn,zz) bind(c,name='qg_obsvec_random_f')
27  use iso_c_binding
28  implicit none
29  type(c_ptr),intent(in) :: odb
30  integer(c_int),intent(in) :: nn
31  real(c_double),intent(inout) :: zz
32  end subroutine qg_obsvec_random_i
33 end interface
34 ! ------------------------------------------------------------------------------
36  integer :: nobs = 0 !< Number of observations
37  integer :: nlev = 0 !< Number of levels
38  real(kind_real),allocatable :: values(:,:) !< Values
39 end type qg_obsvec
40 
41 #define LISTED_TYPE qg_obsvec
42 
43 !> Linked list interface - defines registry_t type
44 #include "oops/util/linkedList_i.f"
45 
46 !> Global registry
47 type(registry_t) :: qg_obsvec_registry
48 ! ------------------------------------------------------------------------------
49 contains
50 ! ------------------------------------------------------------------------------
51 ! Public
52 ! ------------------------------------------------------------------------------
53 !> Linked list implementation
54 #include "oops/util/linkedList_c.f"
55 ! ------------------------------------------------------------------------------
56 !> Setup observation vector
57 subroutine qg_obsvec_setup(self,nlev,nobs)
58 
59 implicit none
60 
61 ! Passed variables
62 type(qg_obsvec),intent(inout) :: self !< Observation vector
63 integer,intent(in) :: nlev !< Number of levels
64 integer,intent(in) :: nobs !< Number of observations
65 
66 ! Set sizes
67 self%nlev = nlev
68 self%nobs = nobs
69 
70 ! Release memory
71 if (allocated(self%values)) deallocate(self%values)
72 
73 ! Allocation
74 allocate(self%values(self%nlev,self%nobs))
75 
76 ! Initialization
77 self%values = 0.0_kind_real
78 
79 end subroutine qg_obsvec_setup
80 ! ------------------------------------------------------------------------------
81 !> Clone observation vector
82 subroutine qg_obsvec_clone(self,other)
83 
84 implicit none
85 
86 ! Passed variables
87 type(qg_obsvec),intent(inout) :: self !< Observation vector
88 type(qg_obsvec),intent(in) :: other !< Other observation vector
89 
90 ! Set sizes
91 self%nlev = other%nlev
92 self%nobs = other%nobs
93 
94 ! Allocation
95 allocate(self%values(self%nlev,self%nobs))
96 
97 end subroutine qg_obsvec_clone
98 ! ------------------------------------------------------------------------------
99 !> Delete observation vector
100 subroutine qg_obsvec_delete(self)
101 
102 implicit none
103 
104 ! Passed variables
105 type(qg_obsvec),intent(inout) :: self !< Observation vector
106 
107 ! Release memory
108 deallocate(self%values)
109 
110 end subroutine qg_obsvec_delete
111 ! ------------------------------------------------------------------------------
112 !> Copy observation vector
113 subroutine qg_obsvec_copy(self,other)
114 
115 implicit none
116 
117 ! Passed variables
118 type(qg_obsvec),intent(inout) :: self !< Observation vector
119 type(qg_obsvec),intent(in) :: other !< Other observation vector
120 
121 if ((other%nlev/=self%nlev).or.(other%nobs/=self%nobs)) then
122  ! Release memory
123  deallocate(self%values)
124 
125  ! Set sizes
126  self%nlev = other%nlev
127  self%nobs = other%nobs
128 
129  ! Allocation
130  allocate(self%values(self%nlev,self%nobs))
131 endif
132 
133 ! Copy data
134 self%values = other%values
135 
136 end subroutine qg_obsvec_copy
137 ! ------------------------------------------------------------------------------
138 !> Copy a local subset of the observation vector
139 subroutine qg_obsvec_copy_local(self,other,idx)
140 
141 implicit none
142 
143 ! Passed variables
144 type(qg_obsvec),intent(inout) :: self !< Observation vector
145 type(qg_obsvec),intent(in) :: other !< Other observation vector
146 integer,intent(in) :: idx(:)
147 
148 ! local variables
149 integer :: i
150 
151 if ((other%nlev/=self%nlev).or.(size(idx)/=self%nobs)) then
152  ! Release memory
153  deallocate(self%values)
154 
155  ! Set sizes
156  self%nlev = other%nlev
157  self%nobs = size(idx)
158 
159  ! Allocation
160  allocate(self%values(self%nlev,self%nobs))
161 endif
162 
163 ! Copy data
164 do i = 1,self%nobs
165  self%values(:,i) = other%values(:,idx(i)+1)
166 enddo
167 
168 end subroutine qg_obsvec_copy_local
169 ! ------------------------------------------------------------------------------
170 !> Set observation vector to zero
171 subroutine qg_obsvec_zero(self)
172 
173 implicit none
174 
175 ! Passed variables
176 type(qg_obsvec),intent(inout) :: self !< Observation vector
177 
178 ! Set observation vector to zero
179 self%values = 0.0
180 
181 end subroutine qg_obsvec_zero
182 ! ------------------------------------------------------------------------------
183 !> Multiply observation vector with a scalar
184 subroutine qg_obsvec_mul_scal(self,zz)
185 
186 implicit none
187 
188 ! Passed variables
189 type(qg_obsvec),intent(inout) :: self !< Observation vector
190 real(kind_real),intent(in) :: zz !< Multiplier
191 
192 ! Multiply observation vector with a scalar
193 self%values = zz*self%values
194 
195 end subroutine qg_obsvec_mul_scal
196 ! ------------------------------------------------------------------------------
197 !> Add observation vector
198 subroutine qg_obsvec_add(self,other)
199 
200 implicit none
201 
202 ! Passed variables
203 type(qg_obsvec),intent(inout) :: self !< Observation vector
204 type(qg_obsvec),intent(in) :: other !< Other observation vector
205 
206 ! Add observation vector
207 self%values = self%values+other%values
208 
209 end subroutine qg_obsvec_add
210 ! ------------------------------------------------------------------------------
211 !> Subtract observation vector
212 subroutine qg_obsvec_sub(self,other)
213 
214 implicit none
215 
216 ! Passed variables
217 type(qg_obsvec),intent(inout) :: self !< Observation vector
218 type(qg_obsvec),intent(in) :: other !< Other observation vector
219 
220 ! Subtract observation vector
221 self%values = self%values-other%values
222 
223 end subroutine qg_obsvec_sub
224 ! ------------------------------------------------------------------------------
225 !> Multiply observation vector
226 subroutine qg_obsvec_mul(self,other)
227 
228 implicit none
229 
230 ! Passed variables
231 type(qg_obsvec),intent(inout) :: self !< Observation vector
232 type(qg_obsvec),intent(in) :: other !< Other observation vector
233 
234 ! Multiply observation vector
235 self%values = self%values*other%values
236 
237 end subroutine qg_obsvec_mul
238 ! ------------------------------------------------------------------------------
239 !> Divide observation vector
240 subroutine qg_obsvec_div(self,other)
241 
242 implicit none
243 
244 ! Passed variables
245 type(qg_obsvec),intent(inout) :: self !< Observation vector
246 type(qg_obsvec),intent(in) :: other !< Other observation vector
247 
248 ! Divide observation vector
249 self%values = self%values/other%values
250 
251 end subroutine qg_obsvec_div
252 ! ------------------------------------------------------------------------------
253 !> Apply axpy on observation vector
254 subroutine qg_obsvec_axpy(self,zz,other)
255 
256 implicit none
257 
258 ! Passed variables
259 type(qg_obsvec),intent(inout) :: self !< Observation vector
260 real(kind_real),intent(in) :: zz !< Multiplier
261 type(qg_obsvec),intent(in) :: other !< Other observation vector
262 
263 ! Apply axpy on observation vector
264 self%values = self%values+zz*other%values
265 
266 end subroutine qg_obsvec_axpy
267 ! ------------------------------------------------------------------------------
268 !> Invert observation vector
269 subroutine qg_obsvec_invert(self)
270 
271 implicit none
272 
273 ! Passed variables
274 type(qg_obsvec),intent(inout) :: self !< Observation vector
275 
276 ! Invert observation vector
277 self%values = 1.0/self%values
278 
279 end subroutine qg_obsvec_invert
280 ! ------------------------------------------------------------------------------
281 !> Generate random observation vector
282 subroutine qg_obsvec_random(c_odb,self)
283 
284 implicit none
285 
286 ! Passed variables
287 type(c_ptr),intent(in) :: c_odb !< Observation data base
288 type(qg_obsvec),intent(inout) :: self !< Observation vector
289 
290 ! Local variables
291 integer :: nval
292 
293 ! Compute total size
294 nval = self%nobs*self%nlev
295 
296 ! Get random values
297 call qg_obsvec_random_i(c_odb,nval,self%values(1,1))
298 
299 end subroutine qg_obsvec_random
300 ! ------------------------------------------------------------------------------
301 !> Compute dot product between observation vectors
302 subroutine qg_obsvec_dotprod(obsvec1,obsvec2,zz)
303 
304 implicit none
305 
306 ! Passed variables
307 type(qg_obsvec),intent(in) :: obsvec1 !< Observation vector 1
308 type(qg_obsvec),intent(in) :: obsvec2 !< Observation vector 2
309 real(kind_real),intent(inout) :: zz !< Dot product
310 
311 ! Local variables
312 integer :: jlev,jobs
313 
314 ! Check sizes
315 if ((obsvec1%nobs/=obsvec2%nobs).or.(obsvec1%nlev/=obsvec2%nlev)) call abor1_ftn('qg_obsvec_dotprod: inconsistent sizes')
316 
317 ! Initalization
318 zz = 0.0
319 
320 ! Loop over values
321 do jobs=1,obsvec1%nobs
322  do jlev=1,obsvec1%nlev
323  zz = zz+obsvec1%values(jlev,jobs)*obsvec2%values(jlev,jobs)
324  enddo
325 enddo
326 
327 end subroutine qg_obsvec_dotprod
328 ! ------------------------------------------------------------------------------
329 !> Compute observation vector statistics
330 subroutine qg_obsvec_stats(self,scaling,zmin,zmax,zavg)
331 
332 implicit none
333 
334 ! Passed variables
335 type(qg_obsvec),intent(in):: self !< Observation vector
336 real(kind_real),intent(inout) :: scaling !< Scaling
337 real(kind_real),intent(inout) :: zmin !< Minimum
338 real(kind_real),intent(inout) :: zmax !< Maximum
339 real(kind_real),intent(inout) :: zavg !< Average
340 
341 ! Local variables
342 integer :: jlev,jobs
343 real(kind_real) :: expo
344 
345 if (self%nobs>0.and.self%nlev>0) then
346  ! Compute statistics
347  if (.not.allocated(self%values)) call abor1_ftn('qg_obsvec_stats: obs vector not allocated')
348 
349  ! Initialization
350  zmin = huge(1.0)
351  zmax = -huge(1.0)
352  zavg = 0.0
353 
354  ! Loop over values
355  do jobs=1,self%nobs
356  do jlev=1,self%nlev
357  zmin = min(self%values(jlev,jobs),zmin)
358  zmax = max(self%values(jlev,jobs),zmax)
359  zavg = zavg+self%values(jlev,jobs)
360  enddo
361  enddo
362 
363  ! Normalization
364  zavg = zavg/real(self%nlev*self%nobs,kind_real)
365 else
366  ! Empty observation vector
367  zmin = 0.0
368  zmax = 0.0
369  zavg = 0.0
370 endif
371 
372 ! Scaling
373 if (abs(zavg)>0.0) then
374  expo = aint(log(abs(zavg))/log(10.0_kind_real))
375  scaling = 10.0**expo
376 else
377  scaling = 1.0
378 endif
379 zmin = zmin/scaling
380 zmax = zmax/scaling
381 zavg = zavg/scaling
382 
383 end subroutine qg_obsvec_stats
384 ! ------------------------------------------------------------------------------
385 !> Get observation vector size
386 subroutine qg_obsvec_nobs(self,kobs)
387 
388 implicit none
389 
390 ! Passed variables
391 type(qg_obsvec),intent(in) :: self !< Observation vector
392 integer,intent(inout) :: kobs !< Observation vector size
393 
394 ! Get observation vector size
395 kobs = self%nobs*self%nlev
396 
397 end subroutine qg_obsvec_nobs
398 
399 ! ------------------------------------------------------------------------------
400 !> Get value from observation vector at location (iob)
401 subroutine qg_obsvec_getat(self,iob,val)
402 
403 implicit none
404 
405 ! Passed variables
406 type(qg_obsvec),intent(in) :: self !< Observation vector
407 integer,intent(in) :: iob !< index into observation vector
408 real(kind_real), intent(out) :: val!< returned value
409 
410 integer :: i1, i2
411 
412 i1 = iob / self%nobs + 1
413 i2 = iob - self%nobs*(i1-1) + 1
414 ! Retrieve obs. value from vector
415 
416 if (i1>self%nlev .or. i2>self%nobs) call abor1_ftn ('qg_obsvec_getat: index is out of bounds')
417 
418 val = self%values(i1,i2)
419 
420 end subroutine qg_obsvec_getat
421 
422 ! ------------------------------------------------------------------------------
423 end module qg_obsvec_mod
qg_obsvec_mod::qg_obsvec_copy
subroutine, public qg_obsvec_copy(self, other)
Copy observation vector.
Definition: qg_obsvec_mod.F90:114
qg_obsvec_mod::qg_obsvec_mul_scal
subroutine, public qg_obsvec_mul_scal(self, zz)
Multiply observation vector with a scalar.
Definition: qg_obsvec_mod.F90:185
qg_obsvec_mod::qg_obsvec_dotprod
subroutine, public qg_obsvec_dotprod(obsvec1, obsvec2, zz)
Compute dot product between observation vectors.
Definition: qg_obsvec_mod.F90:303
qg_obsvec_mod::qg_obsvec_axpy
subroutine, public qg_obsvec_axpy(self, zz, other)
Apply axpy on observation vector.
Definition: qg_obsvec_mod.F90:255
qg_obsvec_mod::qg_obsvec_nobs
subroutine, public qg_obsvec_nobs(self, kobs)
Get observation vector size.
Definition: qg_obsvec_mod.F90:387
qg_obsvec_mod
Definition: qg_obsvec_mod.F90:10
qg_obsvec_mod::qg_obsvec_sub
subroutine, public qg_obsvec_sub(self, other)
Subtract observation vector.
Definition: qg_obsvec_mod.F90:213
qg_obsvec_mod::qg_obsvec_random_i
Definition: qg_obsvec_mod.F90:26
qg_obsvec_mod::qg_obsvec_zero
subroutine, public qg_obsvec_zero(self)
Set observation vector to zero.
Definition: qg_obsvec_mod.F90:172
qg_obsvec_mod::qg_obsvec_delete
subroutine, public qg_obsvec_delete(self)
Delete observation vector.
Definition: qg_obsvec_mod.F90:101
qg_obsvec_mod::qg_obsvec_clone
subroutine, public qg_obsvec_clone(self, other)
Clone observation vector.
Definition: qg_obsvec_mod.F90:83
qg_obsvec_mod::qg_obsvec_mul
subroutine, public qg_obsvec_mul(self, other)
Multiply observation vector.
Definition: qg_obsvec_mod.F90:227
qg_obsvec_mod::qg_obsvec_registry
type(registry_t), public qg_obsvec_registry
Linked list interface - defines registry_t type.
Definition: qg_obsvec_mod.F90:47
qg_obsvec_mod::qg_obsvec_random
subroutine, public qg_obsvec_random(c_odb, self)
Generate random observation vector.
Definition: qg_obsvec_mod.F90:283
qg_obsvec_mod::qg_obsvec_copy_local
subroutine, public qg_obsvec_copy_local(self, other, idx)
Copy a local subset of the observation vector.
Definition: qg_obsvec_mod.F90:140
qg_obsvec_mod::qg_obsvec_setup
subroutine, public qg_obsvec_setup(self, nlev, nobs)
Linked list implementation.
Definition: qg_obsvec_mod.F90:58
qg_obsvec_mod::qg_obsvec_div
subroutine, public qg_obsvec_div(self, other)
Divide observation vector.
Definition: qg_obsvec_mod.F90:241
qg_obsvec_mod::qg_obsvec_stats
subroutine, public qg_obsvec_stats(self, scaling, zmin, zmax, zavg)
Compute observation vector statistics.
Definition: qg_obsvec_mod.F90:331
qg_obsvec_mod::qg_obsvec_getat
subroutine, public qg_obsvec_getat(self, iob, val)
Get value from observation vector at location (iob)
Definition: qg_obsvec_mod.F90:402
qg_obsvec_mod::qg_obsvec
Definition: qg_obsvec_mod.F90:35
qg_obsvec_mod::qg_obsvec_add
subroutine, public qg_obsvec_add(self, other)
Add observation vector.
Definition: qg_obsvec_mod.F90:199
qg_obsvec_mod::qg_obsvec_invert
subroutine, public qg_obsvec_invert(self)
Invert observation vector.
Definition: qg_obsvec_mod.F90:270