OOPS
qg_error_covariance_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2009-2016 ECMWF.
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 ! In applying this licence, ECMWF does not waive the privileges and immunities
6 ! granted to it by virtue of its status as an intergovernmental organisation nor
7 ! does it submit to any jurisdiction.
8 
10 
11 use fckit_configuration_module, only: fckit_configuration
12 use fckit_log_module, only: fckit_log
13 use fft_mod
14 use iso_c_binding
15 use kinds
16 !$ use omp_lib
18 use qg_fields_mod
19 use qg_geom_mod
21 use random_mod
22 
23 implicit none
24 
25 private
30 ! ------------------------------------------------------------------------------
32  real(kind_real) :: sigma !< Standard deviation
33  real(kind_real),allocatable :: sqrt_zonal(:) !< Spectral weights for the spectral of the zonal correlation matrix
34  real(kind_real),allocatable :: sqrt_merid(:,:) !< Square-root of the meridional correlation matrix
35  real(kind_real),allocatable :: sqrt_vert(:,:) !< Square-root of the meridional correlation matrix
36  real(kind_real),allocatable :: norm(:,:) !< Normalization factor
38 
39 real(kind_real),parameter :: eps_ad = 1.0e-10 !< Epsilon value for adjoint tests
40 
41 #define LISTED_TYPE qg_error_covariance_config
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_error_covariance_registry
48 ! ------------------------------------------------------------------------------
49 contains
50 ! ------------------------------------------------------------------------------
51 ! Public
52 ! ------------------------------------------------------------------------------
53 !> Linked list implementation
54 #include "oops/util/linkedList_c.f"
55 ! ------------------------------------------------------------------------------
56 !> Setup error covariance matrix
57 subroutine qg_error_covariance_setup(self,f_conf,geom)
58 
59 implicit none
60 
61 ! Passed variables
62 type(qg_error_covariance_config),intent(inout) :: self !< Error covariance configuration
63 type(fckit_configuration),intent(in) :: f_conf !< FCKIT configuration
64 type(qg_geom),intent(in) :: geom !< Geometry
65 
66 ! Local variables
67 integer :: ix,iy,jy,ky,iz,jz,kz,info
68 real(kind_real) :: horizontal_length_scale,vertical_length_scale
69 real(kind_real) :: distx,condition_number,threshold
70 real(kind_real),allocatable :: struct_fn(:),workx(:)
71 real(kind_real),allocatable :: evalsy(:),worky(:),evectsy(:,:),revalsy(:)
72 real(kind_real),allocatable :: evalsz(:),workz(:),evectsz(:,:),revalsz(:)
73 real(kind_real),allocatable :: norm(:,:)
74 character(len=160) :: record
75 type(qg_fields) :: fld_in,fld_out
76 
77 ! Get parameters
78 call f_conf%get_or_die("standard_deviation",self%sigma)
79 call f_conf%get_or_die("horizontal_length_scale",horizontal_length_scale)
80 call f_conf%get_or_die("vertical_length_scale",vertical_length_scale)
81 call f_conf%get_or_die("maximum_condition_number",condition_number)
82 
83 ! Check nx
84 if (mod(geom%nx,2)/=0) then
85  write(record,*) 'qg_error_covariance_setup: number of zonal gridpoints nx=',geom%nx
86  call fckit_log%error(record)
87  call abor1_ftn('qg_error_covariance_setup: odd number of zonal grid points')
88 endif
89 
90 ! Allocation
91 allocate(self%sqrt_merid(geom%ny,geom%ny))
92 allocate(self%sqrt_vert(geom%nz,geom%nz))
93 allocate(self%norm(geom%ny,geom%nz))
94 allocate(struct_fn(geom%nx))
95 allocate(workx(geom%nx+2))
96 allocate(self%sqrt_zonal(0:geom%nx/2))
97 allocate(evectsy(geom%ny,geom%ny))
98 allocate(evalsy(geom%ny))
99 allocate(worky((geom%ny+3)*geom%ny))
100 allocate(revalsy(geom%ny))
101 allocate(evectsz(geom%nz,geom%nz))
102 allocate(evalsz(geom%nz))
103 allocate(workz((geom%nz+3)*geom%nz))
104 allocate(revalsz(geom%nz))
105 allocate(norm(geom%ny,geom%nz))
106 
107 ! Calculate spectral weights for zonal correlations:
108 ! First we construct the structure function in grid space and FFT it
109 do ix=1,geom%nx
110  distx = geom%x(ix)-geom%x(1)
111  if (distx>0.5*domain_zonal) distx = domain_zonal-distx
112  struct_fn(ix) = exp(-0.5*(distx/horizontal_length_scale)**2)
113 enddo
114 call fft_fwd(geom%nx,struct_fn,workx)
115 workx = workx/real(geom%nx,kind_real)
116 
117 ! We are after sqrt(B) or sqrt(Q),so square-root the coefficents.
118 ! The factor N ensures we get the struture function back if we apply
119 ! B or Q to a unit delta function.
120 threshold = 0.0
121 do ix=0,geom%nx/2
122  threshold = max(threshold,workx(1+2*ix))
123 enddo
124 threshold = threshold/condition_number
125 do ix=0,geom%nx/2
126  self%sqrt_zonal(ix) = sqrt(max(threshold,real(geom%nx,kind_real)*workx(1+2*ix)))
127 enddo
128 
129 ! Generate the lower triangle of the meridional correlation matrix
130 do jy=1,geom%ny
131  do iy=jy,geom%ny
132  evectsy(iy,jy) = exp(-0.5*((geom%y(iy)-geom%y(jy))/horizontal_length_scale)**2)
133  evectsy(jy,iy) = evectsy(iy,jy)
134  enddo
135 enddo
136 
137 ! Calculate the eigen-decomposition
138 call dsyev('V','L',geom%ny,evectsy,geom%ny,evalsy,worky,size(worky),info)
139 if (info/=0) then
140  write(record,*) 'qg_error_covariance_setup: dsyev returns info=',info
141  call fckit_log%error(record)
142  call abor1_ftn ('qg_error_covariance_setup: dsyev failed to sqrt matrix')
143 endif
144 
145 ! Calculate the lower triangle of the symmetric square-root
146 threshold = maxval(evalsy)/condition_number
147 do iy=1,geom%ny
148  evalsy(iy) = sqrt(max(threshold,evalsy(iy)))
149  revalsy(iy) = 1.0/evalsy(iy)
150 enddo
151 do jy=1,geom%ny
152  do iy=jy,geom%ny
153  self%sqrt_merid(iy,jy) = 0.0
154  do ky=1,geom%ny
155  self%sqrt_merid(iy,jy) = self%sqrt_merid(iy,jy)+evectsy(iy,ky)*evectsy(jy,ky)*evalsy(ky)
156  enddo
157  enddo
158 enddo
159 
160 ! Generate the lower triangle of the vertical correlation matrix
161 do jz=1,geom%nz
162  do iz=jz,geom%nz
163  evectsz(iz,jz) = exp(-0.5*((geom%z(iz)-geom%z(jz))/vertical_length_scale)**2)
164  evectsz(jz,iz) = evectsz(iz,jz)
165  enddo
166 enddo
167 
168 ! Calculate the eigen-decomposition
169 call dsyev('V','L',geom%nz,evectsz,geom%nz,evalsz,workz,size(workz),info)
170 if (info/=0) then
171  write(record,*) 'qg_error_covariance_setup: dsyev returns info=',info
172  call fckit_log%error(record)
173  call abor1_ftn ('qg_error_covariance_setup: dsyev failed to sqrt matrix')
174 endif
175 
176 ! Calculate the lower triangle of the symmetric square-root
177 threshold = maxval(evalsz)/condition_number
178 do iz=1,geom%nz
179  evalsz(iz) = sqrt(max(threshold,evalsz(iz)))
180  revalsz(iz) = 1.0/evalsz(iz)
181 enddo
182 do jz=1,geom%nz
183  do iz=jz,geom%nz
184  self%sqrt_vert(iz,jz) = 0.0
185  do kz=1,geom%nz
186  self%sqrt_vert(iz,jz) = self%sqrt_vert(iz,jz)+evectsz(iz,kz)*evectsz(jz,kz)*evalsz(kz)
187  enddo
188  enddo
189 enddo
190 
191 ! Compute normalization factor
192 call qg_fields_create_default(fld_in,geom,.false.)
193 call qg_fields_create_default(fld_out,geom,.false.)
194 self%norm = 1.0
195 do iz=1,geom%nz
196  do iy=1,geom%ny
197  call qg_fields_zero(fld_in)
198  fld_in%gfld3d(1,iy,iz) = 1.0
199  call qg_error_covariance_mult(self,fld_in,fld_out)
200  norm(iy,iz) = 1.0/sqrt(fld_out%gfld3d(1,iy,iz))
201  end do
202 end do
203 self%norm = norm*self%sigma
204 
205 ! Release memory
206 deallocate(struct_fn)
207 deallocate(workx)
208 deallocate(worky)
209 deallocate(revalsy)
210 deallocate(evectsy)
211 deallocate(evalsy)
212 deallocate(workz)
213 deallocate(revalsz)
214 deallocate(evectsz)
215 deallocate(evalsz)
216 call qg_fields_delete(fld_in)
217 call qg_fields_delete(fld_out)
218 
219 end subroutine qg_error_covariance_setup
220 ! ------------------------------------------------------------------------------
221 !> Delete error covariance matrix
223 
224 implicit none
225 
226 ! Passed variables
227 type(qg_error_covariance_config),intent(inout) :: self !< Error covariance configuration
228 
229 ! Release memory
230 deallocate(self%sqrt_zonal)
231 deallocate(self%sqrt_merid)
232 deallocate(self%sqrt_vert)
233 deallocate(self%norm)
234 
235 end subroutine qg_error_covariance_delete
236 ! ------------------------------------------------------------------------------
237 !> Multiply by error covariance matrix
238 subroutine qg_error_covariance_mult(conf,fld_in,fld_out)
239 
240 implicit none
241 
242 ! Passed variables
243 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
244 type(qg_fields),intent(in) :: fld_in !< Input field
245 type(qg_fields),intent(inout) :: fld_out !< Output field
246 
247 ! Local variables
248 type(qg_fields) :: fld_tmp
249 
250 ! Initialization
251 call qg_fields_create_from_other(fld_tmp,fld_in)
252 call qg_fields_copy(fld_tmp,fld_in)
253 
254 ! Apply covariance matrix
255 call qg_error_covariance_sqrt_mult_ad(conf,fld_tmp,fld_out)
256 call qg_fields_copy(fld_tmp,fld_out)
257 call qg_error_covariance_sqrt_mult(conf,fld_tmp,fld_out)
258 
259 end subroutine qg_error_covariance_mult
260 ! ------------------------------------------------------------------------------
261 !> Randomize error covariance
262 subroutine qg_error_covariance_randomize(conf,fld_out)
263 
264 implicit none
265 
266 ! Passed variables
267 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
268 type(qg_fields),intent(inout) :: fld_out !< Output field
269 
270 ! Local variables
271 type(qg_fields) :: fld_tmp
272 
273 ! Initialize temporary field
274 call qg_fields_create_from_other(fld_tmp,fld_out)
275 call qg_fields_random(fld_tmp)
276 
277 ! Apply square-root of the covariance matrix
278 call qg_error_covariance_sqrt_mult(conf,fld_tmp,fld_out)
279 
280 end subroutine qg_error_covariance_randomize
281 ! ------------------------------------------------------------------------------
282 ! Private
283 ! ------------------------------------------------------------------------------
284 !> Multiply by error covariance matrix square-root, zonal part
285 subroutine qg_error_covariance_sqrt_mult_zonal(conf,fld_in,fld_out)
286 
287 implicit none
288 
289 ! Passed variables
290 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
291 type(qg_fields),intent(in) :: fld_in !< Input field
292 type(qg_fields),intent(inout) :: fld_out !< Output field
293 
294 ! Local variables
295 integer :: iy,iz,m,iri
296 real(kind_real) :: zfour(fld_in%geom%nx+2)
297 
298 do iz=1,fld_in%geom%nz
299  do iy=1,fld_in%geom%ny
300  call fft_fwd(fld_in%geom%nx,fld_in%gfld3d(:,iy,iz),zfour)
301  !$omp parallel do schedule(static) private(m,iri)
302  do m=0,fld_in%geom%nx/2
303  do iri=1,2
304  zfour(2*m+iri) = zfour(2*m+iri)*conf%sqrt_zonal(m)
305  enddo
306  enddo
307  !$omp end parallel do
308  call fft_inv(fld_in%geom%nx,zfour,fld_out%gfld3d(:,iy,iz))
309  enddo
310 enddo
311 
313 ! ------------------------------------------------------------------------------
314 !> Multiply by error covariance matrix square-root - meridional part
315 subroutine qg_error_covariance_sqrt_mult_meridional(conf,fld_in,fld_out)
316 
317 implicit none
318 
319 ! Passed variables
320 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
321 type(qg_fields),intent(in) :: fld_in !< Input field
322 type(qg_fields),intent(inout) :: fld_out !< Output field
323 
324 ! Local variables
325 integer :: ix,iz
326 real(kind_real),allocatable :: arr_in(:),arr_out(:)
327 
328 !$omp parallel do schedule(static) private(iz,ix) firstprivate(arr_in,arr_out)
329 do iz=1,fld_in%geom%nz
330  do ix=1,fld_in%geom%nx
331  ! Allocation
332  allocate(arr_in(fld_in%geom%ny))
333  allocate(arr_out(fld_in%geom%ny))
334 
335  ! Initialize
336  arr_in = fld_in%gfld3d(ix,:,iz)
337 
338  ! Apply transform
339  call dsymv('L',fld_in%geom%ny,1.0_kind_real,conf%sqrt_merid,fld_in%geom%ny,arr_in,1,0.0_kind_real,arr_out,1)
340 
341  ! Copy
342  fld_out%gfld3d(ix,:,iz) = arr_out
343 
344  ! Release memory
345  deallocate(arr_in)
346  deallocate(arr_out)
347  enddo
348 enddo
349 !$omp end parallel do
350 
352 ! ------------------------------------------------------------------------------
353 !> Multiply by error covariance matrix square-root - vertical part
354 subroutine qg_error_covariance_sqrt_mult_vertical(conf,fld_in,fld_out)
355 
356 implicit none
357 
358 ! Passed variables
359 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
360 type(qg_fields),intent(in) :: fld_in !< Input field
361 type(qg_fields),intent(inout) :: fld_out !< Output field
362 
363 ! Local variables
364 integer :: ix,iy
365 real(kind_real),allocatable :: arr_in(:),arr_out(:)
366 
367 !$omp parallel do schedule(static) private(iy,ix) firstprivate(arr_in,arr_out)
368 do iy=1,fld_in%geom%ny
369  do ix=1,fld_in%geom%nx
370  ! Allocation
371  allocate(arr_in(fld_in%geom%nz))
372  allocate(arr_out(fld_in%geom%nz))
373 
374  ! Initialize
375  arr_in = fld_in%gfld3d(ix,iy,:)
376 
377  ! Apply transform
378  call dsymv('L',fld_in%geom%nz,1.0_kind_real,conf%sqrt_vert,fld_in%geom%nz,arr_in,1,0.0_kind_real,arr_out,1)
379 
380  ! Copy
381  fld_out%gfld3d(ix,iy,:) = arr_out
382 
383  ! Release memory
384  deallocate(arr_in)
385  deallocate(arr_out)
386  enddo
387 enddo
388 !$omp end parallel do
389 
391 ! ------------------------------------------------------------------------------
392 !> Multiply by error covariance matrix square-root
393 subroutine qg_error_covariance_sqrt_mult(conf,fld_in,fld_out)
394 
395 implicit none
396 
397 ! Passed variables
398 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
399 type(qg_fields),intent(in) :: fld_in !< Input field
400 type(qg_fields),intent(inout) :: fld_out !< Output field
401 
402 ! Local variables
403 integer :: ix
404 type(qg_fields) :: fld_tmp
405 
406 ! Initialization
407 call qg_fields_create_from_other(fld_tmp,fld_in)
408 call qg_fields_copy(fld_tmp,fld_in)
409 
410 ! Multiply by symmetric square-root of vertical correlation matrix
411 call qg_error_covariance_sqrt_mult_vertical(conf,fld_tmp,fld_out)
412 call qg_fields_copy(fld_tmp,fld_out)
413 
414 ! Multiply by square-root of meridional correlation matrix
415 call qg_error_covariance_sqrt_mult_meridional(conf,fld_tmp,fld_out)
416 call qg_fields_copy(fld_tmp,fld_out)
417 
418 ! Multiply by square-root of zonal correlation matrix
419 call qg_error_covariance_sqrt_mult_zonal(conf,fld_tmp,fld_out)
420 call qg_fields_copy(fld_tmp,fld_out)
421 
422 ! Multiply by normalization factor
423 !$omp parallel do schedule(static) private(ix)
424 do ix=1,fld_in%geom%nx
425  fld_out%gfld3d(ix,:,:) = fld_tmp%gfld3d(ix,:,:)*conf%norm
426 end do
427 !$omp end parallel do
428 call qg_fields_copy(fld_tmp,fld_out)
429 
430 ! Multiply by standard deviation
431 fld_out%gfld3d = fld_tmp%gfld3d*conf%sigma
432 
433 end subroutine qg_error_covariance_sqrt_mult
434 ! ------------------------------------------------------------------------------
435 !> Multiply by error covariance matrix square-root - adjoint
436 subroutine qg_error_covariance_sqrt_mult_ad(conf,fld_in,fld_out)
437 
438 implicit none
439 
440 ! Passed variables
441 type(qg_error_covariance_config),intent(in) :: conf !< Error covariance configuration
442 type(qg_fields),intent(in) :: fld_in !< Input field
443 type(qg_fields),intent(inout) :: fld_out !< Output field
444 
445 ! Local variables
446 integer :: ix
447 type(qg_fields) :: fld_tmp
448 
449 ! Initialization
450 call qg_fields_create_from_other(fld_tmp,fld_in)
451 call qg_fields_copy(fld_tmp,fld_in)
452 
453 ! Multiply by standard deviation
454 fld_out%gfld3d = fld_tmp%gfld3d*conf%sigma
455 call qg_fields_copy(fld_tmp,fld_out)
456 
457 ! Multiply by normalization factor
458 !$omp parallel do schedule(static) private(ix)
459 do ix=1,fld_in%geom%nx
460  fld_out%gfld3d(ix,:,:) = fld_tmp%gfld3d(ix,:,:)*conf%norm
461 end do
462 !$omp end parallel do
463 call qg_fields_copy(fld_tmp,fld_out)
464 
465 ! Multiply by square-root of zonal correlation matrix
466 call qg_error_covariance_sqrt_mult_zonal(conf,fld_tmp,fld_out)
467 call qg_fields_copy(fld_tmp,fld_out)
468 
469 ! Multiply by square-root of meridional correlation matrix
470 call qg_error_covariance_sqrt_mult_meridional(conf,fld_tmp,fld_out)
471 call qg_fields_copy(fld_tmp,fld_out)
472 
473 ! Multiply by symmetric square-root of vertical correlation matrix
474 call qg_error_covariance_sqrt_mult_vertical(conf,fld_tmp,fld_out)
475 
477 ! ------------------------------------------------------------------------------
478 end module qg_error_covariance_mod
qg_fields_mod
Definition: qg_fields_mod.F90:9
qg_error_covariance_mod
Definition: qg_error_covariance_mod.F90:9
qg_fields_mod::qg_fields
Definition: qg_fields_mod.F90:51
qg_geom_mod
Definition: qg_geom_mod.F90:9
qg_constants_mod::domain_zonal
real(kind_real), parameter, public domain_zonal
Model domain in zonal direction (m)
Definition: qg_constants_mod.F90:32
qg_error_covariance_mod::qg_error_covariance_setup
subroutine, public qg_error_covariance_setup(self, f_conf, geom)
Linked list implementation.
Definition: qg_error_covariance_mod.F90:58
qg_error_covariance_mod::qg_error_covariance_registry
type(registry_t), public qg_error_covariance_registry
Linked list interface - defines registry_t type.
Definition: qg_error_covariance_mod.F90:47
fft_mod
Fortran module for Eigen FFTs.
Definition: fft_mod.F90:7
qg_error_covariance_mod::qg_error_covariance_config
Definition: qg_error_covariance_mod.F90:31
qg_fields_mod::qg_fields_zero
subroutine, public qg_fields_zero(self)
Set fields to zero.
Definition: qg_fields_mod.F90:218
qg_error_covariance_mod::qg_error_covariance_sqrt_mult_zonal
subroutine qg_error_covariance_sqrt_mult_zonal(conf, fld_in, fld_out)
Multiply by error covariance matrix square-root, zonal part.
Definition: qg_error_covariance_mod.F90:286
qg_error_covariance_mod::qg_error_covariance_delete
subroutine, public qg_error_covariance_delete(self)
Delete error covariance matrix.
Definition: qg_error_covariance_mod.F90:223
qg_geom_mod::qg_geom
Definition: qg_geom_mod.F90:26
oops_variables_mod
Fortran interface to Variables.
Definition: variables_mod.F90:9
qg_error_covariance_mod::qg_error_covariance_sqrt_mult_ad
subroutine qg_error_covariance_sqrt_mult_ad(conf, fld_in, fld_out)
Multiply by error covariance matrix square-root - adjoint.
Definition: qg_error_covariance_mod.F90:437
qg_constants_mod
Definition: qg_constants_mod.F90:9
qg_error_covariance_mod::eps_ad
real(kind_real), parameter eps_ad
Epsilon value for adjoint tests.
Definition: qg_error_covariance_mod.F90:39
qg_error_covariance_mod::qg_error_covariance_randomize
subroutine, public qg_error_covariance_randomize(conf, fld_out)
Randomize error covariance.
Definition: qg_error_covariance_mod.F90:263
qg_error_covariance_mod::qg_error_covariance_sqrt_mult_meridional
subroutine qg_error_covariance_sqrt_mult_meridional(conf, fld_in, fld_out)
Multiply by error covariance matrix square-root - meridional part.
Definition: qg_error_covariance_mod.F90:316
fft_mod::fft_fwd
subroutine, public fft_fwd(kk, pg, pf)
Definition: fft_mod.F90:49
qg_error_covariance_mod::qg_error_covariance_sqrt_mult
subroutine qg_error_covariance_sqrt_mult(conf, fld_in, fld_out)
Multiply by error covariance matrix square-root.
Definition: qg_error_covariance_mod.F90:394
qg_fields_mod::qg_fields_create_from_other
subroutine, public qg_fields_create_from_other(self, other)
Create fields from another one.
Definition: qg_fields_mod.F90:168
qg_fields_mod::qg_fields_random
subroutine, public qg_fields_random(self)
Generate random fields.
Definition: qg_fields_mod.F90:305
fft_mod::fft_inv
subroutine, public fft_inv(kk, pf, pg)
Definition: fft_mod.F90:66
qg_error_covariance_mod::qg_error_covariance_sqrt_mult_vertical
subroutine qg_error_covariance_sqrt_mult_vertical(conf, fld_in, fld_out)
Multiply by error covariance matrix square-root - vertical part.
Definition: qg_error_covariance_mod.F90:355
qg_fields_mod::qg_fields_copy
subroutine, public qg_fields_copy(self, other, bconly)
Copy fields.
Definition: qg_fields_mod.F90:321
qg_fields_mod::qg_fields_delete
subroutine, public qg_fields_delete(self)
Delete fields.
Definition: qg_fields_mod.F90:201
qg_fields_mod::qg_fields_create_default
subroutine, public qg_fields_create_default(self, geom, lbc)
Create fields from geometry (x)
Definition: qg_fields_mod.F90:129
qg_error_covariance_mod::qg_error_covariance_mult
subroutine, public qg_error_covariance_mult(conf, fld_in, fld_out)
Multiply by error covariance matrix.
Definition: qg_error_covariance_mod.F90:239