11 use fckit_configuration_module,
only: fckit_configuration
12 use fckit_log_module,
only: fckit_log
32 real(kind_real) :: sigma
33 real(kind_real),
allocatable :: sqrt_zonal(:)
34 real(kind_real),
allocatable :: sqrt_merid(:,:)
35 real(kind_real),
allocatable :: sqrt_vert(:,:)
36 real(kind_real),
allocatable :: norm(:,:)
39 real(kind_real),
parameter ::
eps_ad = 1.0e-10
41 #define LISTED_TYPE qg_error_covariance_config
44 #include "oops/util/linkedList_i.f"
54 #include "oops/util/linkedList_c.f"
63 type(fckit_configuration),
intent(in) :: f_conf
64 type(
qg_geom),
intent(in) :: geom
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
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)
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')
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))
110 distx = geom%x(ix)-geom%x(1)
112 struct_fn(ix) = exp(-0.5*(distx/horizontal_length_scale)**2)
114 call fft_fwd(geom%nx,struct_fn,workx)
115 workx = workx/real(geom%nx,kind_real)
122 threshold = max(threshold,workx(1+2*ix))
124 threshold = threshold/condition_number
126 self%sqrt_zonal(ix) = sqrt(max(threshold,real(geom%nx,kind_real)*workx(1+2*ix)))
132 evectsy(iy,jy) = exp(-0.5*((geom%y(iy)-geom%y(jy))/horizontal_length_scale)**2)
133 evectsy(jy,iy) = evectsy(iy,jy)
138 call dsyev(
'V',
'L',geom%ny,evectsy,geom%ny,evalsy,worky,
size(worky),info)
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')
146 threshold = maxval(evalsy)/condition_number
148 evalsy(iy) = sqrt(max(threshold,evalsy(iy)))
149 revalsy(iy) = 1.0/evalsy(iy)
153 self%sqrt_merid(iy,jy) = 0.0
155 self%sqrt_merid(iy,jy) = self%sqrt_merid(iy,jy)+evectsy(iy,ky)*evectsy(jy,ky)*evalsy(ky)
163 evectsz(iz,jz) = exp(-0.5*((geom%z(iz)-geom%z(jz))/vertical_length_scale)**2)
164 evectsz(jz,iz) = evectsz(iz,jz)
169 call dsyev(
'V',
'L',geom%nz,evectsz,geom%nz,evalsz,workz,
size(workz),info)
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')
177 threshold = maxval(evalsz)/condition_number
179 evalsz(iz) = sqrt(max(threshold,evalsz(iz)))
180 revalsz(iz) = 1.0/evalsz(iz)
184 self%sqrt_vert(iz,jz) = 0.0
186 self%sqrt_vert(iz,jz) = self%sqrt_vert(iz,jz)+evectsz(iz,kz)*evectsz(jz,kz)*evalsz(kz)
198 fld_in%gfld3d(1,iy,iz) = 1.0
200 norm(iy,iz) = 1.0/sqrt(fld_out%gfld3d(1,iy,iz))
203 self%norm = norm*self%sigma
206 deallocate(struct_fn)
230 deallocate(self%sqrt_zonal)
231 deallocate(self%sqrt_merid)
232 deallocate(self%sqrt_vert)
233 deallocate(self%norm)
295 integer :: iy,iz,m,iri
296 real(kind_real) :: zfour(fld_in%geom%nx+2)
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)
302 do m=0,fld_in%geom%nx/2
304 zfour(2*m+iri) = zfour(2*m+iri)*conf%sqrt_zonal(m)
308 call fft_inv(fld_in%geom%nx,zfour,fld_out%gfld3d(:,iy,iz))
326 real(kind_real),
allocatable :: arr_in(:),arr_out(:)
329 do iz=1,fld_in%geom%nz
330 do ix=1,fld_in%geom%nx
332 allocate(arr_in(fld_in%geom%ny))
333 allocate(arr_out(fld_in%geom%ny))
336 arr_in = fld_in%gfld3d(ix,:,iz)
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)
342 fld_out%gfld3d(ix,:,iz) = arr_out
365 real(kind_real),
allocatable :: arr_in(:),arr_out(:)
368 do iy=1,fld_in%geom%ny
369 do ix=1,fld_in%geom%nx
371 allocate(arr_in(fld_in%geom%nz))
372 allocate(arr_out(fld_in%geom%nz))
375 arr_in = fld_in%gfld3d(ix,iy,:)
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)
381 fld_out%gfld3d(ix,iy,:) = arr_out
424 do ix=1,fld_in%geom%nx
425 fld_out%gfld3d(ix,:,:) = fld_tmp%gfld3d(ix,:,:)*conf%norm
431 fld_out%gfld3d = fld_tmp%gfld3d*conf%sigma
454 fld_out%gfld3d = fld_tmp%gfld3d*conf%sigma
459 do ix=1,fld_in%geom%nx
460 fld_out%gfld3d(ix,:,:) = fld_tmp%gfld3d(ix,:,:)*conf%norm