gen_be_ensrf.f90
References to this file elsewhere.
1 program gen_be_ensrf
2 !
3 !----------------------------------------------------------------------
4 ! Purpose : Perform an Ensemble Square Root Filter (EnSRF) assimilation
5 ! using WRF forecast data.
6 !
7 ! Owner: Dale Barker (NCAR/MMM)
8 ! Please acknowledge author/institute in work that uses this code.
9 !
10 !----------------------------------------------------------------------
11
12 #ifdef crayx1
13 #define iargc ipxfargc
14 #endif
15
16 use da_control, only : stdout,stderr,filename_len
17 use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh
18
19 implicit none
20
21 type info_type
22 character*40 :: id
23 character*40 :: platform
24 integer :: plat
25 integer :: variable
26 end type info_type
27
28 type loc_type
29 real :: x
30 real :: y
31 real :: z
32 end type loc_type
33
34 type ob_type
35 type (info_type):: info
36 type (loc_type) :: loc
37 real :: yo
38 real :: omb_mean
39 real :: oma_mean
40 real :: omb_rms
41 real :: oma_rms
42 real :: sigma_o
43 real :: sigma_b
44 end type ob_type
45
46 integer, parameter :: nv3d = 4 ! #3D variables (u, v, T, q).
47 integer, parameter :: nv2d = 1 ! #2D variables (ps).
48 integer, parameter :: unit = 100 ! Unit number.
49
50 character (len=filename_len) :: filestub ! General filename stub.
51 character (len=filename_len) :: input_file ! Input file.
52 character (len=filename_len) :: output_file ! Output file.
53 character (len=10) :: var ! Variable to search for.
54 character (len=3) :: ce ! Member index -> character.
55
56 integer :: num_members ! Ensemble size.
57 integer :: num_obs ! Number of observations.
58 integer :: ni ! 1st dimension size.
59 integer :: niu ! 1st dimension size (u)
60 integer :: nj ! 2nd dimension size.
61 integer :: njv ! 2nd dimension size (v)
62 integer :: nk ! 3rd dimension size.
63 integer :: nv ! Total number of 2D/3D variables.
64 integer :: nij, nijk, nijkv ! Dimensions.
65 integer :: o, member, i, j, k, v ! Loop counters.
66 integer :: imin, imax, jmin, jmax ! Min/max i,j for localization.
67 integer :: ijkv !
68 integer :: index !
69 integer :: xend
70 real :: num_members_inv
71 real :: num_members_inv1
72 real :: cov_inf_fac ! Covariance inflation factor.
73 real :: cov_loc_rad_m ! Covariance localization radius (m).
74 real :: cov_loc_rad ! Covariance localization radius (gridpts).
75 real :: ds ! Grid resolution (m).
76 real :: x1, y1 ! Grid location of ob.
77 real :: sigma_o2, sigma_b2
78 real :: sigma_hh_inv
79 real :: sigma_xh
80 real :: beta
81 real :: o_minus_b ! Observation minus forecast.
82 real :: kalman_gain
83 real :: x_mean
84 real :: h_mean
85 real :: hb_mean
86 real :: ha_mean
87 real :: sum_f, sum_a
88 real :: hf, ha
89
90 type (ob_type),pointer:: ob(:)
91
92 character(len=10),allocatable :: varr(:) ! Variable name.
93 integer, allocatable :: nkk(:) ! Vertical dimension of field.
94 integer, allocatable :: xstart(:) ! Starting position of variable..
95
96 real, allocatable :: x(:,:), xf(:,:)
97 ! real, allocatable :: xf_cif(:,:)
98 real, allocatable :: h(:), h_devn(:)
99 real, allocatable :: x_devn(:)
100
101 real, allocatable :: uc(:,:) ! u-wind (C grid).
102 real, allocatable :: vc(:,:) ! v-wind (C grid).
103 real, allocatable :: field(:,:) ! Input field.
104 real, allocatable :: dummy(:,:) ! Dummy.
105 real, allocatable :: corr(:,:) ! Correlation.
106 real, allocatable :: xf_mean(:) ! Prior mean.
107 real, allocatable :: xa_mean(:) ! Posterior mean.
108
109 namelist / gen_be_ensrf_nl / filestub, num_members, &
110 cov_inf_fac, cov_loc_rad_m
111
112 stderr = 0
113 stdout = 6
114
115 !---------------------------------------------------------------------------------------------
116 write(6,'(/a)')' [1] Initialize information.'
117 !---------------------------------------------------------------------------------------------
118
119 filestub = 'test'
120 num_members = 56
121 cov_inf_fac = 1.0
122 cov_loc_rad_m = 1500000.0
123
124 open(unit=unit, file='gen_be_ensrf_nl.nl', &
125 form='formatted', status='old', action='read')
126 read(unit, gen_be_ensrf_nl)
127 close(unit)
128
129 write(6,'(a,a)')' Filestub = ', trim(filestub)
130 write(6,'(a,i4)')' Number of ensemble members = ', num_members
131 write(6,'(a,f16.8)')' Covariance Inflation Factor = ', cov_inf_fac
132 write(6,'(a,f16.8)')' Covariance Localization Radius (m) = ', cov_loc_rad_m
133
134 num_members_inv = 1.0 / real(num_members)
135 num_members_inv1 = 1.0 / real(num_members-1)
136
137 !---------------------------------------------------------------------------------------------
138 write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
139 !---------------------------------------------------------------------------------------------
140
141 ! Get grid dimensions from first T field:
142 input_file = trim(filestub)//'.e001'
143 var = 'T'
144 call da_stage0_initialize( input_file, var, ni, nj, nk, ds )
145 niu = ni+1 ! u i dimension is 1 larger.
146 njv = nj+1 ! v j dimension is 1 larger.
147 cov_loc_rad = cov_loc_rad_m / ds
148 nij = ni * nj
149 nijk = nij * nk
150 nijkv = nv3d * nijk + nv2d * nij ! #3D + #2D variables.
151 nv = nv3d + nv2d
152 print *, ni, nj, nk, nv, nij, nijk, nv3d, nv2d, nijkv
153 ! Allocate arrays:
154 allocate ( varr(1:nv) )
155 allocate ( nkk(1:nv) )
156 allocate ( xstart(1:nv) )
157 allocate( x(1:nijkv,1:num_members) )
158 allocate( xf(1:nijkv,1:num_members) )
159 ! allocate( xf_cif(1:nijkv,1:num_members) )
160 allocate( h(1:num_members) )
161 allocate( h_devn(1:num_members) )
162 allocate( x_devn(1:num_members) )
163 do v = 1, nv3d
164 nkk(v) = nk
165 end do
166 do v = nv3d + 1, nv
167 nkk(v) = 1
168 end do
169 allocate( field(1:ni,1:nj) )
170 allocate( uc(1:niu,1:nj) )
171 allocate( vc(1:ni,1:njv) )
172 allocate( dummy(1:ni,1:nj) )
173 allocate( corr(1:ni,1:nj) )
174
175 ! Hardwired:
176 varr(1) = 'U'
177 varr(2) = 'V'
178 varr(3) = 'T'
179 varr(4) = 'QVAPOR'
180 varr(5) = 'PSFC'
181
182 !---------------------------------------------------------------------------------------------
183 write(6,'(/a)')' [3] Read observations.'
184 !---------------------------------------------------------------------------------------------
185
186 input_file = 'observations'
187 open(unit, file = input_file, status='old')
188 !write(56,*)1
189 ! write(56,'(2a10,i6,5f16.8)')'test', 'dale', 3, &
190 ! 20.0, 20.0, 12.0, 250.0, 1.0
191 !stop
192 read(unit,*)num_obs
193 write(6,'(a,i10)')' Number of observations = ', num_obs
194 allocate( ob(1:num_obs) )
195
196 do o = 1, num_obs
197 read(unit,'(2a10,i6,5f16.8)')ob(o) % info % id, ob(o) % info % platform, &
198 ob(o) % info % variable, &
199 ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
200 ob(o) % yo, ob(o) % sigma_o
201 end do
202
203 !---------------------------------------------------------------------------------------------
204 write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
205 !---------------------------------------------------------------------------------------------
206
207 xstart(1) = 1
208 do v = 2, nv
209 xstart(v) = xstart(v-1) + ( ni * nj * nkk(v-1) )
210 end do
211
212 do member = 1, num_members
213
214 write(UNIT=ce,FMT='(i3.3)')member
215 input_file = trim(filestub)//'.e'//ce
216
217 do v = 1, nv
218 do k = 1, nkk(v)
219
220 if ( varr(v) == 'U' ) then
221 call da_get_field( input_file, varr(v), 3, niu, nj, nk, k, uc )
222 do j = 1, nj
223 do i = 1, ni
224 field(i,j) = 0.5 * ( uc(i,j) + uc(i+1,j) )
225 end do
226 end do
227 end if
228
229 if ( varr(v) == 'V' ) then
230 call da_get_field( input_file, varr(v), 3, ni, njv, nk, k, vc )
231 do j = 1, nj
232 do i = 1, ni
233 field(i,j) = 0.5 * ( vc(i,j) + vc(i,j+1) )
234 end do
235 end do
236 end if
237
238 if ( varr(v) == "T" ) then
239 ! Read theta, and convert to temperature:
240 call da_get_trh( input_file, ni, nj, nk, k, field, dummy )
241 end if
242
243 ! Read mixing ratio, and convert to specific humidity:
244 if ( varr(v) == 'QVAPOR' ) then
245 call da_get_field( input_file, varr(v), 3, ni, nj, nk, k, field )
246 field(:,:) = field(:,:) / ( 1.0 + field(:,:) )
247 end if
248
249 ! Read surface pressure:
250 if ( varr(v) == 'PSFC' ) then
251 call da_get_field( input_file, varr(v), 2, ni, nj, nk, k, field )
252 end if
253
254 ! Fill 4D array:
255 index = xstart(v) + (k-1) * nij
256 do j = 1, nj
257 do i = 1, ni
258 xf(index,member) = field(i,j)
259 index = index + 1
260 end do
261 end do
262 end do ! k
263 end do ! v
264 end do !member
265
266 ! Initialize analysis as first guess:
267 x(1:nijkv,1:num_members) = xf(1:nijkv,1:num_members)
268
269 ! Perform initial covariance inflation (to boost input be):
270 call da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
271 ! xf_cif(1:nijkv,1:num_members) = x(1:nijkv,1:num_members) ! Diagnostic
272
273 deallocate( uc )
274 deallocate( vc )
275 deallocate( field )
276 deallocate( dummy )
277
278 !---------------------------------------------------------------------------------------------
279 write(6,'(/a)')' [5] Perform EnSRF:'
280 !---------------------------------------------------------------------------------------------
281
282 do o = 1, num_obs
283 write(6,'(2(a,i8))')' Assimilating observation ', o, ' of ', num_obs
284 x1 = ob(o) % loc % x
285 y1 = ob(o) % loc % y
286 sigma_o2 = ob(o) % sigma_o * ob(o) % sigma_o
287
288 ! Perform observation operator:
289
290 do member = 1, num_members
291 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), h(member) )
292 end do
293
294 ! Mean and deviations of H(x):
295 h_mean = sum(h(1:num_members)) * num_members_inv
296 h_devn(1:num_members) = h(1:num_members) - h_mean
297
298 ! Covariance H Pf H^T:
299 sigma_b2 = sum(h_devn(1:num_members) * h_devn(1:num_members)) * num_members_inv1
300
301 ob(o) % sigma_b = sqrt(sigma_b2)
302
303 ! Store 1/(H Pf HT + R):
304 sigma_hh_inv = 1.0 / ( sigma_b2 + sigma_o2 )
305
306 ! Calculate EnSRF deviation Kalman Gain scaling:
307 beta = 1.0 / ( 1.0 + sqrt( sigma_o2 * sigma_hh_inv ) )
308
309 ! Ensemble mean innovation vector:
310 !!!DALE Uncomment for pseudo observation
311 ! print *, 'Running pseudo observation O-B=1.0:'
312 ! ob(o) % yo = h_mean + 1.0
313 !!!DALE
314 o_minus_b = ob(o) % yo - h_mean
315 write(6,'(a,2f15.5)')' Observation increment = ', o_minus_b
316 write(6,'(a,2f15.5)')' Observation/forecast error variance = ', sigma_o2, sigma_b2
317
318 !-----------------------------------------------------------------------------
319 ! Calculate local support:
320 !-----------------------------------------------------------------------------
321
322 imin = max(nint(ob(o) % loc % x - cov_loc_rad),1)
323 imax = min(nint(ob(o) % loc % x + cov_loc_rad),ni)
324 jmin = max(nint(ob(o) % loc % y - cov_loc_rad),1)
325 jmax = min(nint(ob(o) % loc % y + cov_loc_rad),nj)
326
327 ! write(6,'(a,2f8.2)')' Ob location(x/y) = ', ob(o) % loc % x, ob(o) % loc % y
328 ! write(6,'(a,4i5)')' Min/max i/j = ', imin, imax, jmin, jmax
329
330 do j = jmin, jmax
331 do i = imin, imax
332 call compact_support( x1, real(i), y1, real(j), cov_loc_rad, corr(i,j) )
333 end do
334 end do
335
336 do v = 1, nv
337 do k = 1, nkk(v)
338 do j = jmin, jmax
339 do i = imin, imax
340 ijkv = (v-1) * nijk + (k-1) * nij + (j-1) * ni + i
341
342 x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
343
344 x_devn(1:num_members) = x(ijkv,1:num_members) - x_mean
345
346 ! Calculate PfHT:
347 sigma_xh = sum( x_devn(1:num_members) * h_devn(1:num_members) ) * &
348 num_members_inv1
349
350 ! Apply covariance localization:
351 sigma_xh = sigma_xh * corr(i,j)
352
353 kalman_gain = sigma_xh * sigma_hh_inv
354
355 x_mean = x_mean + kalman_gain * o_minus_b
356
357 x_devn(1:num_members) = x_devn(1:num_members) - &
358 beta * kalman_gain * h_devn(1:num_members)
359
360 x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
361
362 end do
363 end do
364 end do
365 end do
366 end do
367
368 !---------------------------------------------------------------------------------------------
369 write(6,'(/a)')' [4] Calculate diagnostics:'
370 !---------------------------------------------------------------------------------------------
371
372 ! Calculate mean analysis over all grid-points (diagnostic only):
373
374 allocate( xf_mean(1:nijkv) )
375 allocate( xa_mean(1:nijkv) )
376
377 do ijkv = 1, nijkv
378 xf_mean(ijkv) = sum(xf(ijkv,1:num_members)) * num_members_inv
379 xa_mean(ijkv) = sum(x(ijkv,1:num_members)) * num_members_inv
380 end do
381
382 ! Calculate grid-point diagnostics:
383
384 ! call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf_cif, x )
385 call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, x )
386
387 ! Calculate observation diagnostics:
388
389 write(53,*)num_obs
390 do o = 1, num_obs
391
392 ! Ensemble mean-observation diagnostics:
393 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_mean, hb_mean )
394 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xa_mean, ha_mean )
395
396 ob(o) % omb_mean = ob(o) % yo - hb_mean
397 ob(o) % oma_mean = ob(o) % yo - ha_mean
398
399 ! Ensemble spread-observation diagnostics:
400 sum_f = 0.0
401 sum_a = 0.0
402 do member = 1, num_members
403 ! call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_cif(:,member), hf )
404 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf(:,member), hf )
405 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), ha )
406 sum_f = sum_f + ( ob(o) % yo - hf )**2
407 sum_a = sum_a + ( ob(o) % yo - ha )**2
408 end do
409
410 ob(o) % omb_rms = sqrt( sum_f * num_members_inv )
411 ob(o) % oma_rms = sqrt( sum_a * num_members_inv )
412
413 write(53,'(2a10,i8,10f16.8)')ob(o) % info % id, ob(o) % info % platform, &
414 ob(o) % info % variable, &
415 ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
416 ob(o) % yo, ob(o) % omb_mean, ob(o) % oma_mean, &
417 ob(o) % omb_rms, ob(o) % oma_rms, &
418 ob(o) % sigma_o, ob(o) % sigma_b
419 end do
420
421 !---------------------------------------------------------------------------------------------
422 write(6,'(/a)')' [4] Output EnSRF analysis ensemble:'
423 !---------------------------------------------------------------------------------------------
424
425 do member = 1, num_members
426 write(6,'(a,i4)')' Writing ensemble member increments for member', member
427 write(UNIT=ce,FMT='(i3.3)')member
428
429 do v = 1, nv
430
431 ! Output prior ensemble forecasts:
432 output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.e'//trim(ce)
433 open(unit, file = output_file, form='unformatted')
434 write(unit)ni, nj, nkk(v)
435 xend = xstart(v) + nij * nkk(v) - 1
436 write(unit)xf(xstart(v):xend,member)
437 close(unit)
438
439 ! Output posterior ensemble forecasts:
440 output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.e'//trim(ce)
441 open(unit, file = output_file, form='unformatted')
442 write(unit)ni, nj, nkk(v)
443 xend = xstart(v) + nij * nkk(v) - 1
444 write(unit)x(xstart(v):xend,member)
445 close(unit)
446 end do
447 end do
448
449 do v = 1, nv
450
451 ! Output prior ensemble forecasts:
452 output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.mean'
453 open(unit, file = output_file, form='unformatted')
454 write(unit)ni, nj, nkk(v)
455 xend = xstart(v) + nij * nkk(v) - 1
456 write(unit)xf_mean(xstart(v):xend)
457 close(unit)
458
459 ! Output posterior ensemble forecasts:
460 output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.mean'
461 open(unit, file = output_file, form='unformatted')
462 write(unit)ni, nj, nkk(v)
463 xend = xstart(v) + nij * nkk(v) - 1
464 write(unit)xa_mean(xstart(v):xend)
465 close(unit)
466 end do
467
468 deallocate( ob )
469
470 contains
471
472 !-----------------------------------------------------------------------------
473
474 subroutine da_obs_operator( nv, ni, nj, nkk, xstart, ob, x, h )
475
476 implicit none
477
478 integer, intent(in) :: nv, ni, nj
479 integer, intent(in) :: nkk(1:nv)
480 integer, intent(in) :: xstart(1:nv)
481 type (ob_type), intent(in) :: ob
482 ! WHY? was this a target?
483 ! real, intent(in), target :: x(:)
484 real, intent(in) :: x(:)
485 real, intent(out) :: h
486
487 integer :: v, nk
488 integer :: start
489 real :: ob_x, ob_y, ob_z ! Observation location.
490
491 v = ob % info % variable
492 nk = nkk(v)
493 start = xstart( v )
494
495 ! Perform 3D interpolation:
496
497 ob_x = ob % loc % x
498 ob_y = ob % loc % y
499 ob_z = ob % loc % z
500
501 call da_interpolate_3d( ni, nj, nk, ob % loc % x, ob % loc % y, ob % loc % z, &
502 x(start:start + ni * nj * nk ), h )
503
504 end subroutine da_obs_operator
505
506 subroutine da_interpolate_3d( ni, nj, nk, x, y, z, field, h )
507
508 implicit none
509 integer, intent(in) :: ni, nj, nk
510 real, intent(in) :: x, y, z ! Grid location of point.
511 real, intent(in) :: field(1:ni,1:nj,1:nk) ! Field to interpolate.
512 real, intent(out) :: h ! Interpolated value.
513
514 integer :: ii, jj, kk ! Grid locators.
515 real :: dx, dy, dz ! Interpolation weights.
516 real :: h1, h2 ! Interpolation values.
517 real :: h_below ! 2D interp values.
518 real :: h_above ! 2D interp values.
519
520 ii = int(x)
521 jj = int(y)
522 kk = int(z)
523
524 dx = x - real(ii)
525 dy = y - real(jj)
526 dz = z - real(kk)
527
528 ! 2D interpolation on layer below:
529 h1 = ( 1.0 - dx ) * field(ii,jj,kk) + dx * field(ii+1,jj,kk)
530 h2 = ( 1.0 - dx ) * field(ii,jj+1,kk) + dx * field(ii+1,jj+1,kk)
531 h_below = ( 1.0 - dy ) * h1 + dy * h2
532
533 ! 2D interpolation on layer above:
534 h1 = ( 1.0 - dx ) * field(ii,jj,kk+1) + dx * field(ii+1,jj,kk+1)
535 h2 = ( 1.0 - dx ) * field(ii,jj+1,kk+1) + dx * field(ii+1,jj+1,kk+1)
536 h_above = ( 1.0 - dy ) * h1 + dy * h2
537
538 ! Interpolation in vertical:
539 h = ( 1.0 - dz ) * h_below + dz * h_above
540
541 end subroutine da_interpolate_3d
542
543 subroutine compact_support( x1, x2, y1, y2, cov_loc_rad, corr )
544
545 ! Compact support according to (4.10) of Gaspari and Cohn (1999).
546
547 ! Assumes r >=0.
548
549 implicit none
550
551 real, intent(in) :: x1, x2 ! x values of two points.
552 real, intent(in) :: y1, y2 ! y values of two points.
553 real, intent(in) :: cov_loc_rad ! Cut-off correlation > 2c.
554 real, intent(out) :: corr ! Compactly-supported correlation.
555
556 real :: z ! Distance between points.
557 real :: r ! Ratio.
558 real :: r2, r3, r4, r5 ! Powers of r.
559
560 z = sqrt( ( x2 - x1 )**2 + ( y2 - y1)**2 )
561 r = z / cov_loc_rad
562
563 r2 = r * r
564 r3 = r2 * r
565 r4 = r3 * r
566 r5 = r4 * r
567
568 if ( r <= 1.0 ) then
569 corr = -0.25 * r5 + 0.5 * r4 + 0.625 * r3 - 5.0 * r2 / 3.0 + 1
570 else if ( r > 1.0 .and. r < 2.0 ) then
571 corr = r5 / 12.0 - 0.5 * r4 + 0.625 * r3 + 5.0 * r2 / 3.0 - &
572 5.0 * r + 4.0 - 2.0 / ( 3.0 * r )
573
574 else if ( r >= 2.0 ) then
575 corr = 0.0
576 end if
577
578 end subroutine compact_support
579
580 subroutine da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, xa )
581
582 implicit none
583
584 integer, intent(in) :: num_members
585 integer, intent(in) :: ni
586 integer, intent(in) :: nj
587 integer, intent(in) :: nk
588 integer, intent(in) :: nv
589 integer, intent(in) :: nijkv
590 integer, intent(in) :: nkk(1:nv)
591 character*10,intent(in):: varr(1:nv)
592 real, intent(in) :: xf(1:nijkv,1:num_members)
593 real, intent(in) :: xa(1:nijkv,1:num_members)
594
595 character*10 :: name
596 integer :: v, k, j, i, ijkv
597 real :: num_members_inv, nij_inv
598 real :: mnsq_f, mnsq_a
599 real :: domain_mean_f, domain_mean_a
600 real :: domain_stdv_f, domain_stdv_a
601 real :: mean_f(1:ni,1:nj), mean_a(1:ni,1:nj)
602 real :: stdv_f(1:ni,1:nj), stdv_a(1:ni,1:nj)
603
604 num_members_inv = 1.0 / num_members
605 nij_inv = 1.0 / real ( ni * nj )
606
607 write(6,'(a)')' Variable, Level, Mean_F, Mean_A, StDv_F, StdV_A'
608
609 ijkv = 0
610 do v = 1, nv
611 name = varr(v)
612 do k = 1, nkk(v)
613 do j = 1, nj
614 do i = 1, ni
615 ijkv = ijkv + 1
616 mean_f(i,j) = sum(xf(ijkv,1:num_members)) * num_members_inv
617 mean_a(i,j) = sum(xa(ijkv,1:num_members)) * num_members_inv
618 mnsq_f = sum(xf(ijkv,1:num_members)**2) * num_members_inv
619 mnsq_a = sum(xa(ijkv,1:num_members)**2) * num_members_inv
620 stdv_f(i,j) = mnsq_f - mean_f(i,j) * mean_f(i,j)
621 stdv_a(i,j) = mnsq_a - mean_a(i,j) * mean_a(i,j)
622 stdv_f(i,j) = 0.0
623 if ( stdv_f(i,j) > 0.0 ) stdv_f(i,j) = sqrt(stdv_f(i,j))
624 stdv_a(i,j) = 0.0
625 if ( stdv_a(i,j) > 0.0 ) stdv_a(i,j) = sqrt(stdv_a(i,j))
626 end do
627 end do
628
629 domain_mean_f = sum(mean_f(:,:)) * nij_inv
630 domain_mean_a = sum(mean_a(:,:)) * nij_inv
631 domain_stdv_f = sum(stdv_f(:,:)) * nij_inv
632 domain_stdv_a = sum(stdv_a(:,:)) * nij_inv
633
634 write(6,'(a,i6,4f16.8)')trim(name), k, domain_mean_f, domain_mean_a, &
635 domain_stdv_f, domain_stdv_a
636 end do
637 end do
638
639 end subroutine da_get_grid_stats
640
641 subroutine da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
642
643 implicit none
644
645 integer, intent(in) :: nijkv
646 integer, intent(in) :: num_members
647 real, intent(in) :: cov_inf_fac
648 real, intent(inout) :: x(1:nijkv,1:num_members)
649
650 integer :: ijkv
651 real :: num_members_inv
652 real :: x_devn(1:num_members)
653
654 num_members_inv = 1.0 / real(num_members)
655
656 do ijkv = 1, nijkv
657 x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
658 x_devn(1:num_members) = cov_inf_fac * ( x(ijkv,1:num_members) - x_mean )
659 x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
660 end do
661
662 end subroutine da_cov_inflation
663
664 end program gen_be_ensrf
665